Surveying Flux Density in Galaxies with Apparent Large Black Holes at Millimeter/Submillimeter Wavelengths

We present millimeter and submillimeter continuum observations for 36 sources with potentially large black hole shadows at 230 and 345 GHz using the Submillimeter Array. The sources are selected based on the criterion of the large diameter of the black hole shadows. Our motivation is to explore the nature of the accretion flow of the potential candidates of low-luminosity active galactic nuclei (LLAGNs) through photometry at millimeter/submillimeter wavelengths. The detected result serves as the pathfinder of future high-angular resolution observations such as the Event Horizon Telescope. As a result, we successfully detected 17 and eight sources at 230 and 345 GHz, respectively. We reveal that three of the detected sources (IC 310, NGC 1277, and NGC 5846) show significant excess at millimeter/submillimeter wavelengths in comparison with the extrapolation both from low-frequency radio and infrared, which are considered to be attributed to the extended jet and dust component, respectively. One possible explanation is that these excesses are associated with the hot accretion flow onto the supermassive black holes of the LLAGNs. By adopting the advection-dominated accretion flow model's semi-analytic model, we obtained the upper bound of the mass accretion rate. Those are less than 10−2 ṀEDD , where ṀEDD is the Eddington mass accretion rate computed via the Eddington luminosity. This is in good agreement with the expected range of the LLAGNs. The ∼200 mJy flux densities and negative spectral index of IC 1459 within millimeter and submillimeter wavelengths make it a promising candidate for future submillimeter high-angular resolution experiments for imaging the black hole shadow.


Introduction
Supermassive black holes (SMBHs) and their associated inflow/outflow system at the center of the active galactic nuclei (AGNs) are one of the most luminous systems in the universe. It is widely considered that the origin of the corresponding nuclear luminosity and presumably the outflow is attributed to the release of the gravitational energy of accreting material onto the SMBHs. Low-luminosity active galactic nuclei (LLAGNs) are a subclass of AGNs. They account for ∼40% of the nearby galaxies (Luis 2008). It is also known that they tend to host a relatively larger BH (Falcke et al. 2004), and are often categorized as radioloud AGNs (Maoz 2007). They typically do not show the feature of the big blue bump (Luis 2008) and tend to have relatively low luminosity in the optical band. And therefore, it is considered not to host an optically thick and/or geometrically thin accretion disk, which is commonly known as the standard thin disk (Page & Thorne 1974;Shakura & Sunyaev 1976). Instead, accretion materials around the central region of an SMBH in LLAGNs are considered to be relatively hot (compared to that in the case of the standard thin disk), optically thin, and geometrically thick radiatively inefficient accretion flow (RIAF; see, e.g., Yuan & Narayan 2014). The advection-dominated accretion flow (ADAF) is a submodel of RIAF, and it describes the accretion process near the central engine without outflows. In our analysis, we focus on the simplified, semi-analytical ADAF model (Mahadevan 1996) to infer the information from the AGN system such as the mass accretion rate. Note that there are two subcategories that contain the outflow process such as AGN wind (see e.g., Yuan & Yoon 2018): convection-dominated accretion flow (CDAF; Narayan et al. 2000;Quataert & Gruzinov 2000;Igumenshchev & Narayan 2002) and advection-dominated inflow-outflow solution (ADIOS; Blandford & Begelman 1999;Begelman 2012).
The spectral energy distribution (SED) of LLAGNs reveals the rich information of the radiative process near the central engine (e.g., Doi et al. 2011). One of the key predictions of the RIAF/ADAF spectrum that differs from that of standard thin disks is the existence of the millimeter/submillimeter bump (submillimeter bump). The physical nature of the submillimeter bump originates from the thermal synchrotron radiation that comes from hot accreting plasma onto the SMBHs and the inheritance properties of the LLAGNs (Mahadevan 1996). Based on the SED analysis, it is well known that Sgr A * (e.g., Bower et al. 2019, and references therein) exhibits the bump peaked around the millimeter/submillimeter wavelength, which is considered to be attributed to the thermal synchrotron radiation from the innermost region of the hot accretion flow. In addition, the sub-Eddington properties (e.g., Marrone et al. 2006) suggest its accretion flow can be categorized as the RIAF.
Although observations at millimeter/submillimeter wavelengths could be crucial to understanding the nature of hot accretion flows, such observations are limited. Therefore, we conducted Submillimeter Array (SMA) observations toward 25 potential candidates for LLAGNs at 230 and 345 GHz. Our sample consists of sources selected from McConnell & Ma (2013) and Sani et al. (2011). These sources are nearby galaxies, and their black hole masses have been previously determined based on stellar dynamics (see Table 1 for details). This sample comprises galaxies with black hole masses determined from stellar dynamics measurements, ranging from 1.4 × 10 8 to 1.7 × 10 10 M e . Their black hole shadow sizes span from 0.5-26 mas, making them promising targets for ongoing millimeter and submillimeter VLBI instruments. The sample covers a wide range of physical scales, with sources located at distances from 0.7-244 Mpc. Additionally, the sample includes galaxies of various morphological types, including LINER, Seyfert, BL Lac, FR I, and FR II, offering a diverse range of sources to study. This paper is organized as follows. We present our samples and summarize our SMA (see Ho et al. 2004) observations together with the data reduction procedures in Section 2. In Section 3, we present the results and spectral index estimates for subsamples. A discussion and the physical origin of the emissions are provided in Section 4 and summarized in Section 5.

