THE EFFECTS OF ACCRETION FLOW DYNAMICS ON THE BLACK HOLE SHADOW OF SAGITTARIUS A*

Very long baseline interferometry (VLBI) observations of the nearest supermassive black hole, Sagittarius A* (Sgr A*), have provided valuable information concerning the emission region within a few Schwarzschild radii of the black hole (e.g., Doeleman et al. 2008; Fish et al. 2011, 2016; Johnson et al. 2015). The extremely dim luminosity of Sgr A* ($\sim {10}^{-9}{L}_{\mathrm{Edd}}$, where ${L}_{\mathrm{Edd}}$ is the Eddington luminosity) suggests that the accretion flow around the black hole is in the radiatively inefficient accretion flow (RIAF) regime (e.g., Narayan et al. 1995; Manmoto et al. 1997; Özel et al. 2000; Yuan et al. 2003). Compared with a cold, geometrically thin, Keplerian rotating disk (Novikov & Thorne 1973; Page & Thorne 1974), the radiative cooling timescale is much longer than the accretion timescale for an RIAF, resulting in a hot, geometrically thick flow with a sub-Keplerian rotation (e.g., Ichimaru 1977; Narayan & Yi 1994; Narayan et al. 1997; Quataert & Gruzinov 2000).

With a mass of $\sim 4.3\times {10}^{6}$ M_⊙ and at a distance of ∼8.3 kpc (Ghez et al. 2008; Gillessen et al. 2009a, 2009b; Chatzopoulos et al. 2015), Sgr A* is expected to cast a shadow with an angular diameter of ∼50 μas, as predicted by general relativity (Bardeen 1973; Falcke et al. 2000). With ultra-high angular resolution (up to ∼20 μas) and sufficient sensitivity, the forthcoming millimeter/sub-millimeter VLBI observations at $\lambda \lesssim 1.3\,\mathrm{mm}$ ($\nu \gtrsim 230$ GHz), from the Event Horizon Telescope (EHT), should obtain the first image of such a shadow region (Fish et al. 2014; Lu et al. 2016). The emission from the accretion flow in the vicinity of the black hole, before reaching a distant observer, will experience considerable energy shifts due to both the motions of the emitting fluid medium and the strong gravity near the black hole (Broderick & Blandford 2004; Younsi et al. 2012). As such, the observed spectrum, image, and visibility are important indicators to determine the flow dynamics when the accretion flow within several Schwarzschild radii of the black hole is optically thin enough to be observed (at millimeter/sub-millimeter wavelengths in the case of Sgr A*), and are therefore important topics to be explored.

In the pioneering work of Falcke et al. (2000), the authors investigate the observational appearance of the accretion flow and the shadow of Sgr A* at millimeter/sub-millimeter wavelengths. Two kinds of accretion flow dynamics were considered therein: (i) plasma with a free-fall motion, and (ii) plasma modeled as a rigidly rotating shell in Keplerian motion. In the latter case, the flow was assumed to be in Keplerian motion outside the ISCO. Within the ISCO, the flow follows geodesics, with energy and angular momentum specified at the ISCO boundary (Cunningham 1975). Subsequently, by adopting the rigidly rotating Keplerian shell model, Broderick & Loeb (2006) further considered separate populations of thermal and non-thermal electrons, based on a radial power-law-dependence of vertically averaged density and temperature profiles found in Yuan et al. (2003). The resulting spectrum is reasonably consistent with the observational data. Although GRMHD simulations provide further detailed, time-dependent modeling of the accretion flow (e.g., Noble et al. 2007; Mościbrodzka et al. 2009; Dexter et al. 2010; Dolence et al. 2012; Shcherbakov et al. 2012; Shcherbakov & McKinney 2013; Chan et al. 2015a, 2015b; Gold et al. 2016), (semi-)analytic models enable an efficient and flexible means to survey the vast parameter space, constraining the black hole spin, inclination angle, and position angle of Sgr A* (e.g., Broderick et al. 2011a, 2011b, 2016).

Motivated by the sub-Keplerian nature of RIAFs, in this paper, we investigate the resulting spectra, images, and visibilities of a black hole surrounded by accretion flows with differing dynamics. The purpose of this paper is to demonstrate the signature that fluid dynamics alone can produce, by employing a semi-analytic approach. Such an approach will be useful for parameter survey studies for EHT observations and comparison with GRMHD simulation results, as will be applied for our future studies. Of course, these observables are also germane to other factors related to the underlying radiative processes, such as magnetic field configuration and the spatial and energy distribution of electrons in the flow.

