Testing the Linear Relationship between Black Hole Mass and Variability Timescale in Low-luminosity AGNs at Submillimeter Wavelengths

We designed a monitoring observation toward six nearby LLAGNs on a small patch of sky (R.A. ∼1 hr) in a weekly-to-daily observation frequency. Six LLAGNs, Centaurus A (hereafter, Cen A), NGC 4374, NGC 4552, NGC 4579, NGC 4278, and NGC 5077, were selected. These sources are predicted to have variability timescales ranging from a few days to 10 days due to their black hole masses being situated between Sgr A* and M87. This is expected to provide clarity on the submillimeter mass scaling relationship. In addition, these sources are bright in the submillimeter wavelength (≥ 25 mJy) and show similar radio properties to Sgr A*, M81, and M87, that is, a flat-to-inverted spectrum from the centimeter to submillimeter wavelength (Doi et al. 2011).

Six targets were monitored using SMA from 2015 June to 2019 January. We adopted the "piggybacked" observation on the full-track schedule to collect the light curve; therefore, observational frequency varied around 230 GHz. Note that the angular resolution varies from ∼1'' to ∼10'', which is related to the array configuration schedule. The observations were done under good weather conditions with precipitation water vapor <4 mm.

We performed initial flagging and calibration, including bandpass, flux, and gain calibration in IDL software Millimeter Interferometer Reduction (MIR). We then performed imaging and flux measurements in the Miriad package. We carried out flux calibration using a planet moon, mainly Callisto, Titan, or Ganymede. Bright sources, such as 3C 279, 3C 454.3, or 3C 84, were used for the bandpass calibration. The choice of flux and bandpass calibrator depends on the available source in the sky. We performed gain calibration with calibrator 1315-336 for Cen A, 1159 + 292 for NGC4278, 1305-105 for NGC 5077, and 3C 274 for NGC 4374, NGC 4552, and NGC 4579. After the calibration, each receiver and sideband were imaged separately. Source detection was confirmed by both visual and a three S/N threshold applied on an individual image. Once the detection was confirmed, we applied a point-source model fit on the visibility. The measurement of flux and uncertainty from the visibility was then adopted into the light curve. Fitting flux directly on the visibility gets rid of the extra uncertainty that may be introduced by the imaging process. The light-curve data points will be filtered out when a noticeable mismatch in the calibrator's flux, which cannot be explained by the calibrator's variability, is detected between the measurement and the archive. All the targets in the individual map showed a point-source-like structure, and no extended emission has been detected.

In addition to the monitoring observation, we compiled data from the SMA calibrator database ⁸ and ALMA archive for M87, Cen A, NGC 4552, and NGC 4579. In the cases of M87 and Cen A, the SMA calibrator database includes observations from 2003 April and 2005 July to the present, spanning about 18 and 15 yr, respectively. We adopted NGC 4552 and NGC 4579s ALMA observations from the Japanese Virtual Observatory (JVO) image archive. ⁹ We measured the flux density from the FITS images in the JVO image archive by modeling the source as a point source, and we estimated the uncertainty by taking the rms value from the off-source region. Furthermore, we included a measurement of NGC 4579 with a flux of 11 ± 2 mJy observed by the Plateau de Bure Interferometer (PdBI; Krips et al. 2007) around 2002. These archive data significantly expand the light-curve duration and provide useful information for understanding the long-term flux level of our targets.

Light curves of M87, Cen A, NGC 4374, NGC 4278, NGC 5077, NGC 4552, and NGC 4579 are shown in Figure 1. We use symbols with different colors to denote data from either our observation or archive. The full light-curve information is shown in Table 1. We present our unpublished SMA data in Table 2.

Zoom In Zoom Out Reset image size

Table 1. Information of Full Light Curve