Observations and Data Reduction
The observations were conducted in two epochs using the SMA as filler observations. The first epochs were conducted in a single band centered around 225 GHz in early 2016 on February 17, February 19, March 3, March 5, March 13, and March 15. The second epoch was held in late 2017 on September 13 and September 15 with simultaneous observations at 225 and 345 and 8 GHz bandwidth per receiver. Since the targets are widely spread across the hour angle during the Note.
(2) Logarithm of the black hole mass in units of solar mass.
(3) Distance to the source in units of megaparsec (4) Angular diameter of the black hole shadow, i.e., the photon ring of the source, which is derived from Columns 2 and 3. (5) The reference for the black hole mass. (6) AGN classifications. observing windows, the gain calibrators for each target were selected based on the following. This is mainly based on the minimum angular distance to the target from the gain calibrator, which has the maximum brightness in the SMA calibrator catalog. The flux calibrators are planets or planet satellites and are observed usually in the first scan during the observing session. The correlations were performed with the recently-coined application-specific integrated circuit correlator for the first session and the SMA Wideband Astronomical ROACH2 Machine (SWARM) correlator (Primiani et al. 2016) for the second session. The data were calibrated using the MIR-IDL data reduction packages. We follow the standard procedure for restoring system temperature, data flagging, absolute flux, bandpass, and gain calibrations. For each target, the bandpass calibrations are conducted by using the associated gain calibrators. The lower and upper sidebands are calibrated independently, and we combine them before conducting the imaging process.
Once a prior calibration was adopted, we made images for individual sources with hybrid mapping using MIRIAD. We adopt the standardized CLEAN algorithm to reconstruct the images. The imaged area is 256 by 256 pixels and the pixel size is 0 2. The noise of the image is estimated from the standard deviation in the outer region (12″ away from the image center toward the edge) of each image space, and we confirm its distribution is well described by the thermal distribution. The  observed flux density is derived from brightness as the peak pixel in the clean image, as our sources are point sources in the image domain. We quote a nominal ∼10% flux density uncertainty in our measurements (see, e.g., Chen et al. 2014;Liu et al. 2016;White et al. 2018) and incorporate it with the thermal error into the ensemble error. The results of the first epoch are summarized in Table 2, and for the second epoch, the results are in Table 3. The average synthesized beam size is ∼4 and 2 5 for 230 and 345 GHz, respectively.

Photometry
We observed a total of 25 sources at 230 GHz and nine sources at 345 GHz, detecting 17 and eight sources at 230 GHz and 345 GHz, respectively, in two different epochs. The majority of the detected sources have relatively weak flux densities, measuring less than a few millijansky. For the source we detected at 345 GHz, we also detected sources at 230 GHz, which suggests robust detection both at 230 and 345 GHz (see also Table 3). It is noteworthy that all of those sources have extended jet emission at low-frequency VLA images.
We detected several sources multiple times within the epoch through our observations. This suggests confidence in the detection. However, the light curves show significant differences among observing epochs even after accounting for the ensemble uncertainties. This suggests that the origin of the SMA flux may be related to time variations associated with black hole activity.
Fractional variability amplitude is a measure of the intrinsic variability amplitude in excess of the measurement uncertainties. It is commonly used in the study of AGN and is defined as follows: where Var[f] is the variance of the light curve, 〈σ 2 〉 is the mean squared error for individual measurements, and 〈f〉 is the mean flux density. The error on F var for N flux measurements can be defined as We perform the calculation for the source with detection multiple times and present the result in Table 4. Our analysis shows an average F var of 0.18 ± 0.11 across the 15 sources with multiple detections at 230 GHz, which is consistent with previous studies of LLAGN variability at millimeter wavelengths, such as those of Baldi et al. (2015) and Behar et al. (2020). However, it is important to note that our sources are different from those in the literature, and the frequency ranges used in different studies may also affect the variability amplitudes. Thus, direct comparisons between our results and previous studies must be made with caution. Our findings highlight the importance of considering both the range of frequency and the sources when studying millimeter-band variability in LLAGN. It is worth noting that NGC 4649 shows significant variability F var = 0.39 ± 0.09 within 1 week. Future work will aim to extend our analysis to a larger sample of LLAGN and to investigate the physical processes responsible for their variability at millimeter wavelengths.
The merged source images are presented in Figures 1 and 2 for the two epochs. The physical sizes corresponding to the unresolved compact emission for all sources detected in our observations range from 6-2300 pc. When expressed in terms of the Schwarzschild radius unit (2GM BH /c 2 , M BH is the black hole mass), the range is from 3.8 × 10 5 to 1.8 × 10 7 . These sizes were determined by modeling the source structure as a point source, which suggests that the contribution from extended components is minimal under the achieved dynamic range and image resolution of our observations. The spectral energy distributions (SED) for the detected sources at 230 and/ or 345 GHz, along with their counterparts at centimeter and far-IR (FIR) wavelengths, are shown in Figure 3.

Spectral Index Fitting
We performed spectral index fitting on sources detected at both 230 and 345 GHz using Monte Carlo-like fitting Note.
(2) Observation date. procedures. The observed flux densities were resampled, assuming a Gaussian distribution, to account for ensemble errors, which were calculated by incorporating thermal noise and nominal flux calibration errors via uncertainty propagation. Details on these errors can be found in Table 3. Specifically, we computed the spectral index at millimeter wavelengths using a power-law fit, f ∼ ν α , where f represents flux densities, ν is the LO frequency, and α is the spectral index. For each source, we performed 10,000 realizations, generating statistics such as the mean and standard deviation of the spectral index. The fitting results are shown in Figure 4, and the computed spectral index values are listed in Column 2 of Table 4. The spectral index ranges from −2 to 4 in our samples, indicating that further investigation, including low-frequency radio and FIR observations, is needed to determine the origin of the SMA flux.
To determine if there was a millimeter excess in the spectral characteristics of our samples, we used two types of spectral fits: the core spectral fit and the extended spectral fit. These methods were previously used in a study by Behar et al. (2018). The core spectral fit utilized high-angular resolution measurements, such as VLBI, to obtain spectral index fits, while the extended spectral fit used measurements with a larger beam, such as VLA. The results of these fits are presented in Column 3 of Table 4. To calculate the excess of the SMA flux, we obtained the fitted spectral index and compared it to the extrapolated flux from both types of fits. The excess is defined as the ratio between the observed flux and the extrapolated flux, and the results are shown in Column 4 of Table 4. We detected a millimeter excess, a key signature of the RIAF, in three sources: IC 310, NGC 1277, and NGC 5846. With this information, we can quantitatively compare the significance of the millimeter excess in these sources.
As Figure 4 illustrates, NGC 224, NGC 524, and NGC 1600 exhibit an inverted spectral index, which is in good agreement with that expected with the graybody emission associated with dust emission. On the other hand, various spectral indices for the other sources suggest a different origin for the millimeter and submillimeter flux.