This paper is organized as follows. In Section 2, we describe how we model the dynamic of a flow, including a flow with sub-Keplerian rotation mimicking an RIAF. In Section 3, we compare and discuss the spectra and images of a flow subjected to different dynamics. The visibility amplitudes and EHT observations for Sgr A* in 2009 and 2013 are also compared. Finally, a summary and a discussion of future perspectives are presented in Section 4.

In order to investigate the sub-Keplerian nature of an RIAF, we consider flow models with dynamics consisting of a combination of stationary free-fall and Keplerian motions (c.f., Takahashi & Mineshige 2011).

Adopting a $[-+++]$ signature of the background Kerr metric, and $c=G=M=1$ as the normalization of the four-velocity ${u}^{\mu }{u}_{\mu }=-1$, gives

$$\begin{eqnarray}&&{({u}^{t})}^{2}=\displaystyle \frac{1+{g}_{{rr}}{({u}^{r})}^{2}}{{K}_{0}},\end{eqnarray} \tag{ 1 }$$

where ${K}_{0}=-({g}_{{tt}}+{{\rm{\Omega }}}^{2}{g}_{\phi \phi }+2{\rm{\Omega }}{g}_{t\phi })$, ${\rm{\Omega }}={u}^{\phi }/{u}^{t}$ is the the orbital frequency, and ${u}^{\theta }$ is assumed to be small enough that it is negligible. We can now construct a flow model with a four-velocity $({u}^{t},{u}^{r},{u}^{\theta }=0,{u}^{\phi }={\rm{\Omega }}{u}^{t})$ for a given Ω through the following steps: (i) determine u^r, (ii) compute u^t from Equation (1), and then (iii) calculate ${u}^{\phi }$. The only physical requirement is that ${K}_{0}\gt 0$ so that ${({u}^{t})}^{2}\gt 0$ is always true.

It is then straightforward to constrain u^r from the Keplerian value ${u}_{{\rm{K}}}^{r}(a,R=r)$ and the free-fall value ${u}_{\mathrm{ff}}^{r}(a,r,\theta )$, considering the sub-Keplerian case as a combination of the two:

$$\begin{eqnarray}&&{u}_{\mathrm{sub}-{\rm{K}}}^{r}={u}_{{\rm{K}}}^{r}+(1-\alpha )({u}_{\mathrm{ff}}^{r}-{u}_{{\rm{K}}}^{r}),\end{eqnarray} \tag{ 2 }$$

where a is the dimensionless black hole spin parameter, and $0\leqslant \alpha \leqslant 1$. For the Keplerian disk case, R is the distance to the black hole along the equatorial plane. A Keplerian rotating shell is thus described by applying R = r (and neglecting the θ-dependence).

Similarly, Ω may be constrained by the Keplerian value ${{\rm{\Omega }}}_{{\rm{K}}}$ and the free-fall value ${{\rm{\Omega }}}_{\mathrm{ff}}$, considering the sub-Keplerian case as the mixture of the two:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{sub}-{\rm{K}}}={{\rm{\Omega }}}_{{\rm{K}}}+(1-\beta )({{\rm{\Omega }}}_{\mathrm{ff}}-{{\rm{\Omega }}}_{{\rm{K}}}),\end{eqnarray} \tag{ 3 }$$

where $0\leqslant \beta \leqslant 1$. The free parameters α and β control how much the flow deviates from the Keplerian ($\alpha =\beta =1$) and free-fall motions ($\alpha =\beta =0$), specifying the radial and toroidal motion of the flow. Here, we set $\alpha =\beta =0.5$ to represent a sub-Keplerian flow. This approach is generally applicable for all black hole spins. The profiles of u^r and Ω for the cases a = 0 and a = 0.998 are shown in Figure 1. Differences between different flow dynamics become more obvious as the flows approach the black hole. This implies that observations at optically thin (i.e., millimeter/sub-millimeter) wavelengths are crucial to determine the flow dynamics and discriminate between different models.

Zoom In Zoom Out Reset image size

