© S. Flaischlen et al. 2022

1. Introduction

The accretion rate is a crucial parameter for the understanding of accretion disks surrounding young stellar objects (YSOs) and planet formation (e.g., Ercolano & Pascucci 2017). The interaction between the circumstellar disk and the central star is thought to be regulated by magnetic field lines channeling the flow of material along accretion columns from the inner edge of the disk to the stellar surface (Hartmann et al. 2016). Accretion shocks produced when the infalling material hits the photosphere with nearly free-fall velocity produce a variety of observational features like optical and ultraviolet (UV) excess emission, spectral veiling, and strong optical emission lines. From these features, the accretion luminosity, L_acc, can be determined from spectroscopy (e.g., Alcalá et al. 2017) or photometry (e.g., Manara et al. 2012).

The accretion rate tends to decrease over time (e.g., Hartmann et al. 1998) as the disk material gets accreted or dispersed by disk winds, driven by irradiation from the star (Alexander et al. 2014), although several recent results show that accretion rates remain high in a number of older objects (e.g., Ingleby et al. 2014; Rugel et al. 2018; Venuti et al. 2019; Manara et al. 2020). Simulations suggest that X-ray emission is very efficient in photoevaporating the disk material (e.g., Ercolano et al. 2008a,b, 2021; Picogna et al. 2019) and may decrease the accretion rate (e.g., Drake et al. 2009).

Our statistical analysis of 332 accreting young stars in the Orion Nebula Cluster (ONC; Flaischlen et al. 2021, hereafter F21) suggests a weak anticorrelation between residual X-ray luminosities and accretion rates, as expected in the context of X-ray driven photoevaporation and in accordance with previous studies (Telleschi et al. 2007; Drake et al. 2009).

However, a possible caveat in the interpretation of the data stems from the temporal variability of YSOs over a wide range of wavelengths. The X-ray luminosity and the accretion rate are both subject to temporal variability (e.g., Wolk et al. 2004; Venuti et al. 2014, 2021). Since the X-ray data were obtained in January 2003 (Preibisch et al. 2005) and the Hubble Space Telescope (HST) data, on which the accretion rates are based, between October 2004 and April 2005 (Robberto et al. 2013), possible effects of the variability on timescales of a few years on the correlation analysis cannot be ruled out.

The accretion rate variability manifests itself mainly in the optical and near-UV bands (e.g., Venuti et al. 2014, 2015; Robinson & Espaillat 2019; Schneider et al. 2020). The observed brightness variations stem from a multitude of mechanisms, including accretion bursts (e.g., Venuti et al. 2014), inner disk warps (e.g., Frasca et al. 2020), accretion columns creating hot spots (e.g., Cody et al. 2014) moving in and out of view while the YSO rotates (e.g., Cody & Hillenbrand 2018; Schneider et al. 2020) and phases of stable or unstable accretion regimes (e.g., Sousa et al. 2016).

In recent years, an abundance of effort has been made regarding variability analyses of YSOs. They were monitored with ground-based facilities (e.g., Costigan et al. 2014; Venuti et al. 2014, 2021) and space-based instruments such as the Convection, Rotation and Planetary Transits (CoRoT) satellite (e.g., Venuti et al. 2017), the Kepler space telescope in its K2 mission (e.g., Pouilly et al. 2020; Rebull et al. 2020), the Transiting Exoplanet Survey Satellite (TESS; e.g., Thanathibodee et al. 2020), and the Microvariability & Oscillations of STars (MOST) space telescope (e.g., Sousa et al. 2016; Siwak et al. 2018).

The UV excess is a rather direct tracer of the accretion rate (e.g., Schneider et al. 2020). In this study, we are interested in characterizing the accretion rate variability through the UV excess of a particular sample of young stars in the ONC in order to analyze how this variability influences our statistical analysis of the relation between X-ray luminosities and accretion rates on timescales of weeks to a few years.

To this aim, we performed a multiyear and multicolor photometric monitoring of the ONC with the wide-angle camera (the Wendelstein Wide Field Imager) at the 2 m Fraunhofer Telescope located on Mount Wendelstein in Germany (Kosyra et al. 2014). The monitoring provided u′, g′, and r′ band photometry for several sources over the course of ∼4 yr, from which accretion rates could be estimated for a time span of ∼2 yr. The data were used to analyze how the accretion rate variability may influence the correlation analysis between X-ray activity and accretion and to characterize the accretion rate variability.

In Sect. 2, we describe the observation and data reduction steps of our WWFI data. Section 3 illustrates our method for deriving stellar parameters and accretion rates from the obtained photometry. In Sect. 4, we characterize the accretion rate variability and analyze how it influences the correlation analysis between X-ray luminosities and accretion rates.

2. Observation and data reduction

2.1. Target region

The ONC is among the nearest regions of ongoing star formation (Bally 2008), with a distance of ∼403 pc (Kuhn et al. 2019). Therefore, it serves as a benchmark for star formation and has been studied in detail over the years. As such, stellar parameters for a large fraction of the young stellar population are available in the literature (Manara et al. 2012). Especially relevant for this work are the spectral type determinations from Da Rio et al. (2012) and Hillenbrand et al. (2013). Together with our photometric observation detailed below, the spectral types allow an estimate of the accretion rates.

2.2. Data reduction and photometric calibration

We observed the ONC in a multiyear and multicolor monitoring with the 2 m Fraunhofer Telescope at the Wendelstein Observatory located on Mount Wendelstein in Germany (Kosyra et al. 2014). The telescope hosts the Wendelstein Wide Field Imager (WWFI), a wide field camera consisting of four 4096 px × 4109 px charge-coupled device (CCD) detectors arranged in a 2 × 2 mosaic pattern. With 27.6′×28.9′, the field of view covers the region of the X-ray data (17′×17′) used in our previous analysis of the relation between X-ray luminosities and accretion rates in F21, which was based on data from the Chandra Orion Ultradeep Project (COUP; Getman et al. 2005).