Empirical Classifications of LLAGNs
It has been considered that the total flux density observed at the radio band is associated with the jet, while those at FIR originated in the dust. By comparing the observed SMA flux with respect to the measurements at radio and FIR bands, we are able to investigate the origin of the emission at millimeter/ submillimeter wavelength. As a result, based on the SED characteristics, the observed SMA sources can be categorized mainly into four groups, and we discuss the property of each source that we detected at both 230 and 345 GHz in the following subsection. We show the SED diagram of each representative of the four types in Figure 5 and the associated detailed classifications in Table 4.
The second category contains the sources where the SMA flux cannot be explained by extrapolation from the radio components at centimeter wavelengths or the FIR extrapolations. Thus, the sources in the group indicate the presence of the submillimeter bump, and we need an additional feature to explain the observed SMA flux. One possible explanation is that the components originated from the central engine, such as accretion flow onto the central SMBHs. IC 310, NGC 1277, and NGC5846 are good examples of this group.
The third category includes the sources where the SMA flux density is significantly larger than the extrapolation from the Table 4 Empirical Categories, Fitted Mass Accretion Rate, and Fitted Spectral Index (1) Source name (2) Fitted spectral index for sources at 230 and 345 GHz simultaneously observed. See Section 3.2 for a detailed discussion.
(3) The spectral index is fitted using high-angular resolution measurements, and the values in parentheses indicate the spectral index obtained using low-angular resolution measurements. (4) The estimated 230 GHz flux excess from the SMA can be calculated by comparing the observed value to the extrapolated line. This is done by taking the ratio between the two values. The value outside of the parentheses is compared to a larger beam, while the value inside the parentheses is compared to a higher angular resolution spectral fit. (5) Fractional variability with constraints for sources with multiple observations. Refer to Section 3.1 for the formula details. (6) The mass accretion rate in terms of the ratio to Eddington accretion rate. (7) The fitted mass accretion rate in units of solar mass per year. (8) Bolometric luminosity evaluated via X-ray measurements. (9) The empirical classified categories provided in this paper. See Section 4 for a detailed discussion.
FIR flux measurements. At the same time, we cannot see the peak flux in the FIR band. In this case, the peak frequency in FIR could be lower than the observed lowest frequency of the millimeter wavelengths. If this is the case, the estimations from the extrapolation of the FIR data set could be the lower limit. Therefore, there is a possibility that the observed SMA total flux density could be explained by the contribution from the dust emission. For the sources in this group, we note that because of the limited measurements at the radio band for some sources, it is not clear whether the extensions can explain the SMA flux from the radio or not. NGC 4751 belongs to this group. The fourth category is a similar classification to the third category, while the SMA flux density for the sources in the class is significantly smaller than the extrapolation from the FIR flux measurements. Therefore, we think the SMA flux is more likely associated with dust than jet emissions.
We classify the sources via multiple factors, such as the millimeter and submillimeter flux relationships to the FIR or radio. While the fitted spectral index among millimeter and submillimeter ranges provided valuable information for the classification, we were able to obtain it only for a limited number of sources. Based on the samples we found, Type Is possesses a steep spectral index, while Type IVs have a positive one. Type II is slightly positive, but we only have one sample for this type.
In the context of radio-quiet AGN, such as Seyfert galaxies, it is generally accepted that these host a standard thin disk (Fausnaugh et al. 2018), which stands in contrast to the hot accretion flow. Alternative scenarios have been suggested to explain the millimeter flux, such as the corona wind associated  Table 2. with the accretion disk or star-forming cores (see, e.g., Panessa et al. 2019). Our sample, which includes NGC 5252 and Cygnus A is particularly interesting in this regard. These cases could be explained by the aforementioned alternative scenarios.
We perform the dust SED fits as in Matsushita et al. (2009) for obtaining the dust temperature parameters for the sources classified as the third and fourth categories The dust SED is assumed to have the form of S(ν) ∝ ò(ν)B(ν, T), where ò(ν) is the emissivity function, and B(ν, T) is the Plank function with given frequency ν and temperature T. ò(ν) takes the form of ( ) , where ν c is the critical frequency when the opacity becomes 1. The fitting utilized our SMA and FIR data, which are listed in the caption of Figure 3. The average dust  , 4, 5, 8, 16, 32, 64, 128, 256, 512, 1024, and 2048 times the associated image noise. The star symbol for each image shows the central AGN. The images are merged from the upper and lower sidebands to enhance the sensitivity. The image statistics are referred to in Table 3. temperature of the third and fourth categories is within the range of 15-25 K. These values are consistent with the lower end of the dust temperature in elliptical galaxies (Kokusho et al. 2019). The fitting results and the classifications are discussed in Section 4.3 and summarized together with the fitted ADAF mass accretion rate in Table 4.