In what follows, we adopt the best-fit parameters (a = 0, inclination angle $i=68^\circ$) constrained by the model fitting for 1.3 mm (230 GHz) VLBI observations in Broderick et al. (2011a), which are also consistent with Broderick et al. (2016). We assume hybrid populations of thermal and non-thermal electrons that obey a relativistic Maxwellian distribution and a power-law distribution in energy, respectively. The non-thermal electrons follow a power-law distribution with a spectral index of 1.25 (which corresponds to an energy index of 3.5), and a lower cutoff Lorentz factor of 10². In addition, similarly to Broderick et al. (2011a), the spatial distributions of the temperature and density of these populations are described by

$$\begin{eqnarray}&&{n}_{e,\mathrm{th}}={n}_{e,\mathrm{th}}^{0}\ {r}^{-1.1}{e}^{-{z}^{2}/2{\rho }^{2}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{T}_{e}={T}_{e}^{0}\ {r}^{-0.84},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{n}_{e,\mathrm{nth}}={n}_{e,\mathrm{nth}}^{0}\ {r}^{-2.02}{e}^{-{z}^{2}/2{\rho }^{2}},\end{eqnarray} \tag{ 6 }$$

where $\rho =r\sin \theta$, and $z=r\cos \theta$. The normalizations of ${n}_{e,\mathrm{th}}^{0}$, T⁰_e, and ${n}_{e,\mathrm{nth}}^{0}$, which are equivalent to specifying the accretion rate, are chosen to fit the observed spectrum; therefore, they are functions of black hole spin and inclination angle. In such a description, the electron population decreases rapidly when approaching the funnel region, $\rho \to 0$.

The resulting emission is considered as synchrotron emission from the hybrid populations, with the magnetic field strength determined via, in physical units,

$$\begin{eqnarray}&&\displaystyle \frac{{B}^{2}}{8\pi }=0.1\ {n}_{e,\mathrm{th}}\displaystyle \frac{{m}_{{\rm{p}}}{c}^{2}{r}_{{\rm{g}}}}{6r},\end{eqnarray} \tag{ 7 }$$

where ${r}_{{\rm{g}}}={GM}/{c}^{2}$. Recent 1.3 mm observations with the EHT have detected linear polarization on scales of several ${r}_{{\rm{g}}}$, indicating the existence of a magnetic field structure that is predominantly stochastic, but also partially ordered on event horizon scales (Johnson et al. 2015). However, although the magnetic field configuration is associated with differential motion within the flow, via the coupling between the plasma and the magnetic field, the field configuration resulting from different flow dynamics is unclear. Here, we adopt the angle-averaged emissivity ${j}_{\mathrm{th}}$ in Mahadevan et al. (1996) for computing thermal synchrotron emission, and the absorption coefficient ${\alpha }_{\mathrm{th}}$ is obtained via Kirchoff's law. For non-thermal emission, the emissivity ${j}_{\mathrm{nth}}$ and absorption coefficient ${\alpha }_{\mathrm{nth}}$ are taken from Dexter et al. (2012) by applying a fixed pitch angle of $\pi /3$. This enables us to ensure all parameters—except the flow dynamics—remain unchanged. The resulting images and spectra of the three models, as summarized in Table 1, are calculated using the general-relativistic radiative transfer code Odyssey (Pu et al. 2016).

Table 1.  Accretion Flow Model Parameters

a = 0, $i=68^\circ$, $\alpha =\beta =0.5$, ${n}_{e,\mathrm{th}}^{0}=2.5\times {10}^{7}\ {\mathrm{cm}}^{-3}$, ${T}_{e}^{0}=1.5\times {10}^{11}\ {\rm{K}}$, ${n}_{e,\mathrm{nth}}^{0}=8\times {10}^{4}\ {\mathrm{cm}}^{-3}$
Model	ur	${\rm{\Omega }}={u}^{\phi }/{u}^{t}$
Kep ($r\gt {r}_{0}$)	${u}_{{\rm{K}}}^{r}(r)=0$	${{\rm{\Omega }}}_{{\rm{K}}}{(r)=({r}^{3/2}+a)}^{-1}$
Kep (r ≤ r0)	${u}_{{\rm{K}}}^{r}(r)=-{[2/(3{r}_{0})]}^{1/2}{[({r}_{0}/r)-1]}^{3/2}$	${{\rm{\Omega }}}_{{\rm{K}}}(r)=(\lambda +{aH}){[{r}^{2}+2r(1+H)]}^{-1}$
sub-Kep	${u}_{\mathrm{sub}-{\rm{K}}}^{r}(r,\theta )={u}_{{\rm{K}}}^{r}+(1-\alpha )({u}_{\mathrm{ff}}^{r}-{u}_{{\rm{K}}}^{r})$	${{\rm{\Omega }}}_{\mathrm{sub}-{\rm{K}}}(r,\theta )={{\rm{\Omega }}}_{{\rm{K}}}+(1-\beta )({{\rm{\Omega }}}_{\mathrm{ff}}-{{\rm{\Omega }}}_{{\rm{K}}})$
free-fall	${u}_{\mathrm{ff}}^{r}(r,\theta )=-{[2r({r}^{2}+{a}^{2})]}^{1/2}{{\rm{\Sigma }}}^{-1}$	${{\rm{\Omega }}}_{\mathrm{ff}}(r,\theta )=2{ar}{{ \mathcal A }}^{-1}$

Note: ${\rm{\Sigma }}\equiv {r}^{2}+{a}^{2}{\cos }^{2}\theta$, ${\rm{\Delta }}\equiv {r}^{2}-2r+{a}^{2}$, $\lambda \equiv ({r}_{0}^{2}-2a\sqrt{{r}_{0}}+{a}^{2})/(\sqrt{{r}_{0}^{3}}-2\sqrt{{r}_{0}}+a)$, $H\equiv (2r-a\lambda )/{\rm{\Delta }}$, ${ \mathcal A }\equiv {({r}^{2}+{a}^{2})}^{2}-{a}^{2}{\rm{\Delta }}{\sin }^{2}\theta$, where r₀ is the ISCO radius. In every model, ${u}^{\theta }=0$ for all r.

Download table as:  ASCIITypeset image

The synchrotron spectra of the three dynamical flow models are presented in Figure 2 as solid lines, a result of the competition between the emission coefficients and absorption coefficients of the hybrid electron populations in the source function $S=({j}_{\mathrm{th}}+{j}_{\mathrm{nth}})/({\alpha }_{\mathrm{th}}+{\alpha }_{\mathrm{nth}})$. For comparison, the spectrum includes only thermal synchrotron emission ($S={j}_{\mathrm{th}}/{\alpha }_{\mathrm{th}}$), as shown by dots, for the case of a Keplerian flow. For all model spectra, thermal synchrotron emission dominates the millimeter/sub-millimeter range close to the spectral peak, and non-thermal synchrotron emission is necessary to explain the observed spectrum at lower frequencies.

Zoom In Zoom Out Reset image size

In order to understand the importance of general-relativistic effects, the empty circles shows the profile for the case where the frequency shift

$$\begin{eqnarray}&&\displaystyle \frac{\nu }{{\nu }_{0}}=\displaystyle \frac{{p}_{\alpha }{u}^{\alpha }{| }_{\infty }}{{p}_{\alpha }{u}^{\alpha }{| }_{0}},\end{eqnarray} \tag{ 8 }$$

of the radiation is neglected (i.e., set to 1) in the general-relativistic radiative transfer calculation. In the above equation, ν is the photon frequency, ${p}_{\alpha }$ is its covariant four momentum, ${u}^{\alpha }$ is the four velocity of a fluid particle, and the subscripts "0" and "$\infty$" denote quantities evaluated, respectively, in the local co-moving frame of the flow and of a distant observer. In the lower-frequency region, the surrounding accretion flow remains optically thick and the observed emission is dominated by the photosphere located far from the black hole. Therefore, the frequency shift correction is unimportant in this region. The self-absorbed synchrotron emission in this part of spectrum is dominated by $S={j}_{\mathrm{nth}}/{\alpha }_{\mathrm{th}}$ (Özel et al. 2000). The degeneracy in the dynamics of the three models results in similar spectral profiles.

As the flow becomes optically thin at several hundred GHz (i.e., mm/sub-mm wavelengths), the emission from the flow within a few ${r}_{{\rm{g}}}$ of the black hole becomes observable and the black hole shadow is revealed. In general, the frequency shift reduces the observed luminosity, shifting the thermal peak to lower frequencies. The energy shift, due to Doppler effects arising from the toroidal motion of the flow, results in a larger luminosity in the case of a Keplerian flow compared to other cases with less significant toroidal motion (see Figure 1 for the differences between different flow dynamics). The profiles of Keplerian and free-fall cases may be taken as the limiting cases of a sub-Keplerian RIAF, and therefore represent the boundaries of all possible sub-Keplerian flows (i.e., all possible values of α and β). Interestingly, observational data that show variations at different observation epochs are also located within these boundaries. This indicates that the averaged behavior of non-stationary RIAF flow dynamics at different observation epochs may cause the observed variations. Simultaneous Atacama Large Millimeter/sub-millimeter Array (ALMA) observations at 345, 690, and 890 GHz will be extremely useful in constraining these models in the optically thin window ($\gtrsim 100$ GHz).

The corresponding black hole shadow images at 1.3 mm are dominated by thermal synchrotron emission, as shown in the left column of Figure 3. The position angle of the BH spin is set to be 150° east of north, close to the best-fit parameter of a Keplerian rotating-flow in Broderick et al. (2016). For the case of a Schwarzschild black hole (a = 0), frame-dragging effects due to the rotation of spacetime are absent; consequently, the luminosity contrast between the left- and right-hand sides of the image (divided by the projection of the rotation axis of the flow) are purely determined by the motion of the flow. In Figure 3, the image of the Keplerian model has the largest luminosity contrast (top left panel), and the image of the free-falling flow is axisymmetric (bottom left panel). For the case of a Kerr black hole, the resulting shadow image is asymmetric, even for a flow in free-fall, with the brighter side corresponding to photons that are Doppler-boosted into the observer's line of sight by frame-dragging (see, e.g., Figure 1 of Falcke et al. 2000; compare also with the Ω profiles in Figure 1). Consequently, the resulting crescent structure from the brighter side is the combined effect of both the fluid dynamics and black hole spin in general.

Zoom In Zoom Out Reset image size

In Figure 3, we also show the distributions of the visibility amplitude and phase, which are Fourier-transformed from the model image and observables of the VLBI, along with uv-coverages of future EHT observations with telescopes at eight different sites (see, e.g., Fish et al. 2014; Lu et al. 2016). Despite different dynamics, for all cases, the visibility amplitudes perpendicular to the spin axis are smoother and more extended, compared to the direction parallel to the spin axis. Such characteristics may be used for determining the orientation of the black hole spin, as suggested in Medeiros et al. (2016). In addition, these three models exhibit significantly different visibility distributions, not only at baselines longer than 4 Gλ, as traced in future full-array observations, but also at shorter baselines already observed in previous campaigns (Doeleman et al. 2008; Fish et al. 2011; Johnson et al. 2015) with US stations at California (CARMA), Arizona (SMT), and Hawaii (SMA, JCMT, and CSO). These baselines are indicated, in red, in Figure 3. Therefore, the dynamics of an RIAF are an important factor in interpreting previous EHT observations, as well as future observations.

For a more detailed comparison, in Figure 4, we show the visibility amplitudes of three models on EHT baselines, as well as the data of EHT observations in 2009 (Fish et al. 2011) and 2013 (Johnson et al. 2015). The angular-broadening, due to scattering effects (Johnson & Gwinn 2015), was corrected with the scattering kernel in Bower et al. (2006). Interestingly, year-to-year variations in the visibility amplitude can be well-explained simply through the dynamics, without changing other parameters, such as the black hole spin and inclination angle. One of the most notable differences between the 2009 and 2013 data is the behavior at long US baselines between the US mainland (SMT, CARMA) and Hawaii, with baseline lengths of ∼3 Gλ⁴ ; the amplitude decreased with baseline length in 2009, whereas it increased in 2013.

Zoom In Zoom Out Reset image size

The behavior of the 2009 data can be well-explained with the Keplerian model, as already demonstrated in previous studies (Broderick et al. 2011a, 2011b, 2016). In the case of the Keplerian model, because the black hole shadow is partially covered and smoothed by emission from the approaching side of the toroidal flow, the visibility amplitude has a Gaussian-like behavior in the NEE–SWW direction, where long US baselines are distributed (see top middle panel of Figure 3). On the other hand, rising amplitudes in the 2013 data suggest that the black hole shadow should be more clear than in previous years, which can be well-reproduced with sub-Keplerian or free-fall models. This indicates that year-to-year structural variations in Sgr A* on event-horizon scales potentially reflect the non-stationary dynamics of an RIAF, while simultaneously suggesting that consideration of the flow dynamics is essential in understanding and interpreting data from millimeter/sub-millimeter VLBI observations.

We considered three different accretion flow models: (i) a Keplerian shell rotating outside the ISCO and infalling, with constant angular momentum, inside the ISCO (e.g., Falcke et al. 2000; Broderick & Loeb 2006); (ii) a sub-Keplerian flow (to mimic the flow motion of an RIAF); and (iii) a flow in free-fall, with zero angular momentum at infinity (e.g., Falcke et al. 2000). It was demonstrated that the sub-Keplerian nature of an RIAF is important for modeling black hole images at millimeter/sub-millimeter wavelengths, and hence, crucial for interpreting the observed VLBI visibility and variability of Sgr A*.

Our semi-analytic model of RIAF dynamics are constructed through a combined description of both Keplerian and free-fall fluid motion, and are, therefore, sub-Keplerian everywhere. As first noted by Abramowicz & Zurek (1981), if a hydrodynamical flow has a radial profile that is sub-Keplerian everywhere, plasmas in the vicinity of the black hole are then able to flow into higher latitudes near the horizon, ${r}_{{\rm{h}}}$, resulting in a quasi-spherical geometry. That is, the ratio of height H to cylindrical radius R follows a profile $(H/R){| }_{r\gg {r}_{{\rm{h}}}}\sim 1$ and $(H/R){| }_{{r}_{{\rm{h}}}}\sim 1$, which is consistent with the assumption of the spatial populations described in Equations (4) and (6). It should be noted that we assume our RIAF model characterizes a stationary flow state for given values of α and β (which respectively control u^r and ${\rm{\Omega }}={u}^{\phi }/{u}^{t}$ of the flow), and that we assume these values may vary from observation to observation, due to the sufficiently long dynamical timescale between observations. In addition, magnetic effects are ignored. The simple model is useful for future exploration of the parameter space, providing a quick, flexible overview of the most likely flow dynamic ranges. We plan to perform further detailed analysis, using the closure phase data recently published in Fish et al. (2016). Their data will prove useful to constrain the dynamics of the accretion flow, as well as the black hole spin, its inclination angle, and position angle.

Although GRMHD simulations with an initial equilibrium torus ($(H/R){| }_{r\gg {r}_{{\rm{h}}}}\sim 1$) have successfully simulated sub-Keplerian accretion flows when a stationary state is reached (e.g., McKinney et al. 2012; Narayan et al. 2012), further examination of how $(H/R){| }_{{r}_{{\rm{h}}}}$ is affected by the presence of a magnetic field will be important to understanding the spatial distribution of plasma in the vicinity of the event horizon. For example, it has been discussed in McKinney et al. (2012) that the accumulated polar magnetic flux can compress the flow height when approaching the horizon, resulting in $(H/R){| }_{{r}_{{\rm{h}}}}\ll 1$ (see top panel of Figure 11 in McKinney et al. 2012, and Figure 1 in O'Riordan et al. 2016). Such an understanding is directly related to the shadow image, when the flow becomes optically thin. In addition, compared to that of a simulation in which an initial equilibrium torus is adopted, it is expected that the resulting fluid dynamics of a simulation in which mass is continuously injected onto the grid would be more closer to a free-fall scenario. Such boundary conditions have been applied in many hydrodynamical or MHD simulations with a pseudo-Newtonian potential (Igumenshchev & Abramowicz 2000; Igumenshchev & Narayan 2002; McKinney & Gammie 2002; Igumenshchev et al. 2003). Future comparison between different GRMHD simulation results and our model parameters α and β will help to improve the modeling of RIAF flows.

We thank Ziri Younsi for helpful discussions and proofreading. We also thank Vincent L. Fish and Roman Gold for comments. H.-Y. Pu and K. Asada are grateful to the Greenland Telescope (GLT) team in ASIAA for the support and encouragement. H.-Y. Pu and K. Asada are supported by the Ministry of Science and Technology (MOST) of Taiwan, under the grant MOST 103-2112-M-001-038-MY2. K. Akiyama is financially supported by the program of Postdoctoral Fellowships for Research Abroad at the Japan Society for the Promotion of Science, and grants from the National Science Foundation (NSF). This research has made use of NASA's Astrophysics Data System.