For the monitoring, the u′, g′, and r′ filters were used, which are approximately comparable with the photometric system of the Sloan Digital Sky Survey (Fukugita et al. 1996). Since the central wavelength of the u′ filter is close to the Balmer Jump, it is particularly sensitive for accretion effects. Figure 1 shows a WWFI composite image of the ONC taken with WWFI. The uncertainty of the u′ photometry in the center region is large due to the very bright nebulosity.

Fig. 1. WWFI u′g′r′ color composition of the ONC (image credit: Wendelstein Observatory). Overlaid are the stars monitored in all three bands. The yellow symbols indicate stars that have a match with the X-ray sources listed in the COUP catalog (J/ApJS/160/319/coup, Getman et al. 2005). The star symbols further indicate sources that were monitored for at least five epochs and allowed an estimate for the mass accretion rate.

The monitoring began on November 24, 2014, and comprises 42 epochs until February 17, 2019. In order to fill up the gaps of the CCD mosaic, a 13-step spiral dithering pattern of five exposures for each filter was utilized. With single exposure times of 30 s (u′), 20 s (g′) and 10 s (r′), sources with a brightness in the range of 11.6–20.8 mag in the u′ band and 12.1–21.0 mag (11.6–20.8 mag) in the g′ band (r′ band) could be detected with a sufficient signal to noise ratio and without saturation. The mean seeing amounts to 1.60″ ± 0.29″ for the u′ band, 1.39″ ± 0.28″ for the g′ band, and 1.19″ ± 0.29″ for the r′ band. Due to weather conditions and technical maintenance, not all bands could be observed at each epoch. Table A.1 in the appendix lists at which date how many sources were extracted for each filter.

The raw data were reduced using the WWFI pipeline (see Kluge et al. 2020 for a more detailed description of the typical reduction steps). The data reduction includes bias subtraction, bad pixel, cosmic ray and satellite masking as well as flat fielding. The instrumental magnitudes were converted to AB magnitudes and the zero-points were calibrated by comparing the magnitudes with Pan-STARRS (Flewelling et al. 2020) sources in the field of view for the g′ and r′ filters. In order to reduce the influence of variations due to different weather and instrumental conditions, the zero-points for each frame were further calibrated by aligning the magnitudes of non-variable sources with the corresponding magnitudes in a reference frame. Since there is no equivalent to the u′ filter in the Pan-STARRS catalog, we calibrated the zero-point by fitting the color locus of non-accreting stars in the g′−r′ versus u′−g′ color-color diagram to the line describing the photospheric colors determined via synthetic photometry on empirical spectra of diskless class III objects using the u′ passband from the WWFI setup. In Sect. 3.2.1, we explain in more detail how we obtained this line. Figure 2 shows the color-color diagram for all targets and epochs.

Fig. 2. g′−r′ vs. u′−g′ color-color diagram of YSOs monitored with WWFI as gray dots. The bigger dots indicate the objects displaying accretion according to Manara et al. (2012). The gray cross shows the typical uncertainty of the photometry. The thick line is the sequence tracing the colors of the photosphere, obtained by synthetic photometry on spectra of diskless class III objects observed with X-shooter (Manara et al. 2013, 2017b). The colors of these template spectra are shown as green squares. The arrow traces the reddening vector for AV = 1 mag, assuming the reddening law of Cardelli et al. (1989) and a galactic reddening parameter of RV = 3.1.

Finally, the obtained magnitudes and photometric uncertainties for the three filters for each target and epoch were compiled into a single catalog, where we have limited ourselves to entries with photometry available simultaneously in all three bands and an uncertainty in the u′ band magnitude < 0.1 mag. The catalog contains 169 stars with an average of ∼14.5 u′, g′, and r′ observations per star.

2.3. Data selection

In order to ensure that the catalog contains only photometry of accreting stars, the source positions were matched with the catalog of accreting stars from Manara et al. (2012), which are complete down to the hydrogen burning limit. This step reduced our sample size to 61. For 60 sources, we found parallaxes in the Gaia Early Data Release 3 (EDR3) catalog (Gaia Collaboration 2016, 2021). The mean parallax of the sources is $\langle \overline{\omega }\rangle =2.48\pm 0.06\mathrm{mas}$, which corresponds to a mean distance of ∼403 pc, in agreement with the distance of ${403}_{-6}^{+7}\mathrm{pc}$ determined by Kuhn et al. (2019). We excluded three sources as likely background stars according to the criteria that their distance is larger than the one reported by Kuhn et al. (2019) within 3σ and kept the one source without a listed parallax value (V1118 Ori) and assigned to it the distance of 403 pc, with the typical uncertainty of ≈ ± 9 pc deduced from our sample.

We complemented the catalog with the effective temperatures determined by Da Rio et al. (2012). Where available, we took more recent spectral types obtained by Hillenbrand et al. (2013) and converted them to T_eff according to the temperature scales of Kenyon & Hartmann (1995) and Luhman et al. (2003) for K and M type stars, respectively. During this process, we excluded one additional source due to its uncertain spectral classification. The final catalog contains data for 57 sources. It is arranged in two tables: the first (Table 1) lists the sources with the collected parameters from the literature. The second (Table 2) comprises the photometric data from our survey and the epochs as well as the object names according to Table 1.

3. Method

The accretion rates can be estimated from the observed u′, g′, and r′ photometry and the effective temperatures known from the literature with an approach similar to the one used by Da Rio et al. (2010) and Manara et al. (2012), where the latter called it the “2CD” method.

The method assumes that the position of the sources in the g′− r′ versus u′− g′ color-color diagram depends solely on three parameters: the effective temperature, the accretion luminosity, and the extinction. Without accretion, the colors would fall on a line determined by the photosphere of the stars at a position set by the effective temperature and scatter along the reddening vector. Accretion shifts the position towards bluer values to the colors of a purely accreting source.