Inferences of Mass Accretion Rate via the Semi-analytical Model of ADAF
Modeling the SED is one of the powerful tools in order to probe the mass accretion process onto the SMBH in the RIAF state. It has been demonstrated toward our galactic center for decades (Event Horizon Telescope Collaboration 2022), and similar approaches are considered for nearby LLAGNs (e.g., M 87). Based on the previous studies (Event Horizon Telescope Collaboration 2019), it has been suggested that the dominant emitting component can be changed from jet to accretion flow around millimeter/submillimeter wavelength, and observations in this range are extremely important to study the jet-accretion flow dichotomy, while it is an uncharted frequency range and required detailed studies.
To model the SED of LLAGNs, we adopt the analytical model of ADAF (Mahadevan 1996). This model depicts the energy distribution of the central engine for the low mass accretion rate genre in the Newtonian limit. The mass accretion rate of ADAF is typical 10  Brown et al. 2011;IC 310: Condon et al. 2002;Kadler et al. 2012;NGC 224: Condon et al. 2002;Stil et al. 2008;NGC 524: Nagar et al. 2005;Capetti et al. 2009;Brown et al. 2011;NGC 1600: Cohen et al. 2007Brown et al. 2011;NGC 1399: Shurkin et al. 2007NGC 1277: Park et al. 2017NGC 4526: Vollmer et al. 2004;Brown et al. 2011;NGC 4649: Vollmer et al. 2004;Dunn et al. 2010;Brown et al. 2011;IC 4296: Ricci et al. 2005;Murphy et al. 2009;Brown et al. 2011;NGC 5252: Edelson & Edelson 1987; PGC 049940: Allison et al. 2014;NGC 5846: Filho et al. 2004;Nagar et al. 2005;Brown et al. 2011 = . The dominant energy source is synchrotron radiation from radio to submillimeter range, inverse Compton for the submillimeter to FIR, and freefree radiation beyond FIR frequencies. It takes the following parameters to perform the calculation: the black hole mass, mass accretion rate, viscosity parameter α, and plasma β as the primary model parameters. Once the parameters are given, the model calculates the electron temperature at the first stage via equipartition assumptions. The turnover frequencies are derived afterward together with the full SED profile. Since we are conducting in the power domain, the distance to the central engine is also centric for the calculations. We used the parameters described in Table 1 for the black hole mass and distances. The semi-analytical approach provided in Mahadevan (1996) is applied to estimate the SED of an ADAF. While the model in Mahadevan (1996) includes synchrotron, bremsstrahlung, and Compton processes and is capable of estimating the spectra from 10 8 -10 20 Hz, here we focus on the spectral range ∼10 8−13 near the peak location of the synchrotron contribution. In addition to the black hole mass and accretion rate, the other fiducial parameters for the spectral computation, such as the plasma β(= 0.5) and viscosity α (= 0.3), follow the values considered in Hirotani & Pu (2016).
Given the black hole mass together with the detected flux at 230/345 GHz, we derived the mass accretion rate. The mass accretion rate and the associated categories as mentioned in Section 4.3 are listed in Table 4. The derived mass accretion rate for our samples ranges from 10 −4 to 10 −2  M Edd except for two sources: Cygnus A and NGC 6251, which are the sources that the ADAF model gives the lower limit estimates for the mass accretion rates. The large mass accretion rate presumably originated from the overestimation of the flux at millimeter/ submillimeter bands due to the contamination from the strong jet components. In the meantime, there may be similar contamination for the other LLAGNs, and the derived mass accretion rate could be an upper bound for many cases. Nevertheless, the derived range of from 10 −4 to 10 −2  M Edd is in good agreement with the RIAF prediction. We discuss individual cases more carefully in Section 4.3.
The ADAF model treats the mass accretion rate as a constant in the spatial domain, i.e., it does not vary as a function of the distance from the SMBH. On the other hand, in the broad picture, the RIAF model embeds the spatial dependent mass accretion rate by comprising the convection and wind, e.g., the ADIOS model. Thus, we should treat the mass accretion rate constrained by the ADAF model as an upper bound since the other effects usually decrease it.
In the ADAF model, the millimeter range is expected to be dominated by thermal synchrotron radiation If the millimeter frequency range is less than the synchrotron peak frequency, the spectral index is expected to be around 0.4 and to be optically thick. On the other hand, if the millimeter frequency range is higher than the synchrotron peak frequency, the spectral index can range from −2.5 to −1. These suggest that the SED can provide valuable insights into synchrotrondominated regions of ADAF and help us better understand the physics of these environments.

Multiwavelength Comparisons
In this section, we discuss the implications of the ensemble behaviors of derived mass accretion rate values via ADAF SED together with the X-ray-based measurements. The proxy for the bolometric luminosity is the X-ray measurements in the 2-10 keV range, which dominate the bolometric luminosity. And therefore, it is presumed that these strong activities originated in the vicinity of the central engine in which the accretion process plays a critical role. By comparing the bolometric power and ADAF power defined as  = , which is derived from the SED modeling, we see that the bolometric power is, in general, smaller than that of ADAF. It implies that the output power is sustainable via the accretion process. Figure 6 shows the bolometric power to the ADAF power in our samples. The efficiency (η) is defined as the ratio between the bolometric power to the ADAF power. The average value of efficiencies in our sample is around 10 −5 . It is consistent with the expectations in the hot accretion flow scheme (Yuan & Narayan 2014). The maximum value of η from our samples is 8 × 10 −4 , and the overall values of η < 0.01%, which is consistent with the expected range for the ADAF, indicating that most of the accreting energy is just advecting onto the SMBH.
Next, we compare the jet power with respect to the accretion power. Typically, jet power is estimated by PdV works based on the jet cavity seen in X-ray imaging. P is the pressure and dV is the variation of the observed volume. However, as lacking X-ray images with high resolution for most of our samples, the jet power is difficult to acquire. Alternatively, we can estimate the power of the jet cavity using radio flux as a proxy. To do this, we use the empirical scaling relations between jet power (p jet ) and synchrotron power (p radio ) derived in Cavagnolo et al. (2010): ( ) p p 8.6 10 jet radio 40 0.7 =´. We convert the observable radio flux into jet power using this relation, with the criterion being the radio flux observed at the minimum frequencies, provided they are greater than 1.5 GHz. The efficiencies (p jet /p ADAF ) measured in the above approximation are generally lower than 1%, which is again consistent with the expectation of low efficiencies from the hot accretion flow.