Source	Start Date	End Date	Consist	Mean Flux	Var a
	(MJD)	(MJD)	(#)	(Jy)	(%)
M87	52736.0	59789.0	SMA: 83, SMACAL: 231	1.571 ± 0.261	15
Cen A	53564.0	59267.0	SMA: 44, SMACAL: 61, ALMA: 9	7.743 ± 1.766	22
NGC 4278	57174.0	58484.0	SMA: 69	0.036 ± 0.007	19
NGC 4374	57174.0	58484.0	SMA: 83	0.133 ± 0.024	18
NGC 4552	56803.0	58484.0	SMA: 22, ALMA: 2	0.023 ± 0.005	16
NGC 4579	52379.0	58553.0	SMA: 19, ALMA: 5, PdBI: 1	0.037 ± 0.031	84
NGC 5077	57174.0	58484.0	SMA: 75	0.120 ± 0.033	27

Note.

^a The mean fraction of variation is defined as Var$=\sqrt{{\sigma }^{2}-{\delta }^{2}}/\langle F\rangle$, where σ² is the variance of the light curve, ${\delta }^{2}=1/N{\sum }_{i}^{N}{\delta }_{i}^{2}$ is the mean square of flux uncertainties, and 〈F〉 is the mean of flux.

Download table as:  ASCIITypeset image

The composition of each light curve is shown in Table 1. In addition, we have included M87 (3C 274) in this study. The variability properties of M87 were investigated in Bower et al. (2015) using an 11 yr light curve with 147 data points. In this paper, we have extended the time baseline to 18 yr with 314 data points. We further included M87's flux measurement from our observation, which provides a finer cadence for probing M87's variability.

Cen A and NGC 4374 show long-term stability, which is consistent with our understanding of LLAGN variability properties (Dexter et al. 2014; Bower et al. 2015). NGC 4552 and NGC 4579 also show stability during their whole light curve; however, note that we only have a small number of data across a short period. This suggests that the result in the following analysis section should be treated with caution. On the other hand, flux levels gradually increase for NGC 4278 and NGC 5077, suggesting that these two sources may behave more like a regular AGN. For M87, we observed a long declining trend starting from 2015, but it reversed itself and reached the mean again recently, which was not included in the previous study, making it worth reanalyzing.

We have applied a first-order structure function (SF) analysis to the light curves shown in the previous section. The structure function SF(Δt_i ) of a light curve F_j with N numbers of data points in the Δt_i time bin is defined as

$$\begin{eqnarray}&&\mathrm{SF}{\left({\rm{\Delta }}{t}_{i}\right)}^{2}=\displaystyle \frac{1}{N({\rm{\Delta }}{t}_{i})}\displaystyle \sum _{j}{\rm{\Delta }}{F}_{j}^{2}({\rm{\Delta }}{t}_{i}).\end{eqnarray} \tag{ 1 }$$

The SF results are shown in Figure 2. Error bars, deriving from calculations in Simonetti et al.'s (1985) Appendix, show one σ of measurement noise variance plus the combined statistic error in each bin. However, the SF method has been shown to suffer from the dependent sampling issue, that is, each of the SF points is correlated to others, and the numbers of independent sampling are finite. For example, a 1000 day observation duration can only have five independent measurements for a process that operates on a timescale of 200 days. This sampling issue will be more significant at the long timescale. In this study, we estimated this uncertainty by considering the number of independent measurements (N_ind) in each bin. The SF in each bin has a fractional uncertainty of  $1/\sqrt{{N}_{\mathrm{ind}}}$, which is displayed as the gray shaded region in Figure 2.

Zoom In Zoom Out Reset image size

The SF provides observations to flux variation between different time lags; therefore, it can reflect the slope properties of the PSD on the short timescale as well as detect the variability saturation on the long timescale. A typical SF is composed of three parts: a flat component in the shortest timescale representing the noise level of the measurement; an increasing trend with a constant slope on a short timescale; and a slope changing to a flat plateau, which reflects the saturation of the variability on a long timescale (see, e.g., Figure 8 in Wielgus et al. 2022). In our SF results, M87, Cen A, and NGC 4374 show an apparent changing slope on the timescale of 10 to 100 days. On the other hand, a suspected steep-to-flat characteristic is detected in NGC 4278's SF, while no detection is found for NGC 5077. The lack of a flat plateau for NGC 4278 and NGC 5077 suggests that there is no stable flux level in their long-term variability, as we can observe in their light curves. With only a few data points (∼20), the SFs of NGC 4552 and NGC 4579 are ambiguous. Therefore, we believe the saturated timescale of a few to 10 days in these cases could be a spurious detection.

The SF has been widely applied in light-curve variability studies (e.g., Collier & Peterson 2001; Czerny et al. 2003; Dexter et al. 2014; Kozłowski 2016a) to avoid ground-based observational bias. It provides a useful insight into the qualitative comparison between sources when the actual PSD is relatively difficult to derive. However, Eracleous et al. (2010) pointed out that the limitation of SF is related to three points. First, properly handling the estimation of the timescale and its uncertainty could be challenging because SF data points are nonindependent. Second, any specific total sampling duration may lead to a spurious saturated timescale. Most important, the SF can still be affected by the observational gap in the light curve.

We further compare the Gaussian process model fitting result (the DRW in the next section) with the SF measurement. We present the SF of a simulated DRW light curve with parameters inferred from the DRW fitting process. In Figure 2, the DRW light curve is labeled in cyan. The dashed line represents the simulated light curve with a regular sampling cadence of 1 day, while the dotted line shows the actual sampling of the SMA observation, and both simulated light curves present the actual duration of the real data. Additionally, we exhibit the DRW-inferred parameters, such as the timescale (vertical red dashed line) and the doubled long-term variance (vertical orange dashed line). The analysis suggests that for M87 and Cen A, the SF (black line) is consistent with the DRW model (cyan dashed line), and the DRW-inferred timescale is situated at the first local maximum, which is generally considered as being where the SF slope breaks. In contrast, the comparison becomes marginally consistent for sources such as NGC 4374 and NGC 4278 and inconsistent for NGC 5077 and NGC 4579. We hypothesize that the SF result is strongly biased by the sampling window effect by simulating the actual sampling DRW light curve (cyan dotted line). Therefore, we conclude that the SF is not sufficient to constrain the timescale for most of the sources.

To quantify the variable behavior of these sources, we have modeled the light curves as the result of a stochastic DRW process. DRW is a specific Markovian stationary Gaussian process. The application of DRW as a mathematical model to describe the optical variability of quasars was first proposed by Kelly et al. (2009). The validation of DRW to describe the submillimeter light curve of Sgr A* was then investigated by Dexter et al. (2014) and applied to other submillimeter sources by Bower et al. (2015). DRW contains only three parameters: the mean of the light curve μ, the long-term standard deviation σ, and the characteristic timescale τ. The structure function of the DRW model is

$$\begin{eqnarray}&&{\mathrm{SF}}^{2}({\rm{\Delta }}t)=2{\sigma }^{2}(1-{e}^{-| {\rm{\Delta }}t/\tau | }).\end{eqnarray} \tag{ 2 }$$

The τ determines where the light curve is decorrelated, which corresponds to the saturation of the variability. In the frequency domain, the PSD of the DRW model can be given as

$$\begin{eqnarray}&&\mathrm{PSD}(f)=\displaystyle \frac{4{\sigma }^{2}\tau }{1+{\left(2\pi f\tau \right)}^{2}},\end{eqnarray} \tag{ 3 }$$

which is composed of a red noise spectrum (α_PSD =−2) in the higher-frequency regime (f_τ ≫ 1) and a flat white noise spectrum at low frequencies (f_τ ≪ 1).

We have followed Kelly et al. (2009) and Dexter et al. (2014) to fit light curves with the following procedure. The likelihood function of a DRW fitting the observed light curve is given by a product of Gaussian describing the residuals from each light-curve point as

$$\begin{eqnarray}&&P(\{{x}_{i}\}| \mu ,\tau ,\sigma )=\displaystyle \prod _{i=1}^{n}{\left[2\pi ({{\rm{\Omega }}}_{i}+{\sigma }_{i}^{2})\right]}^{-\tfrac{1}{2}}\exp \left[-\displaystyle \frac{1}{2}\displaystyle \frac{{\left({\hat{x}}_{i}-{x}_{i}^{* }\right)}^{2}}{{{\rm{\Omega }}}_{i}+{\sigma }_{i}^{2}}\right],\end{eqnarray} \tag{ 4 }$$

where σ_i represents measurement uncertainties and ${x}_{i}^{* }$ is the difference between each data point and mean,

$$\begin{eqnarray}&&{x}_{i}^{* }={x}_{i}-\mu .\end{eqnarray} \tag{ 5 }$$

The quantities Ω_i and ${\hat{x}}_{i}$ are iterated by the following equations,

$$\begin{eqnarray}&&{\hat{x}}_{i}={a}_{i}{\hat{x}}_{i-1}+\displaystyle \frac{{a}_{i}{{\rm{\Omega }}}_{i-1}}{{{\rm{\Omega }}}_{i-1}+{\sigma }_{i-1}^{2}}({x}_{i-1}^{* }-{\hat{x}}_{i-1})\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{i}={{\rm{\Omega }}}_{1}(1-{a}_{i}^{2})+{a}_{i}^{2}{{\rm{\Omega }}}_{i-1}\left(1-\displaystyle \frac{{{\rm{\Omega }}}_{i-1}}{{{\rm{\Omega }}}_{i-1}+{\sigma }_{i-1}^{2}}\right),\end{eqnarray} \tag{ 7 }$$

and regulated by a_i , which is the expected correlation between data points,

$$\begin{eqnarray}&&{a}_{i}=\exp [-({t}_{i}-{t}_{i-1})/\tau ].\end{eqnarray} \tag{ 8 }$$

The iterative procedure starts from the initial conditions

$$\begin{eqnarray}&&{\hat{x}}_{1}=0,\,{{\rm{\Omega }}}_{1}={\sigma }^{2}.\end{eqnarray} \tag{ 9 }$$

Note that σ² is the variance of the light curve.

This procedure allows us to calculate the best-fit DRW solution from a given vector of parameters. Therefore, the fit quality can be evaluated by examining the distribution of the residual χ and the reduced chi-square ${\chi }_{n}^{2}$, calculated as

$$\begin{eqnarray}&&{\chi }_{\nu }^{2}=\displaystyle \frac{1}{\nu }\displaystyle \sum _{i=1}^{n}\displaystyle \frac{{\left({\hat{x}}_{i}-{x}_{i}^{* }\right)}^{2}}{{{\rm{\Omega }}}_{i}+{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 10 }$$

where ν is the degree of freedom. A proper fit should meet the criteria that the residuals are independent, characterized by a Gaussian distribution with a mean of 0 and a variance of 1, and consequently, the ${\chi }_{\nu }^{2}$ should approach 1. Figure 3 illustrates the best-fit DRW solution, the light curve, and the residual χ for seven sources. The values of ${\chi }_{\nu }^{2}$ are listed in the final column of Table 3. One can observe that the best-fit DRW solution will approach the mean of the light curve if the time interval between two data points is significantly greater than the τ, which agrees with the fitting procedure in Equations (6) and (7). We draw attention to this characteristic of the fitting outcomes, particularly for NGC 4374, since the DRW fits a smaller timescale than the sampling cadence. This results in a DRW light curve with less variability than the actual light curve observed.

Zoom In Zoom Out Reset image size

We then applied a Bayesian approach to compute the posterior probability of τ, μ, and σ from the likelihood. Following Dexter et al. (2014), we assume uniform priors on all three parameters. The Metropolis–Hastings algorithm was adopted to sample τ, μ, and σ. The corner plot of the parameters solution is presented in Figure 4 using Cen A as an example. The inferred parameters were estimated by taking the median of the sample population (the red line in Figure 4). We adopted a 2σ uncertainty, which corresponds to the 95% confidence interval (CI) in the population shown in the vertical black dashed line.

Zoom In Zoom Out Reset image size

We have applied the DRW model to measure the τ of our six new samples and to reanalyze M87. Figure 5 showcases the DRW fit solution of the probability distribution of τ, which is the most relevant parameter in the light curve. A thorough examination of the probability distribution was carried out by plotting the 68%, 95%, and 99.7% CIs for the sample. Note that a good measurement is only valid when both sides of the distribution converge at τ range from τ > Δt and τ < T, where Δt is the minimum observation separation in the light curve, and T is the total light-curve duration. We observe a long tail extending to a short or long timescale, suggesting that τ of the source is not well captured. However, in these cases, the timescale has a high likelihood of being located where the Metropolis–Hastings algorithm conducted intense sampling. As a result, we present the median and 95% CI of the distribution to reflect the highest possible parameter value in these sources.

Zoom In Zoom Out Reset image size

The inferred timescale τ, the inferred light-curve mean μ, and the inferred standard distribution σ are presented in Table 3. The reduced χ² of the best-fit results are calculated in the last column. We note that the correction for the cosmological time dilation is not necessary in these cases. The correction according to Equation (17) in Kelly et al. (2009) is ∼1% for the most distant object NGC 5077 (z ∼ 0.00936; from NASA/IPAC Extragalactic Database), which is smaller than the inferred uncertainty of the timescale.

A long-enough monitoring observation is critical for a DRW study, which provides sufficient probe to the long-term variance, i.e., the white noise part of the PDS. Combining archive data with our monitoring observation, the total light-curve duration for each source is at least 10 times longer than the τ we found, suggesting that the DRW result should be secure in our cases (Kozłowski 2017). On the other hand, the stationary character of the underlying physical process becomes a concern for such a long monitoring observation. In other words, the effect of the underlying process change should be considered.

We found short timescales to be constrained (or with a robust trend to be constrained) in all of our cases except NGC 4552. With a proper sampling cadence (both a cadence much shorter than the timescale and a duration much longer than the timescale), we observed that the degree of how well a timescale can be constrained relies on a certain number of data points. For example, in Dexter et al. (2014), Bower et al. (2015), and this study, sources, such as Sgr A*, M87, M81, and Cen A, whose light curve has over 100 data points, have been accumulated, and this seems to provide a reliable constraint on the timescale. NGC 4374, NGC 4278, and NGC 5077, with only around 80 data points accumulated, exhibit significant uncertainty in either the upper or lower limit. Furthermore, timescale τ becomes poorly constrained in the case of NGC 4579 and NGC4552. Therefore, we suggest that 100 data points seem to be a lower-limit number for measuring a reliable timescale in this experiment. Below we discuss the analysis details for individual cases.

3.3.1. Cen A

3.3.2. NGC 4374

3.3.3. NGC 4278 and NGC 5077

3.3.4. NGC 4579 and NGC 4552

3.3.5. M87

3.3.6. Limitation and Future Work

We have compiled black hole masses from the literature in Table 4. Most of the masses were measured using the multiple dynamical method, which provided a precise mass estimation. By accurately measuring stellar proper motion, Sgr A* has two independent, high-precision measurements with a mass of 4.297 ± 0.012 × 10⁶ M_⊙ using the GRAVITY instrument (GRAVITY Collaboration et al. 2022), and a mass of 3.951 ± 0.047 × 10⁶ M_⊙ using Keck (Do et al. 2019). We adopted a mass of 4.297 ± 0.012 × 10⁶ M_⊙ for Sgr A* and note that, for our purposes, the discrepancy in the mass measurement is irrelevant. We adopted a mass of 6.5 ± 0.7 × 10⁹ M_⊙ for M87, which was determined by fitting the physical parameters of lensed emissions near the event horizon (Event Horizon Telescope Collaboration et al. 2019b). This value is consistent with the stellar dynamical method, which suggests a mass of 6.6 ± 0.4 × 10⁹ M_⊙ (Gebhardt et al. 2011). M81 has two mass measurements using both the stellar and dynamic methods; therefore, we followed Kormendy & Ho (2013) to adopt a mass of ${6.5}_{-1.5}^{+2.5}\times {10}^{7}{M}_{\odot }$, which took the mean of two measurements. Cen A has a mass of 5.5 ± 3 × 10⁷ M_⊙ using the stellar dynamic method (Cappellari et al. 2009). For NGC 4374, two groups analyzed the same data but identified different black hole masses. We followed Walsh et al. (2010) who clarified this discrepancy and adopted a mass of ${8.5}_{-0.8}^{+0.9}\times {10}^{8}{M}_{\odot }$. NGC 5077 has a mass of ${6.8}_{-2.8}^{+4.3}\times {10}^{8}{M}_{\odot }$ using the gas dynamic method (de Francesco et al. 2008) and NGC 4552 has a mass of (4.7 ± 0.5) × 10⁸ determined by the stellar dynamic method (Sahu et al. 2019). In addition to dynamical methods, we followed the work done by Wang & Zhang (2003) to adopt NGC 4278 and NGC 4579's black hole masses that are estimated by the M_BH–σ relationship with a 20% uncertainty on the stellar velocity dispersion. However, we note that systematic errors can be large for all black hole mass estimates.

Table 4. Black Hole Masses Compiled from the Literature

Source	MBH	Method	References	BH Shadow Size a
	(M⊙)			(μas)
Sgr A*	(4.297 ± 0.012) × 106	Stellar proper motion	GRAVITY Collaboration et al. (2022)	∼51.3
M87	(6.5 ± 0.7) × 109	Lensed emission	Event Horizon Telescope Collaboration et al. (2019b)	∼38.2
M81	$({6.5}_{-1.5}^{+2.5})\times {10}^{7}$	Stellar and gas	Kormendy & Ho (2013)	∼1.8
Cen A	(5.5 ± 3) × 107	Stellar	Neumayer (2010)	∼1.4
NGC 4278	(3.1 ± 0.5) × 108	M–σ relation	Wang & Zhang (2003)	∼1.8
NGC 4374	$({8.5}_{-0.8}^{+0.9})\times {10}^{8}$	Gas	Walsh et al. (2010)	∼5.0
NGC 4552	(4.7 ± 0.5) × 108	Stellar	Sahu et al. (2019)	∼2.8
NGC 4579	(1.1 ± 0.2) × 108	M–σ relation	Wang & Zhang (2003)	∼0.6
NGC 5077	$({6.8}_{-2.8}^{+4.3})\times {10}^{8}$	Gas	de Francesco et al. (2008)	∼1.5

Note.

^a We assume a fiducial radius of 5R_s for the potential observable black hole shadow size. The angular diameter of the black hole shadow is calculated using 5R_s = 10GM/Dc², where D is the distance to the observer. The distances of the black holes are adopted from Sgr A*: GRAVITY Collaboration et al. (2022); M87: Event Horizon Telescope Collaboration et al. (2019b); M81: Freedman et al. (1994); Cen A: Harris et al. (2010); and others from NASA/IPAC Extragalactic Database (NED).

Download table as:  ASCIITypeset image

To qualitatively understand the physical meaning of the timescale and the connection with the black hole mass, we presented the black hole mass–characteristic timescale correlation among nine sources in Figure 6. New LLAGN samples in this study are labeled in purple, while sources measured in previous studies are labeled in different colors. Here we compare these measured timescales with the physical timescales in the accretion system. The orbital timescale t_orb and the viscous timescale t_visc, which scale with the black hole mass M_BH and Schwarzschild radius R_s , can be expressed in the following equations:

$$\begin{eqnarray}&&{t}_{\mathrm{orb}}\approx 0.5266\left(\displaystyle \frac{M}{{10}^{8}{M}_{\odot }}\right){\left(\displaystyle \frac{R}{3{R}_{s}}\right)}^{\tfrac{3}{2}}\,\mathrm{day},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{visc}}\approx 17.55\left(\displaystyle \frac{M}{{10}^{8}{M}_{\odot }}\right){\left(\displaystyle \frac{R}{3{R}_{s}}\right)}^{\tfrac{3}{2}}{\left(\displaystyle \frac{\alpha }{0.03}\right)}^{-1}{\left(\displaystyle \frac{H/R}{0.6}\right)}^{-2}\,\mathrm{day}.\end{eqnarray} \tag{ 12 }$$

The use of parameter α ∼ 0.03 − 0.05 and the scale height H/R ∼ 0.6 is based on previous simulation results (Narayan et al. 2012).

Zoom In Zoom Out Reset image size

Two physical timescales are labeled in Figure 6 with the innermost stable circular orbit (ISCO) radius (3R_s = 6GM/c²) for a Schwarzschild black hole. We found that no source falls under the orbital ISCO timescale, and most sources have a characteristic timescale comparable with the orbital or viscous timescale in the inner region of the accretion disk. Note that the new M87 result agrees with the orbital ISCO timescale since the period of ISCO varies from 5 days to about a month, relevant to the spin of the black hole from 1 to 0. NGC 4374 is situated near the orbital ISCO line. However, as previously discussed in the analysis section, the algorithm has a tendency to choose the unconstrained long tail. As a result, we anticipate that the well-constrained timescale would be larger.

Bower et al. (2015) found a linear M_BH–τ correlation among three LLAGNs, Sgr A*, M81, and M87. Our new results agree with the trend they found. Besides the lower-limit measurement for NGC 4552, the whole samples show a significant trend that the timescale of these sources increases along with their black hole mass. This can be understood by approximating the timescale scaling around the black hole to be $\tau \propto {M}_{\mathrm{BH}}{R}_{s}^{3/2}$ (a simplification of Equations (11) and (12)). One can expect a linear black hole mass–minimum timescale correlation to establish when the emission originates from the same R_s , where R_s is in the event horizon scale. In contrast, Cen A shows a much larger timescale in this population and deviates from the linear relationship. The VLBI experiment result (discussed in the next paragraph) suggests that our timescale of a few to 10 days of Cen A comes from a more outer region of the accretion disk. From our result, the black hole mass–minimum timescale correlation in the submillimeter light curve indicates that the submillimeter emission originates from the event horizon scale for these sources.

This interpretation is consistent with our understanding of the VLBI experiment results and the multiwavelength studies. The synchrotron radiation dominates the emission that originates from the accretion flow at the black hole vicinity, and the synchrotron self-absorption leads to an optically thick photosphere that shrinks along with the frequency until the spectral peak. For instance, we saw the intrinsic size of the emission core shrink along with the higher-frequency VLBI observation for M87 and Sgr A* (Hada et al. 2011; Johnson et al. 2018). The black hole image was then resolved by the EHT in the 2017 campaign using the 230 GHz global VLBI. The sizes of the rings are ∼5.5R_s and ∼4R_s , respectively (Event Horizon Telescope Collaboration et al. 2019a; Event Horizon Telescope Collaboration et al. 2022a). By contrast, an edge-brightened jet structure was observed in the 2017 EHT campaign for Cen A, and the interpreted core size is ∼32μas (∼120R_s ; Janssen et al. 2021). Based on the jet model, the optically thin frequency is expected to be located between the wavelengths from terahertz to far-infrared and further, as agreed by the Cen A core SED study (Abdo et al. 2010). In the cases of M81 and NGC 4374, Jiang et al. (2018) and Jiang et al. (2021) report the unresolved core observation with VLBA at 88 GHz. Combined with their flat-to-inverse slope spectrum, the VLBI measurement suggests the upper limit of the core size at 230 GHz to be <80R_s and <30R_s for M81 and NGC 4374. On the other hand, the submillimeter emission core size of NGC 4278, NGC 4579, and NGC 5077, lacking a high-angular-resolution observation, is unknown. However, one can interpret that the submillimeter emission is optically thin for NGC 4278 and NGC 4579 from their SED study (see Mason et al. 2013; Bandyopadhyay et al. 2019). Although the SED study is incomplete, it generally provides properties of the nonthermal emission and further supports our timescale result.

To investigate the M_BH–τ scaling relationship suggested in Figure 6, we perform both linear and power fits on our sample set. This sample set includes τ values from Sgr A* and M81 obtained from the literature and excludes the lower-limit constraint of NGC 4552. We exclude Cen A from our analysis of the mass scaling relationship because its submillimeter properties, as suggested by VLBI experiments and SED studies, differ from those of our other samples. When multiple timescale measurements are available, we use the best measurement for our discussion of the mass scaling relationship. For Sgr A*, we adopt the 8 hr measurement reported in Wielgus et al. (2022), as it provides high-quality data that were previously unattainable. Additionally, it is worth mentioning that Wielgus et al. (2022) have extensively discussed the comparison of timescale measurements obtained from different sets of light curves, which indicates the evolution of the physical state of Sgr A*. For M87, we use our 21 day measurement, as it provides the densest sampling cadence and longest time baseline in our observation.

Following Kelly (2007), we adopt a Bayesian method and apply a Metropolis–Hastings sampler to account for both M_BH and τ uncertainties in the linear regression process in both linear space and logarithmic space. The best-fit parameters, i.e., the slope and the intercept, are estimated from the median of the probability distribution, and we present 1σ as the fit uncertainty.

One can expect a τ–M_BH linear relationship ($\tau \propto {M}_{\mathrm{BH}}{R}_{s}^{3/2}$) from the physical interpretation. However, we found a linear fit result of $\tau ={2.912}_{-1.030}^{+0.783}\,(\mathrm{days})\,({M}_{\mathrm{BH}}/{10}^{8}{M}_{\odot })$ with a coefficient of determination of R² ∼ 0.23, suggesting that the relationship between τ and M_BH is not simply linear. On the other hand, the power-law fit provides a result of

$$\begin{eqnarray}&&\tau \approx {3.67}_{-1.56}^{+1.07}\,\mathrm{days}\,{\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{10}^{8}{M}_{\odot }}\right)}^{{0.55}_{-0.09}^{+0.09}},\end{eqnarray} \tag{ 13 }$$