The accretion luminosity and the extinction can be obtained simultaneously from the observed magnitudes and the T_eff values when the photospheric colors and the colors of the accretion source are known. In the following, we describe our implementation of the approach in more detail and explain how we obtained the required photospheric colors.

3.1. Flux density modeling

The flux density of the sources in question is assumed to be a composition of the pure photospheric flux density and the accretion flux density. Without extinction, the flux density of a source located at a distance d can be expressed through

$$\begin{array}{cc}& F_{\lambda }=I_{\mathrm{phot},\lambda }·\frac{\pi R_{\ast }^2}{d^2}+I_{\mathrm{acc},\lambda }·\frac{A_{\mathrm{acc}}}{d^2}\hfill \\ \hfill & ={\left(\frac{R_{\ast }}{R_{\odot }}\right)}^2·(F_{\mathrm{phot},\lambda }+\eta ·F_{\mathrm{acc},\lambda }),\hfill \end{array}$$(1)

where R_* is the radius of the source with the intensity I_{phot,  λ} and A_acc is the area of the accreting regions with the intensity I_acc,  λ. The factor (R_*/R_⊙)² scales the fluxes according to the radius, where it is assumed that F_{phot,  λ} and F_acc,  λ refer to an area similar to π $R_{\odot }^2$, and η ≔ A_acc/$(\pi R_{*}^2)$ describes the fraction of the projected area of the photosphere covered by the accreting regions.

3.2. Color modeling

Using synthetic photometry (see Casagrande & VandenBerg 2014 for a detailed description of the topic), the flux density can be related to the expected magnitudes m_ζ for the ζ filter (where ζ stands for u′, g′, and r′ or any other filter) via

$$\begin{array}{c}\hfill m_{\zeta }=-2.5log{\overline{F}}_{\zeta }-Z_{\zeta },\end{array}$$(2)

with the zero-point Z_ζ and

$$\begin{array}{c}\hfill {\overline{F}}_{\zeta }=\frac{\int \lambda F_{\lambda }{T}_{\zeta ,\lambda }\mathrm{d}\lambda }{\int \lambda T_{\zeta ,\lambda }\mathrm{d}\lambda },\end{array}$$(3)

where T_ζ, λ is the passthrough function of the respective filter. The magnitude according to Eq. (2) is given by

$$\begin{array}{c}\hfill m_{\zeta }(\eta )=-2.5log({\overline{F}}_{\mathrm{phot},\zeta }+\eta ·{\overline{F}}_{\mathrm{acc},\zeta })-5log\left(\frac{R_{\ast }}{R_{\odot }}\right)-Z_{\zeta }.\end{array}$$(4)

The respective colors are independent of R_*:

$$\begin{array}{c}\hfill m_{\zeta }(\eta )-m_{\mu }(\eta )=-2.5log\left(\frac{{\overline{F}}_{\mathrm{phot},\zeta }+\eta ·{\overline{F}}_{\mathrm{acc},\zeta }}{{\overline{F}}_{\mathrm{phot},\mu }+\eta ·{\overline{F}}_{\mathrm{acc},\mu }}\right)-Z_{\zeta }+Z_{\mu }.\end{array}$$(5)

In the following subsections, we explain how we obtained the photospheric fluxes described by ${\overline{F}}_{\mathrm{phot},\zeta }$ and the accretion flux ${\overline{F}}_{\mathrm{acc},\zeta }$.

3.2.1. Photospheric flux

The photospheric flux ${\overline{F}}_{\mathrm{phot},\zeta }$ required for Eq. (5) was determined from empirically obtained spectra of young, diskless class III objects serving as photospheric templates. They were observed with X-shooter by Manara et al. (2013, 2017a). We excluded the spectra of the targets named 2MASS J11195652-7504529, [LES2004] ChaI 601, [LES2004] ChaI 717, and V1251 Cen due to a low signal-to-noise ratio (S/N).

For each of the remaining template spectra, we calculated F_{phot,  λ} according to Eq. (1) from the observed flux density by scaling the distance of the template source, taken from the latest Gaia EDR3 release (Gaia Collaboration 2016, 2021), to the distance of the ONC. Furthermore, the radius of the template source was scaled to the solar radius. We obtained the radius of the template source from the luminosities determined by Manara et al. (2013, 2017a).

With F_{phot,  λ}, we calculated ${\overline{F}}_{\mathrm{phot},\zeta }$ with Eq. (3) using the passthrough functions T_ζ, λ of the u′, g′, and r′ filters used by the WWFI setup (Kosyra et al. 2014). In order to increase the S/N, the template spectra were median smoothed with a typical box-width of ∼ 3 Å. For convenience, we converted the fluxes to magnitudes using Eq. (2) and denoted them m_{phot,  ζ}.

Since the relation between these magnitudes and T_eff is complicated and to avoid a possible bias by assuming a specific shape, we utilized a nonparametric regression approach in order to approximate the relation between the magnitudes and the effective temperature. Similar to Manara et al. (2017a), we used a local second degree polynomial regression with a Gaussian kernel implemented in the Python toolkit pyqtfit.nprmethods.LocalPolynomialKernel for this task. The 2σ confidence level was obtained using bootstrapping.

In the final step, the fitting results were sampled at 20 T_eff values, from which the fitting results of the nonparametric regression can be recreated by interpolating with cubic splines. This allows ${\overline{F}}_{\mathrm{phot},\zeta }(T_{\mathrm{eff}})$ to be obtained for every given T_eff. A plot of the relation between m_{phot,  ζ} and T_eff as well as the fitting results are shown in Appendix D.

3.2.2. The accretion colors