Angular Resolution Differences
While all the targets are observed as pointlike sources for each instrument at various frequencies, the angular resolution of each instrument we used is different. It may lead to a mismatch of the corresponding sources and could impact our interpretations of SEDs. For instance, we use data obtained by the Herschel Telescope for FIR, which provides a typical angular resolution of 5″ to 37″ (Pilbratt et al. 2010). If the source lacks the Herschel data, we use the Infrared Astronomical Satellite (IRAS; Neugebauer et al. 1984) instead. The typical angular resolution ranges from 30″ and 2′ for 12 μm and 100 μm, respectively. Therefore, the SMA beam is smaller, and we may overestimate FIR flux. On the other hand, we typically use a Very Long Baseline Array (VLBA) and Very Large Array (VLA) for the centimeter to millimeter wavelengths data set, in which the typical angular resolution scales around a few milliarcseconds and a few arcseconds, respectively. The SMA beam is comparable with the VLA; hence, the extrapolation Figure 5. SEDs together with the ADAF and dust fit of four galaxies, including IC 1459, IC 310, NGC 4751, and NGC 1600, represent Type I-IV AGNs, respectively. The y-axis displays luminosity, and the x-axis represents the observed frequency. Each subplot shows (1) The components of the ADAF spectrum with a given central SMBH mass, where for each panel from top to bottom, the blue profiles are fit with varying mass accretion rates of 10 −2 , 10 −3 , and 10 −4 times the Eddington accretion rate.
(2) The light green distribution represents the dust fit utilizing the FIR and SMA measurements. The first category includes the sources with flux densities observed by SMA exceeding the fitted spectral index extrapolated from FIR flux measurements. In the meantime, the SMA flux densities are in good agreement with the extrapolated flux using the total flux density obtained with VLA (and VLBI) observations at centimeter wavelength. It suggests the SMA flux is not associated with the dust emission, but it is reasonable to assume it is related to the jet emission. This group contains well-known nearby AGNs with prominent jets, such as IC 1459, NGC 6251, and Cygnus A.
does not affect much in terms of angular resolution. However, the SMA beam is almost three orders of magnitude larger than the VLBIs.
In the following, the beam mismatch is considered in the prevailing circumstances when we discuss the contribution from the extra components in comparison with centimeter-tomillimeter/submillimeter measurements. On the other hand, we discuss the property of individual sources by assuming the sources have pointlike structures in comparison with FIR-tomillimeter/submillimeter measurements. Note that this situation will be improved with better SED modeling with the matched and high-angular resolution for each instrument with future observations.

Individual Source Discussions
In this subsection, we discuss the sources with 230 and/or 345 GHz detections. The category mentioned in Section 4.1 along with the fitted mass accretion rate is briefly summarized.

Type I: IC 1459
IC 1459 is one of the nearby galaxies hosting a considerable supermassive black hole. With previous VLBI observations, this source is categorized as the gigahertz peaked spectrum source, which is usually considered a young radio source that could be evolved into FR II (Tingay & Edwards 2015). IC 1459 exhibits that the SMA flux positively deviates from the extrapolations from the FIR measurements by assuming graybody radiation (See Figure 3) and therefore excludes the possibility that the galactic dust emission dominates at millimeter/submillimeter bands.
In the lower-frequencies radio band, there is a discrepancy between the trend of the total flux measured at the arcsecond scale by the ATCA and the trend of the milliarcsecond scale by the VLBI. As seen in the sub-figure for the IC 1459 in Figure 3. While the ATCA image shows a point-source-like structure, there should be a substructure within the arcsecond scale. If we compare the SMA flux with the flux extrapolated from the measurements of the ATCA, it shows a good agreement with each other. It may suggest that the extended components, which are resolved with VLBI, can be responsible for the differences between the SMA flux density and estimated flux density with VLBI measurements by assuming a single powerlaw spectrum. In this case, the millimeter/submillimeter flux is more likely to be the jetted components. Note that the SMA 230 and 345 GHz spectral index shows a different tendency with the VLA/ATCA/VLBI measurements, a nearly flat spectrum.
Alternatively, if we assume that the SMA flux is associated with the ADAF, the upper bound of the mass accretion rate of IC 1459 can be estimated to be 5.7 × 10 −2 M e yr −1 by adopting the semi-analytical ADAF model. This result is consistent with the prediction for LLAGNs and is in good agreement with Balmaverde et al. (2008), who reported a similar mass accretion rate of 5.3 × 10 −2 M e yr −1 with the Chandra data. The submillimeter bump behavior is also well explained by an ADAF model.

Type I: NGC 6251
NGC 6251 is a BL Lac object with Doppler-boosted radio emission (Chiaberge et al. 2003). It is worth noting that sources classified as BL LAC/Seyfert may not exhibit a big blue bump at optical wavelengths, and this classification may have limitations that warrant further discussion. Additionally, the outflow inclination (as described by BL LAC) and the  Table 4. See Section 4.2.2 for further discussion. bolometric behaviors (as described by LLAGN) are not mutually exclusive properties. Understanding the diverse properties of BL LAC/Seyfert sources is of great importance and merits further investigation into their underlying physics.
NGC 6251 is one of the nearby giant elliptical galaxies, which has an exceptionally straight and long jet from sub-pc to sub-Mpc scale (Tseng et al. 2016). Since the FIR observations do not show a clear peak, one may assume the dust emission is responsible for the total flux density at 345 GHz. However, 230 GHz is less likely to share the same origination as 345 GHz in this case due to the different spectral slopes. Consequently, we think the origin of the flux at millimeter/submillimeter from the dust is unlikely, at least for the 230 GHz.
On the other hand, there is no discrepancy between VLBI and VLA core flux measurements, and the millimeter/ submillimeter flux follows the trend of the radio observations in the SED. The spectra derived at these frequencies show a single power-law nature with negative spectral indices (see Figure 3). It indicates that optically thin synchrotron radiation associated with the jet components is probably responsible for the total flux density in the submillimeter regime. While the SMA flux is interpretable with the synchronous jet, it is also utilized to derive the upper bound for the possible accretion flow composite. The upper limit for the mass accretion rate provided by the ADAF model gives the value of 1.3 × 10 −1 M e yr −1 .