with a coefficient of determination of R² ∼ 0.61, and we report a 1σ intrinsic scatter of ${0.74}_{-0.3}^{+0.54}$ dex for the data.

Based on our limited sample, our experiment suggests that there is a scaling relationship between the black hole mass and timescale; however, the correlation is weak and there is significant intrinsic scatter, implying that there may be a third parameter involved that may not be easy to measure, such as the optical depth. We highlight that this relationship is largely driven by Sgr A* and M87. By considering different measurements that reflect the different status of the sources' underlying physics for Sgr A* and M87, we found that the relationship can appear more steep or more shallow and tighter or looser. This evidence suggests that the robustness of the relationship requires further samples to be added to the M_BH–τ correlation, especially the timescale measurement for black hole mass equal to or smaller than Sgr A*, and black hole mass larger than M87.

In Figure 7, we present the linear regression result of our submillimeter measurement as well as the optical and X-ray linear fit from Burke et al. (2021) and González-Martín & Vaughan (2012). The optical data, containing 67 AGNs, have the best-fit model of $\tau ={107}_{-12}^{+11}\,\mathrm{days}\,{({M}_{\mathrm{BH}}/{10}^{8}{M}_{\odot })}^{{0.38}_{-0.04}^{+0.05}}$, which considers errors in both parameters. X-ray data, including 42 light-curve measurements, suggest a linear correlation fit of $\tau ={0.02}_{-0.01}^{+0.02}\,\mathrm{days}\,{({M}_{\mathrm{BH}}/{10}^{6}{M}_{\odot })}^{1.09\pm 0.21}$, but without considering the black hole mass error.