In order to compare Eq. (5) with the observed data, the flux associated with the accretion ${\overline{F}}_{\mathrm{acc},\zeta }$ is required. We chose the models of Manara et al. (2013), which consist of the intensity of a slab of pure hydrogen (bound-free and free-free emission from H and H⁻) assuming local thermodynamic equilibrium. The models depend on three parameters: the temperature (T_slab), the optical depth at a reference wavelength of λ = 300 nm (τ₃₀₀), and the electron density (n_e).

We calculated the flux density F_acc,  λ according to Eq. (1) from the intensity I_acc,  λ of the slab and obtained the magnitudes in the WWFI filters for several model parameters using Eqs. (2) and (3). We denoted them m_acc,  ζ. In addition, we determined the luminosity L_acc,  0 of the model slab by integrating the flux density and setting a distance of d = 403 pc.

In order to select an appropriate accretion model for our data, we dereddened the colors of our sample with the extinction values determined by Manara et al. (2012). The resulting color-color diagram is shown in Fig. 3 and illustrates how the colors are distributed when neglecting extinction. Furthermore, Fig. 3 shows the predicted displacement of the sources from the line describing the photosphere according to Eq. (5). For a given T_eff value, a track in the color-color diagram is determined by varying η. The track starts for η = 0 at the line describing the pure photospheric emission and converges at the colors of the pure accretion spectrum for increasing η. The plot also shows lines of constant ratio between the accretion luminosity and the total luminosity, L_acc/L_tot. The calculation of these lines is detailed in Appendix E.

Fig. 3. Color-color diagram of the WWFI sample, dereddened with the extinction values determined by Manara et al. (2012) and plotted as 2D Gaussians according to their photometric uncertainties. The values are colored according to their Teff value. Overlaid are the mean colors for each star as solid dots. The cross in the upper left corner indicates their typical range of variability. The thick line shows the photospheric colors for non-accreting stars. The thin lines represent tracks along which the sources are displaced from the photospheric track for increasing accretion. They converge in a point representing the pure accretion spectrum, described by a hydrogen slab model. The dotted lines indicate constant accretion to total luminosities ratios, Lacc/Ltot.

The model with the parameters T_slab = 11 000 K, τ₃₀₀ = 5.00, and n_e = 10¹⁵ cm⁻³ lead to tracks describing the extinction-corrected colors best. The corresponding magnitudes are m_{acc,  u′} = 10.70 mag, m_{acc,  g′} = 10.76 mag and m_{acc,  r′} = 10.60 mag, while the accretion luminosity amounts to log(L_acc,  0/L_⊙) = 1.04. We use these parameters for our subsequent analysis.

3.3. Accretion luminosities and extinctions from observations

We applied Eq. (5) to each u′g′r′ observation and derived η and A_V using T_eff, the photospheric flux, and the accretion model. The observed position in the color-color diagram is shifted along the reddening vector, until it intersects the unique track calculated from T_eff. This procedure provides A_V as well as η. The yet missing radius R_* can be calculated from the observed magnitudes, A_V, and η with Eq. (4). From the definition of η, the accretion luminosity follows from

$$\begin{array}{c}\hfill L_{\mathrm{acc}}=\eta ·{\left(\frac{R_{\ast }}{R_{\odot }}\right)}^2·L_{\mathrm{acc},0}.\end{array}$$(6)

In the next subsection, we explain in more detail how we estimated the uncertainties of A_V, R_*, and L_acc.

3.4. Uncertainties

In order to estimate the uncertainties, we used a Monte Carlo approach similar to Manara et al. (2012): the colors were randomly displaced according to their photometric uncertainty assuming Gaussian. The same was done with T_eff, assuming an uncertainty corresponding to one stellar subclass. A limitation of our approach is the possibility that the spectral types of the stars in our sample changed by more than one stellar subclass in the time span between their determination and our photometric monitoring.

The 1σ uncertainty of the photospheric magnitudes is assumed to be the 1σ confidence interval. The accretion luminosity and the extinction were then determined with the method described above. This step is repeated 10⁴ times and the 1σ uncertainty is defined as the standard deviation of the obtained set of values. With this method, we estimated uncertainties for A_V, R_*, and L_acc.

If there was more than one solution, that is the line intersects the track at more than one point, we chose the solution with the lower accretion luminosity. For some iterations, no solution could be found. This is the case when the colors are not inside the range of the combination of the model accretion spectrum and the extinction. Similar to Manara et al. (2012), we assigned the results a confidence level according to the amount of successful intersections: if there was a successful intersection for all 10⁴ iterations, we assigned the result a confidence level of 3σ and a lower confidence level of 2σ and 1σ if there were only successful intersections for more than 95% and 68% of the iterations, respectively.

For the further analysis, we only used sources with a confidence level ≥1σ and A_V ≥ 0 mag. Furthermore, we limited ourselves to stars for which at least five epochs are available.

3.5. Accretion rates and stellar parameters

The A_V values obtained with our method range from 0.00 mag to 3.04 mag, with a mean value of 0.70 mag. For each star, we calculated the mean of the extinction values and compared the results with the extinctions obtained by Manara et al. (2012). We found a good match based on the distribution of the differences ΔA_V = ⟨A_{V,  this work}⟩ − A_{V,  Manara}, which displays a median value of −0.09 mag and a standard deviation of 0.57 mag.

With the given T_eff values and the radii obtained with our method, the bolometric luminosities L_* can be calculated via L_*/L_⊙ = (R_*/R_⊙)² ⋅ (T_eff/T_⊙)⁴. From the bolometric luminosities and the known spectral types, we estimated the stellar masses M and ages τ by interpolating PARSEC 1.2S isochrones and mass tracks from Bressan et al. (2012) on the Hertzsprung-Russel diagram (HRD), shown in Fig. 4. Under the assumption that the gravitational infall releases its energy at a magnetospheric radius of ∼5 stellar radii (Shu et al. 1994), the accretion rate can be estimated from L_acc (e.g., Hartmann et al. 2016):

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{acc}}\approx \frac{L_{\mathrm{acc}}{R}_{\ast }}{0.8GM}.\end{array}$$(7)