Type I: IC 4296
IC 4296 is one of the brightness sources in our sample in terms of the 230 GHz flux. It is a nearby elliptical galaxy hosting a typical FR I morphology (Grossová et al. 2019). The observed SMA total flux is substantially lower than the total flux density extrapolated from the FIR measurements via a graybody (Golombek et al. 1988), indicating that contribution from the dust would not be sufficient.
The radio core and symmetric two-sided jets are clearly detected at arcsecond scale resolution with previous centimeterradio observations (Grossová et al. 2019). The observed SMA total flux density is slightly lower than the estimated flux density with extrapolation using the multifrequency radio measurements toward the core. And therefore, it is natural to assume that a significant contribution to the SMA total flux density comes from the jet. Similar to the other sources, conditional on the ADAF model fits, IC 4296 can derive the upper bound of the mass accretion rate, and it is 2.9 × 10 −2 M e yr −1 . This value is smaller than the Bondi accretion rate of 3 M e yr −1 reported by Balmaverde et al. (2008) using Chandra data, which measures a more outer region. This discrepancy can potentially be explained by the decreasing mass accretion rate toward the inner region, as predicted by the RIAF model. Therefore, while this does not necessarily contradict the classification of IC 4296 as an LLAGN, it does suggest that a more generalized RIAF model may be necessary for this source, given that the outflow or dissipation of the mass accretion rate could be critical for explaining the discrepancies of the mass accretion rate.

Type I: NGC 5252
NGC 5252 is classified as a Seyfert galaxy that displays the presence of an Extended Narrow Line Region. Weak emissions associated with the optical nuclei have been revealed through VLA observations. However, the expected contribution from a dust component, based on graybody emission extrapolated from IRAS observations (Meléndez et al. 2014), is lower than the observed millimeter flux of NGC 5252 by an order of magnitude. Additionally, the millimeter flux has an excess of 40% compared to the expected value from a centimeter trend. Therefore, we speculate that the detected SMA total flux is associated with the jet emission and has a similar origin to the radio emission. As a result, we do not provide a mass accretion rate fit for this source due to its distinct behavior.

Type I: NGC 1052
NGC 1052 is one of the giant elliptical galaxies, also classified as a low-ionization nuclear emission-line region (LINER) galaxy (Heckman 1980). Previous VLBI observations unveiled the presence of the radio jet and counter-jet, together with the presence of ionized torus responsible for free-free absorption at centimeter wavelengths (Kameno et al. 2001). This source is known as the gigahertz peaked spectrum source, which typically has a peaked flux around 1-10 GHz in its spectrum.
The power converted from the SMA flux of NGC 1052 exceeds the expectations from the FIR graybody extrapolations by order of magnitudes. And therefore, dust origination is highly unlikely. This agrees with the recent studies based on VLTI/MIDI measurements NGC 4649, also known as M 60, is located in the Virgo Cluster. Its central black hole poses a large mass of 2.1 ×10 9 M e . We detect a weak flux density at 230 GHz with ∼10 mJy with the SMA. Given the features that the millimeter flux is greater than the expectations of radio trends and lower than that of FIR (see Figure 3), it implies that the flux can be originated from either dust and/or a signature of a submillimeter bump. The SMA flux provides an upper limit on the mass accretion rate, which we estimate to be ∼10 −2 M e yr −1 using ADAF models. This value is consistent with the Bondi accretion rate reported in Vattakunnel et al. (2010), where the Chandra measurements yielded a value of 1.1 × 10 −2 M e yr −1 .