Zoom In Zoom Out Reset image size

While various types of AGNs may be classified into different accretion models, such as the standard thin disk model and the radiatively inefficient accretion flow, resulting in different wavelengths being attributed to various emission mechanisms, a comparison of variability properties between wavelengths still offers a unique opportunity to gain insight into the disk configuration and the origin of variability. The submillimeter and X-ray wavelengths show a significantly shorter timescale compared to the optical wavelength, suggesting that submillimeter and X-ray variability can provide information about the inner regions of the black hole system. On the other hand, the optical timescale is comparable to that of the outer disk's dynamics or thermal timescale. These observations of the timescales at different wavelengths are in agreement with our current understanding of the structure of the accretion disk/flow. Moreover, slopes are different in the X-ray, submillimeter, and optical M_BH–τ correlations, which come with 1, 0.6, and 0.4, respectively. In the inner part of the accretion flow, X-ray correlation indicates that plasma around the black hole (corona) correlates tightly with black hole mass. In contrast to X-ray, our submillimeter correlation may be impacted by additional physical parameters, such as the optical depth or black hole spin. The optical correlation in the outer region suggests that other factors have a significant impact on the variability or that the variability origin differs from that seen in the X-ray or submillimeter data. However, it is important to note that the X-ray sample shows more scatter than the optical and submillimeter samples, which was attributed to the absorption along the line of sight proposed by González-Martín (2018), indicating a more complex situation. Meanwhile, the submillimeter results are less scattered but remain uncertain due to being based on a smaller population.