Fig. 4. Hertzsprung-Russel diagram of the observed sources. Overlaid are the PARSEC 1.2S isochrones and mass tracks from Bressan et al. (2012). The dots represent the mean bolometric luminosity for each star as a function of Teff. The vertical bars indicate the minimum and maximum range of L* for each star. The cross at the bottom shows the typical uncertainty. In our analysis, we only regarded sources with 0.2 M⊙ ≤ M ≤ 2.0 M⊙ and 5.5 ≤ log(τ[yr]) ≤ 7.3 in order to exclude outliers in the Hertzsprung-Russel diagram.

Similar to Sect. 3.4, the uncertainties of the masses, ages, and accretion rates were estimated with a Monte Carlo approach, propagating the uncertainties from the photometry, the effective temperatures, and the photospheric magnitudes. The results are summarized in Table 3.

3.6. The final sample

For the further analysis, we calculated the ratio L_acc/L_* and compared it with relation (1) derived by Manara et al. (2017a) in order to exclude values that are likely below the “noise level” due to photospheric emission. In addition, we limited the sample according to 0.2 M_⊙ ≤ M ≤ 2.0 M_⊙ and 5.5 ≤ log(τ[yr]) ≤ 7.3 in order to exclude outliers in the HRD.

Since we are primarily interested in the influence of the accretion rate variability on the relation between X-ray emission and accretion, we limited our sample further to stars with a match with the X-ray data from the COUP catalog (Getman et al. 2005). This left us with a total of 23 stars with ∼12 mass accretion rate values per average, good u′g′r′ photometry for ≥5 nights, Gaia distances consistent with the ONC, evidence of accretion, and a confidence level ≥1σ. For each star, we also obtained the mean values of the stellar parameters and propagated their uncertainties. We report the resulting values in Table 4.

4. Results

4.1. Accretion rate variability

In order to probe the variability of the obtained accretion rates, we calculated the range of variability of the accretion rates $\mathrm{\Delta }log\left({\dot{M}}_{\mathrm{acc}}\right)$, defined as the difference between the maximum and the minimum logarithmic accretion rate for each target. In addition, we calculated the interquartile range (IQR) Δ_IQR log (Ṁ_acc), which is defined as the difference between the upper and the lower quartile.

To estimate the uncertainties, we used again a Monte Carlo approach: the mass accretion rates were varied according to their uncertainty assuming Gaussian and Δlog(Ṁ_acc) as well as Δ_IQR log(Ṁ_acc) was calculated for each source. This step was repeated 10⁴ times and the uncertainty estimated as the 1σ standard deviation drawn from the obtained distributions. The resulting values are listed in Table 4.

Figure 5 shows a box plot of the mass accretion rates for each source, subtracted by their median. The range of the variability varies between 0.13 dex and 1.89 dex, with a median value of ∼0.54 dex. We estimated the uncertainty of the median variability with the same method as described above and obtained ∼0.10 dex. The accretion rate variability is lower than the scatter found in F21 for the Ṁ_acc − M relation of 0.80 dex, defined as the standard deviation of the residua obtained by subtracting the regression line from the observed Ṁ_acc values. The values agree within 1σ with the typical mass accretion rate variability of ∼0.5 dex determined for week-timescales (e.g., Venuti et al. 2014) and is lower than the value of ∼0.65 dex in the time range of several months found by Nguyen et al. (2009) using Hα equivalent widths. Our results support previous findings that the major contribution to the accretion variability is found on shorter timescales than the year-timescale. The IQR varies between 0.04 dex and 0.65 dex, with a median value of 0.30 dex and a typical uncertainty of ∼0.04 dex.

Fig. 5. Boxplot showing the accretion rate variability of the targets. The vertical solid lines indicate the range of the accretion rates and the height of the boxes the interquartile range. The circles show the accretion rates for each target and epoch. The values are normalized to the median values and color-coded according to their effective temperature.

A similar consideration regarding the extinction yields Δ_IQRA_V = (0.22 ± 0.04) mag, a variability compatible within 1σ with the variability of the mass accretion rates. The IQR variability of the bolometric luminosity reads Δ_IQRlog(L_*/L_⊙) = 0.06 ± 0.01, in accordance with the typical uncertainty of ∼0.07 dex we obtained for the (logarithmic) bolometric luminosity with our Monte Carlo method. We also calculated the IQR variability of the obtained masses and found Δ_IQRM = (0.02 ± 0.01) M_⊙.

In order to probe the accretion rate variability for all the timescales available in our sample, we followed the approach of Costigan et al. (2012, 2014) and calculated for each object the difference between each observation timestamp and all other observation timestamps. Then, we plotted the respective mass accretion rate differences as a function of the time differences. The result is shown in Fig. 6. Furthermore, we determined the mean of the mass accretion rate differences for bins of equal logarithmic duration and propagated the uncertainties. Although the data are irregularly and undersampled, one can deduce that the majority of the mass accretion rate variability is given by the days-to-weeks timescale.

Fig. 6. Differences of the logarithmic mass accretion rates for two epochs TA and TB as a function of the logarithmic time differences TB − TA for each measurement. The values are plotted as squares, where for each distinct object a different color was chosen according to its effective temperature. Overlaid are the mean values of the logarithmic mass accretion rate differences for different time bins as magenta columns, together with their estimated 1σ uncertainties.

4.2. Relation between accretion rates and stellar mass

We performed a linear regression in order to analyze the relation between the (logarithmic) mean masses and the mean (logarithmic) mass accretion rates. To this aim, we used the fully Bayesian LINMIXERR approach based on the method developed by Kelly (2007) that takes uncertainties in both variables into account as well as intrinsic scatter. The resulting relation reads