Type I: Cygnus A
This source is considered as a representative source for FR II-type radio galaxies, as established by Fanaroff & Riley (1974). As such, there is a wealth of information available from previous observations. Lo et al. (2021) detected a low polarization fraction at the 230 GHz of Cygnus A, and interpreted it with the front and concentrated magnetized plasma toward the central region. One of the reasonable explanations is due to the magnetized plasma associated with the ADAF/RIAF. And therefore, it is worth noting that Cygnus A, while being a high excitation FR II source, a low polarization nature at the millimeter/submillimeter range makes it debatable whether it can be classified as an LLAGN. However, its powerful jet activities, P jet ∼ 3.9 × 10 45 erg s −1 (McNamara et al. 2011), make it a source worthy of further investigation.
The SMA flux densities exceed the fitted spectral index extrapolated from FIR flux measured with Hershel. It indicates that the SMA flux is not originated from dust associated with the host galaxy (See Figure 3 for the references).
If we compare the extrapolated flux using the total flux density obtained with VLBI observations at centimeter wavelength, the SMA flux densities at 230 and 345 GHz are an order of magnitudes higher than the extrapolated flux. However, the SMA flux is in good agreement with that estimated total flux density extrapolated with VLA measurements. These imply that millimeter flux can be dominated by the contamination from the radio jet, and Cygnus A can be classified as type I in our category.
Therefore, we speculate that there is extended jet emission, which resolved out at milliarcsecond resolution of centimeter-VLBI but contributes to the total flux density at 230/345 GHz at arcsecond resolution. In this case, the core flux poses a single power law up to 230/345 GHz, indicating less evidence for the synchrotron cooling beyond 6 mas Boccardi et al. (2016).
Alternatively, there is a possibility that a different component is only visible at 230/345 GHz. It is evident for SED in Sgr A * , and referred to as millimeter/submillimeter bump (Dexter et al. 2010). One of the possible origins is ADAF, as indicated by some of the millimeter flux surveys toward LLAGNs (Doi et al. 2005(Doi et al. , 2011. If we adopt a semi-analytic ADAF model to our measurements as the first-order approximation, the mass accretion rate of Cygnus A is estimated to be greater than 5.9 × 10 −1 M e yr −1 . It is slightly above the Eddington limit for the LLAGN, while it can sufficiently explain the observed flux. Hence, it could be that Cygnus A is less likely to be classified as LLAGN, and/or the adopted model is insufficient.

Type I: PGC 049940
PGC 049940 is identified as the brightest cluster galaxy of the Abell 1836, displaying FR II-type morphology with the bright central core at 1.6 GHz (Liu & Zhang 2002;Stawarz et al. 2014). The millimeter flux follows the extrapolation from the low frequencies and is more prominent than that from the FIR. The radio jet and PGC 049940 can be classified as type I in our category. Alternatively, the mass accretion rate estimation given by the ADAF model is an upper limit considering the jet contamination, which gives the 8.6 × 10 −2 M e yr −1 . Our analysis is consistent with the mass accretion rate estimations in Stawarz et al. (2014), which are analyzed via optical measurements assuming the radiative efficiency (the similar definition of η in Section 4.2.2) of 1%.

Type I: NGC 1399
NGC 1399 is the central galaxy of the Fornax cluster and shows a weak radio core at 1.4 GHz (Jones et al. 1997). The millimeter flux is at least two orders of magnitude higher than the dust profile fitted from the FIR as shown in Figure 3, suggesting that dust emission is negligible. However, due to the lack of radio observations at other frequencies, it is difficult to determine whether the radio emission is affected by extended features such as jets. NGC 1399 can be classified as type I in our category. The mass accretion rate constrained by the ADAF model is 2.8 × 10 −2 M e yr −1 , which is lower than that of PGC 049940 and implies a lower radiative efficiency.

Type II: IC 310
IC 310 is a head-tail radio galaxy (Gendron-Marsolais et al. 2020) located in the Perseus galaxy cluster. High energy γ-ray emission was detected in IC 310 (Aleksić et al. 2014). black hole magnetospheres (Hirotani & Pu 2016) together with the ADAF are used as one of the models to explain the high energy radiation associated with it.
First, the SMA flux density is higher than the estimated total flux density by extrapolating from the FIR flux at 100 μm with dusty graybody emission. Note that we do not see the peak in the SED within the FIR regime associated with the dust. And therefore, peak frequency could be lower than the observed ranges. In this case, there is a possibility that the observed SMA flux at 345 GHz is associated with the dust. However, even in this case, the dust cannot explain SMA flux at 230 GHz. Thus, we presume the SMA flux is not originated from the dusty components.
Secondary, very different from Cygnus A and IC 1459, there is no difference in the radio flux densities observed at arcsecond and milliarcsecond scale resolutions. If we simply extrapolate the radio flux up to the millimeter/submillimeter regime, the observed SMA flux at 230 and 345 GHz are higher by factors of 5.1 and 7.6, respectively. Therefore, it requires an additional component to correspond to the emission and clearly indicates the existence of the millimeter/submillimeter bump.
If we adopt a simple ADAF model, the mass accretion rate is estimated to be 1.4 × 10 −1 M e yr −1 , respectively. It is in good agreement with the prediction for LLAGNs. Therefore, we think this is a genuine case, we directly observed the RIAF components.

Type II: NGC 1277
NGC 1277 is one of the galaxies, which are suggested to host ultramassive BHs of the order of 10 10 M e at the center. It has been shown that NGC 1277 has a steep spectrum from 1.5-5 GHz and suggests an inverted spectrum toward the millimeter/submillimeter with the PdBI measurements at 115 GHz (Sijbring 1993;Scharwächter et al. 2016). Our independent measurements with the SMA at 230 GHz revealed an inverted spectrum from 115-230 GHz. It demonstrates that NGC 1277, similar to IC 310, potentially harbors a submillimeter bump in the SED. The possibility that the millimeter/ submillimeter emission originates from the dust emission is ruled out by utilizing the Herschel measurements as claimed in Scharwächter et al. (2016). The ADAF fit of the mass accretion rate gives the value of 3.7 × 10 −1 M e yr −1 .

Type II: NGC 5846
NGC 5846 is also an elliptical galaxy located in the Virgo Cluster. Since the FIR measurements indicate the turnover feature around millimeter to the smaller wavelengths, it is natural to assume the graybody radiation from FIR in order to derive the upper bound for the contributions of the dust. The total flux measured via SMA is larger than the extrapolated total flux density from the FIR measurements assuming dust emission. Hence, we think the contribution from the dust to the observed SMA total flux density is not significant.
Similar to IC 310, the estimated mass accretion rate of 1.8 × 10 −2 M e yr −1 for NGC 5846 based on the semianalytical ADAF model is consistent with the prediction for LLAGNs and in good agreement with Balmaverde et al. (2008), who reported a similar Bondi accretion rate of 1.3 × 10 −2 M e yr −1 with the Chandra data. The excess of the total flux density at 230 GHz compared to the extrapolated value from the radio VLA/VLBI observations suggests the presence of the submillimeter bump, and its origin is neither dust nor jet. Instead, it could be associated with the accretion flow, supporting the possibility of the ADAF model.

Type III: NGC 4751
NGC 4751 is weakly detected at 230 GHz as of ∼13 mJy with the SMA. The millimeter flux is in good agreement with the extrapolation from the FIR (see Figure 3). The fitted dust temperature as described in Section 4.1 is 25 K. Together with lack of the detection at centimeter wavelength, we cannot role out the possibility that the origin of the emission is associated with dust. In this case, we can use the millimeter flux as the upper limit coming from the accretion flow, and the mass accretion rate hinted by the ADAF model gives the value of 3.1 × 10 −2 M e yr −1 as an upper bound.

Type IV: NGC 4526
NGC 4526 is a spiral galaxy located in the Virgo Cluster. VLA measurements reveal that it is a point source with an extended structure along the east-west direction (Nyland et al. 2017). Our SMA study shows that NGC 4526 is unresolved at 230 GHz and weakly detected at ∼16 mJy. It is possibly explained by the dust emission as the hints from the FIR flux density extension, and the fitted dust temperature as described in Section 4.1 is 20 K. Hence, NGC 4526 is labeled as Category IV in our analysis, and the upper bound of the mass accretion rate constrained by the ADAF model gives ∼10 −2 M e yr −1 . The value from our ADAF model is approximately an order of magnitude larger than the mass accretion rate at the Bondi radius estimated by Vattakunnel et al. (2010) using Chandra measurements, which was reported to be 1.1 × 10 −3 M e yr −1 . This discrepancy may be explained by the possibility of contamination of the millimeter flux by the dust components.

Type IV: NGC 1600
NGC 1600 is one of the nearby elliptical galaxies, which possesses a very large SMBH up to 1.7 × 10 10 M e based on stellar dynamics (Thomas et al. 2016). Our SMA observations provide the 10.5 and 58.8 mJy at 230 and 345 GHz toward NGC 1600, respectively. The positive spectral index implies that the source behaves optically thin in the millimeter/ submillimeter bands range, and the fitted dust temperature is 15 K, both suggesting that the source is dust-dominated. With our SMA observations, we derived an upper limit for the mass accretion rate to be 3.7 × 10 −1 M e yr −1 by assuming the ADAF model.

Type IV: NGC 224
NGC 224 is one of the nearest galaxies, known as M31, and the Andromeda galaxy. It is known that this galaxy hosts the most quiescent SMBH at its center, similar to the one at our galactic center (Sgr A * ). Based on the Jansky VLA observations, Yang et al. (2017) reported that the total flux density of NGC 224 is around 10 μJy between 10 and 20 GHz with a spectral index of −0.5, which suggests optically thin synchrotron radiation. Smith et al. (2021) also show that the flux density is less than 10 mJy beam −1 around its central region from the James Clerk Maxwell Telescope at 850 μm. Our SMA observations toward M 31 confirm the above prior studies and detect this source with the flux of ∼20 and ∼50 mJy at 230 and 345 GHz, respectively. Similar to NGC 1600, the positive spectral index among millimeter and submillimeter flux and the fitted dust temperature of 20 K, both indicate that the source is dust-dominated in this frequency range. Nevertheless, conditional on the ADAF model, the mass accretion rate is provided as the upper limit of 3.1 × 10 −3 M e yr −1 .

Type IV: NGC 524
NGC 524 is another nearby galaxy with a compact radio emission at its center. The nature of the SED is similar to that of NGC 1600. While it has been detected with VLA (Nagar et al. 2005) and VLBA (Filho et al. 2004), the detected flux is relatively weak. We classify this source as Type IV as it carries a positive spectral index at the millimeter and submillimeter range and a fitted dust temperature of 15 K, regarded as dustdominated. The flux densities of NGC 524 are 18 and 48 mJy at 230 and 345 GHz. Conditional on the ADAF model, the mass accretion rate is provided as the upper limit of 1.8 × 10 −3 M e yr −1 .

Flux Densities and Diameters of Black Shadow
The flux densities obtained via the 230 and 345 GHz wave bands with SMA provide a pathfinder to investigate whether it is viable to follow up the sources with instruments with higher angular resolution. We show the observed SMA fluxes at 230/ 345 GHz and the expected diameter of black hole shadow in Figure 7. Here we presumed the expected diameter of the black hole shadow to be five times the angular size of the Schwarzschild radius, which is in reasonable agreement with the anticipated photon ring size (Johannsen & Psaltis 2010).
NGC 1600 is expected to possess the apparently most significant black hole shadow in our sample, and the flux densities of ∼50 mJy at 345 GHz can be the potential target of the future EHT. However, the positive spectral properties at the millimeter/submillimeter bands suggest that we cannot rule out the possibility that dusty components dominated the flux. Hence it may not be resolvable with microarcsecond scale instruments. NGC 224 shares a similar argument with NGC 1600.
IC 1459 has the next largest black hole shadow, and its size is expected to be 11.4 mas. Notably, SMA detected a point source with the flux density of ∼200 mJy, which is within the sensitivity range of the current EHT. We anticipate that the flux density origin is associated with the innermost jet, similar to the M87. Therefore, IC 1459 is one of the prospective sources for future EHT observations.
For the sources whose shadow size is lower than 10 mas, a systematic study of LLAGN on the nature of accretion flow and/or innermost jet at the black hole shadow scale measurements is possible, while it is less likely to observe the black hole shadow with the current instruments.
In Section 4.2.2, we discussed the efficiencies between accretion power and bolometric power, which may be overestimated due to the use of the total flux density from the SMA observations in the semi-analytical model. This is because the total flux density may include contamination from the extended jet components. For instance, the ratio between EHT flux to the SMA or ALMA core flux of M87 (Kuo et al. 2014;Bower et al. 2018; Event Horizon Telescope Collaboration 2019) is around 0.66 Jy/1.32 Jy ; 50%, indicating that the SMA beam may be affected by the jet components. This claim is supported by the observation that the flux decays when the synthetic beam shrinks to the region closer to the central engine. In order to obtain more accurate estimates, higher angular resolution observations are necessary. Future millimeter/submillimeter VLBI observations are expected to provide better constraints on the efficiency of the accretion process and the bolometric power.

Summary
We utilized the SMA to conduct millimeter and submillimeter continuum observations toward 36 nearby sources that potentially have large black hole shadows. The summary of our results is as follows.
We have 17 and eight sources detected at 230 and 345 GHz, respectively. For the sources detected at 345 GHz, all are simultaneously detected at 230 GHz as well. The sources are categorized into four subsets by the multiwavelength comparisons, essentially archived low-frequency radio and midinfrared measurements. The detected flux of Type I is considered to be explainable by the nonthermal extended jet emission of the LLAGN and is not by the dust. There are nine sources in this category. The millimeter/submillimeter bump feature is considered to be prominent for Type II sources since neither nonthermal jet nor dust can explain the flux. In this category, it is feasible to have the hot accretion flow around the central engine. There are three sources in Type II: IC 310, NGC 1277, and NGC 5846. The remaining Types III and IV are considered to have dust-contributed nature in their system. There are one or three sources for each group. By utilizing the semi-analytical ADAF model and the exogenous black hole mass and distance, we fit the mass accretion rate to each detected source. The fitted mass accretion rate is, in general, less than 10 −2  M Edd , which is consistent with the predictions of the LLAGN property.
The largest diameter of the black hole shadow of 26 mas within our samples, NGC 1600, presumably has the dustdominated characteristic at millimeter and submillimeter wavelengths. It yields imaging of the black hole shadow toward this source less plausible. On the contrary, IC 1459, which possesses the feature of being optically thin at these wavelengths and considerable sufficient flux density, is the alleged candidate for future submillimeter VLBI measurements to conduct the imaging of the black hole shadow.
The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. W.P. L. and K.   Table 1, and is in the units of microarcsecond (μas). The blue circles denote the sources with detections at 230 GHz, while the green diamonds indicate the measurements for 345 GHz. The error bars for the 230 and 345 detecting points denote the 1σ uncertainties. The red squares are the sources without detections at 230 GHz, and the value for the flux densities is the 3σ level as the detection upper bound.