$$\begin{array}{c}\hfill log\left(\frac{⟨{\dot{M}}_{\mathrm{acc}}⟩}{M_{\odot }{\mathrm{yr}}^{-1}}\right)=(-8.05\pm 0.15)+(1.07\pm 0.43)·log\left(\frac{M}{M_{\odot }}\right),\end{array}$$(8)

with a linear correlation coefficient of r = 0.51 ± 0.19, indicating a positive correlation and in agreement with the slope of 1.07 ± 0.22 we found in our analysis of 332 accreting sources in F21.

In Fig. 7, the regression result is displayed as a dotted line. We assume that the agreement of the regression results from our smaller sample with the results of the larger sample used in F21 indicates that it is not strongly affected by selection effects and therefore representative for the population of accreting stars in the ONC.

Fig. 7. Logarithmic mean mass accretion rates vs. logarithm of the masses shown as solid dots. The solid bars indicate the maximum and the minimum accretion rates for each source. The dotted line shows the result of the regression obtained with the Bayesian LINMIX_ERR method. The blue cross shows the typical uncertainty of the values, the red bar the typical variability range.

4.3. Relation between X-ray activity and accretion

The main goal of this work is to check whether the accretion variability introduces a bias in the detected anticorrelation between X-ray luminosities L_X and mass accretion rates Ṁ_acc, which can be interpreted as a signature of X-ray driven photoevaporation. Under the assumption that the mean accretion rates over the timescale of ∼2 yr are more characteristic than a single “snapshot” of the accretion rates obtained during a single observation, we repeated the analysis using the mean values. We utilized a partial regression analysis for this task. The method is described in more detail in F21 and can be summarized as follows: since the accretion rates are correlated with the stellar mass (as was shown in numerous works, e.g., Alcalá et al. 2017 and the previous section) and the X-ray luminosities as well (e.g., Preibisch et al. 2005; Telleschi et al. 2007), the values must be corrected for the common mass dependence in order to avoid spurious correlations. This is done by obtaining the linear regression result of the Ṁ_acc − M and the L_X − M relations and calculating the residuals. If there was no intrinsic relation between L_X and Ṁ_acc, the linear regression of the residuals would indicate no correlation.

To minimize the effect of X-ray variability, we tried to correct the influence of flaring activity from the X-ray luminosities using the “characteristic count rates” determined by Wolk et al. (2005) as the quiescent levels of X-ray emission outside the time windows with significant flaring activity in the light curves of the COUP sources. The corrected values are denoted “characteristic” X-ray luminosities, L_X,  char. Additionally, we corrected the values for the more recent Gaia distances listed in Table 1. We assumed a typical uncertainty of L_X,  char of 0.15 dex derived from the spectral fits used to obtain the X-ray luminosities (Preibisch et al. 2005). Next, we performed the linear regression with LINMIXERR and found the relation

$$\begin{array}{c}\hfill log\left(\frac{L_{\mathrm{X}},char}{\mathrm{erg}{\mathrm{s}}^{-1}}\right)=(30.33\pm 0.15)+(1.54\pm 0.45)·log\left(\frac{M}{M_{\odot }}\right),\end{array}$$(9)

with a linear correlation coefficient of r = 0.65 ± 0.16, indicating a significant positive correlation. The slope is flatter than the one obtained in our previous analysis with the larger sample (2.08 ± 0.16).

We calculated the residuals using Eqs. (8) and (9), propagated the uncertainties and performed the linear regression with LINMIXERR. We found the following relation:

$$\begin{array}{c}\hfill log\left(\frac{⟨{\dot{M}}_{\mathrm{acc}}⟩}{{\dot{M}}_{\mathrm{acc}}(M)}\right)=0.00\pm 0.11+(-0.22\pm 0.39)·log\left(\frac{L_{\mathrm{X},\mathrm{char}}}{L_{\mathrm{X},\mathrm{char}}(M)}\right).\end{array}$$(10)

Also, the negative slope suggests an anticorrelation, the large uncertainty of ±0.39 indicates that the sample size is too small in order to draw conclusions about the correlation of the residuals. The linear correlation coefficient of r = −0.19 ± 0.32 points in this direction as well. This result is kind of expected since the observed anticorrelation in the larger sample of F21 is weak, as suggested by the models of X-ray driven photoevaporation.

Since it is known that the accretion rates as well as the X-ray luminosities tend to decrease with increasing age (e.g., Preibisch & Feigelson 2005; Telleschi et al. 2007; Manara et al. 2012), it can be expected that the scatter in the relation can be reduced by correcting not only for the mass but also for the isochronal age. Therefore, we fitted planes to the quantities with a standard Levenberg-Marquardt regression and regarded both the mean mass and the mean isochronal age simultaneously. Uncertainties were not taken into account in this step. Then, we performed the partial regression with the residuals calculated from these relations and found:

$$\begin{array}{cc}& log\left(\frac{⟨{\dot{M}}_{\mathrm{acc}}⟩}{{\dot{M}}_{\mathrm{acc}}(M,\tau )}\right)\hfill \\ \hfill & =0.00\pm 0.09+(-0.37\pm 0.24)·log\left(\frac{L_{\mathrm{X},\mathrm{char}}}{L_{\mathrm{X},\mathrm{char}}(M,\tau )}\right)\hfill \end{array}$$(11)

with a linear correlation coefficient of r = −0.36 ± 0.23. A scatter plot of the residuals together with the regression result is shown in Fig. 8. Equation (11) shows a 1σ significant anticorrelation between the residual accretion rates and the residual X-ray luminosities, but the sample is too small to draw a stronger conclusion. The interpretation of this result is further difficult due to the uncertainties of the isochronal age estimates (e.g., Preibisch 2012; Soderblom et al. 2014) and their correlation with the accretion rates (Da Rio et al. 2014).

Fig. 8. Residual accretion rates vs. residual X-ray luminosities for the ONC taking both the mean mass and the mean isochronal age dependence into account. The dotted line shows the best fit obtained with LINMIX_ERR. The solid bars indicate the maximum and the minimum residual accretion rates for each source.

4.4. Influence of accretion variability on a larger sample

In order to take advantage of a larger sample, we harkened back to our previous analysis in F21 and tried to incorporate our findings about the accretion variability in this study. To do so, we set the uncertainty of the accretion rates to twice the typical range of the accretion rate variability of ∼0.54 dex we found in our monitoring survey and repeated the partial linear regression analysis. It is not likely that the sample accommodates targets with substantially larger variability. Young stars of FU Ori type are known to produce changes in the accretion rate for several orders of magnitude due to bursting events (Hartmann et al. 2016), but with typically ∼10⁻⁵ M_⊙ yr⁻¹ (e.g., Banzatti et al. 2015), the accretion rates of these targets are more than one order of magnitude higher than the highest value measured in the HST sample used in F21.

Nevertheless, we intended to include the possibility of such rare bursting events in our analysis and chose a Monte Carlo approach: conservatively, we assumed that 5% of the targets in our sample were subject to strong variability. Therefore, we assigned 16 of the 322 sources a scatter of 5 dex and repeated the partial regression analysis. We iterated this process 1000 times and randomly shuffled the 16 sources to which the 5 dex scatter has been assigned for every iteration step. Each time, we stored the linear correlation coefficients drawn from the LINMIXERR posterior distributions, consisting of 12 600 values. From their distribution, we determined the mean and the standard deviation to obtain r = −0.74 ± 0.25. We can conclude with 95% confidence that −1.00 ≤ r ≤ −0.24, meaning that the probability is ∼95% that the null hypothesis (no or a positive correlation) is “not” true. Thus, the anticorrelation is significant, even under the conservative assumption of strong accretion variability contamination of our sample.

5. Summary and conclusions

We presented the results of a multiyear photometric monitoring of accretion rates of young stars in the ONC. The observations were carried out with the WWFI instrument and comprises photometry in the u′, g′, and r′ filters. Accretion rates were estimated from the observed displacement of the source positions in a color-color diagram by modeling the colors as a combination of an empirically determined, photospheric contribution and an accretion model.

The results allowed us to study how the accretion rate varies on timescales up to a few years. We found a typical interquartile range of ∼0.3 dex. We showed that the accretion rate variability has likely not introduced a bias in our previous study of the relation between X-ray activity and accretion rates (Flaischlen et al. 2021), where we reported a weak anticorrelation between the quantities, supporting the theoretical models of X-ray driven photoevaporation (Ercolano et al. 2008b, 2009; Drake et al. 2009; Owen et al. 2010; Picogna et al. 2019).

Acknowledgments

We wish to thank the referee for helpful suggestions. We thank the LMU master physics student Benedikt Mayr for his contribution to the WWFI data reduction. This work made use of data obtained at the Wendelstein Observatory. The 2 m telescope project is funded by the Bavarian government and by the German Federal government through a common funding process. Part of the 2 m instrumentation including some of the upgrades for the infrastructure and the 40 cm telescope housing were funded by the Cluster of Excellence “Origin of the Universe” of the German Science foundation DFG. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under DFG project number 325594231 in the context of the Research Unit FOR 2634/1: “Planet Formation Witnesses and Probes: TRANSITION DISKS”. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 823823 (DUSTBUSTERS). This research was partly supported by the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094-390783311. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the SIMBAD database and the VizieR catalog services operated at Strasbourg astronomical Data Center (CDS).

References

Appendix A: Observation log

Table A.1 lists how many targets were observed per epoch and filter.

Appendix B: Derived filter properties

For AB magnitude systems like the SDSS system WWFI uses, the zero-point Z_ζ is given by

$$\begin{array}{c}\hfill Z_{\zeta }=2.5log\left(\frac{{\lambda }_{p,\zeta }^2}{c}\right)+48.6,\end{array}$$(B.1)

with the pivot wavelength λ_p,  ζ, which can be calculated via

$$\begin{array}{c}\hfill {\lambda }_{p,\zeta }^2=\frac{\int \lambda T_{\zeta ,\lambda }\mathrm{d}\lambda }{\int {\lambda }^{-1}{T}_{\zeta ,\lambda }\mathrm{d}\lambda },\end{array}$$(B.2)

where T_ζ, λ is the instrumental response function of the respective filter and the system it is installed on (Casagrande & VandenBerg 2014). We used T_ζ, λ as given by Kosyra et al. (2014).

We also calculated the ratio of the extinction through each filter and the visual extinction, A_ζ/A_V. To this aim, we assumed a constant flux, an approximation justified in Appendix C, and the reddening law of Cardelli et al. (1989) as well as a galactic reddening parameter of R_V = 3.1, which was found to be appropriate for our sample (Da Rio et al. 2010). Using Eq. 2, we calculated Δm_ζ,  AV := m_ζ,  AV − m_{ζ,  AV  =  0 mag} for A_V values between 0 and 5 mag in steps of 0.1 mag. Then, we fitted a straight line of origin through the Δm_ζ,  AV - A_V relation and deduced A_ζ/A_V from the slope. The calculated parameters of the filter system are summarized in Table B.1.

Appendix C: Accuracy of the extinction estimates

In the described approach, we assumed that extinction moves the points in the color-color diagram along lines parallel to the reddening vector. However, the reddening vectors are not exactly parallel. This would be the case if the flux density was constant over the filter’s bandwidth, which is not the case in reality. In order to get more precise solutions, the term 10^−0.4A_λ has to be multiplied to the total flux density defined by Eq. 1. We repeated the steps of section 3 for a set of A_V values between 0 and 10 mag, in steps of 0.01 mag. This renders basically Table D.1 into a cube, from which the photospheric magnitudes can be interpolated with reasonable precision for any A_V in the chosen interval. Similarly, the colors associated with the accretion model can be obtained.

For a given T_eff, the missing two values η and A_V can be obtained by comparing the calculated colors with the observed ones and regarding the relations as a system of two nonlinear equations, which can be solved with a Quasi-Newton method, for instance. We found that the A_V and L_acc values obtained with this more precise method differ not significantly from the values gained with the simpler method for A_V < 5 mag. Since it is computationally much less expensive, we chose to use the simpler method for our analysis.

Appendix D: The photospheric magnitudes

Figure D.1 shows the magnitudes m_{phot,  u′}, m_{phot,  g′}, and m_{phot,  r′} obtained from the photospheric template spectra as described in Sec. 3.2.1 as a function of T_eff. The fitting results are listed in Table D.1.

Fig. D.1. Results of the synthetic photometry with the photospheric template spectra observed by Manara et al. (2013, 2017a). Solid lines show the best fits and the shaded regions the 95 % confidence interval. The lines can be recreated by cubic spline interpolation of the values listed in Table D.1. Shown are the u′g′r′ magnitudes as a function of Teff.

Appendix E: Lines of constant accretion through total luminosity

In Fig. 3, lines of constant ratio L_acc/L_tot are shown. In order to calculate them, one can replace the factor (R_*/R_⊙)² in Eq. 6 with L_*/L_⊙ ⋅ (T_⊙/T_eff)⁴ and solve for η. Substituting L_* in the resulting expression with L_tot − L_acc yields

$$\begin{array}{c}\hfill \eta \left(T_{\mathrm{eff}}\right)=\frac{1}{(L_{\mathrm{tot}}/L_{\mathrm{acc}}-1)·L_{\mathrm{acc},0}/L_{\odot }}·{\left(\frac{T_{\mathrm{eff}}}{T_{\odot }}\right)}^4.\end{array}$$(E.1)

All Tables

All Figures

Fig. 1. WWFI u′g′r′ color composition of the ONC (image credit: Wendelstein Observatory). Overlaid are the stars monitored in all three bands. The yellow symbols indicate stars that have a match with the X-ray sources listed in the COUP catalog (J/ApJS/160/319/coup, Getman et al. 2005). The star symbols further indicate sources that were monitored for at least five epochs and allowed an estimate for the mass accretion rate.

In the text

Fig. 2. g′−r′ vs. u′−g′ color-color diagram of YSOs monitored with WWFI as gray dots. The bigger dots indicate the objects displaying accretion according to Manara et al. (2012). The gray cross shows the typical uncertainty of the photometry. The thick line is the sequence tracing the colors of the photosphere, obtained by synthetic photometry on spectra of diskless class III objects observed with X-shooter (Manara et al. 2013, 2017b). The colors of these template spectra are shown as green squares. The arrow traces the reddening vector for AV = 1 mag, assuming the reddening law of Cardelli et al. (1989) and a galactic reddening parameter of RV = 3.1.

In the text

Fig. 3. Color-color diagram of the WWFI sample, dereddened with the extinction values determined by Manara et al. (2012) and plotted as 2D Gaussians according to their photometric uncertainties. The values are colored according to their Teff value. Overlaid are the mean colors for each star as solid dots. The cross in the upper left corner indicates their typical range of variability. The thick line shows the photospheric colors for non-accreting stars. The thin lines represent tracks along which the sources are displaced from the photospheric track for increasing accretion. They converge in a point representing the pure accretion spectrum, described by a hydrogen slab model. The dotted lines indicate constant accretion to total luminosities ratios, Lacc/Ltot.

In the text

Fig. 4. Hertzsprung-Russel diagram of the observed sources. Overlaid are the PARSEC 1.2S isochrones and mass tracks from Bressan et al. (2012). The dots represent the mean bolometric luminosity for each star as a function of Teff. The vertical bars indicate the minimum and maximum range of L* for each star. The cross at the bottom shows the typical uncertainty. In our analysis, we only regarded sources with 0.2 M⊙ ≤ M ≤ 2.0 M⊙ and 5.5 ≤ log(τ[yr]) ≤ 7.3 in order to exclude outliers in the Hertzsprung-Russel diagram.

In the text

	Fig. 5. Boxplot showing the accretion rate variability of the targets. The vertical solid lines indicate the range of the accretion rates and the height of the boxes the interquartile range. The circles show the accretion rates for each target and epoch. The values are normalized to the median values and color-coded according to their effective temperature.
In the text

Fig. 6. Differences of the logarithmic mass accretion rates for two epochs TA and TB as a function of the logarithmic time differences TB − TA for each measurement. The values are plotted as squares, where for each distinct object a different color was chosen according to its effective temperature. Overlaid are the mean values of the logarithmic mass accretion rate differences for different time bins as magenta columns, together with their estimated 1σ uncertainties.

In the text

	Fig. 7. Logarithmic mean mass accretion rates vs. logarithm of the masses shown as solid dots. The solid bars indicate the maximum and the minimum accretion rates for each source. The dotted line shows the result of the regression obtained with the Bayesian LINMIX_ERR method. The blue cross shows the typical uncertainty of the values, the red bar the typical variability range.
In the text

	Fig. 8. Residual accretion rates vs. residual X-ray luminosities for the ONC taking both the mean mass and the mean isochronal age dependence into account. The dotted line shows the best fit obtained with LINMIX_ERR. The solid bars indicate the maximum and the minimum residual accretion rates for each source.
In the text

	Fig. D.1. Results of the synthetic photometry with the photospheric template spectra observed by Manara et al. (2013, 2017a). Solid lines show the best fits and the shaded regions the 95 % confidence interval. The lines can be recreated by cubic spline interpolation of the values listed in Table D.1. Shown are the u′g′r′ magnitudes as a function of Teff.
In the text