X-RAY DIAGNOSTICS OF GIANT MOLECULAR CLOUDS IN THE GALACTIC CENTER REGION AND PAST ACTIVITY OF Sgr A*

In order to calculate X-ray emissions from illuminated clouds, we constructed a new framework based on Monte Carlo simulation. This framework, which is called MONACO (standing for MONte Carlo simulation for Astrophysics and COsmology), is a general-purpose code applicable to photon-tracking simulations of many astrophysical objects. We utilize the Geant4 simulation toolkit library (Agostinelli et al. 2003; Allison et al. 2006) for tracking particles and building geometry, while physical processes are our original implementations. Geant4, which is widely used for experimental particle physics and nuclear science, provides sophisticated handling of complicated geometries such as radiation detectors.

Figure 1 shows a schematic concept of the Monte Carlo simulation. In the simulation, we build a cloud geometry system containing physical properties such as chemical composition and number density of the cloud; then, in each Monte Carlo trial, a photon with energy E₀ and direction $\bm {\Omega }_0$ is generated at time t₀ and position $\bm {x}_0$. The simulation tracks the particle until it escapes from the system, recording information about the last interaction of the photon $(E_1, \bm {\Omega _1}, t_1, \bm {x}_1)$, which an observer regards as emission. A detailed description of the Monte Carlo simulation and its framework is given in Appendix A.

Zoom In Zoom Out Reset image size

We have developed our original implementation of physical processes adaptable to Monte Carlo simulation frameworks based on Geant4. In neutral matter such as molecular clouds, relevant interactions of X-rays are photoelectric absorption and scattering by bound electrons in hydrogen and helium. Scattering by electrons in heavy elements is negligible because of their small abundances. While the simulation tracks a photon by Monte Carlo methods, an interaction occurs according to its mean free path; then the status of the photon is changed according to differential cross section of the interaction.

As cross section data of photoelectric absorption, we adopt data based on the Evaluated Photon Data Library 97 (EPDL97),⁵ which is provided along with the Geant4 toolkit as the data of Low-Energy Electromagnetic processes version 6.9. A K-shell fluorescent photon is generated with a probability of a fluorescence yield, following a photoelectric effect. As relevant atomic properties, line energies of K-shell fluorescence, K-shell fluorescence yields, and Kβ-to-Kα ratios are taken from Thompson et al. (2001), Krause et al. (1979), and Ertuǧral et al. (2007), respectively. In addition, Kα₂-to-Kα₁ ratios are fixed to 0.5.

Scattering by bound electrons can be divided into three channels according to the final state of the target electron: Rayleigh scattering (elastic), Raman scattering (excitation), and Compton scattering (ionization). The difference in the target between a free electron and a bound electron is significant for the fine structure of the iron line complex including the Compton shoulder. The shape of the Compton shoulder is smeared by an electron binding effect that is known as the Doppler broadening effect due to the finite momentum of the electron bound to an atom or molecule. For atomic hydrogen, cross sections and differential cross sections of the three channels are provided by analytical formulae (Sunyaev & Churazov 1996). For molecular hydrogen, we adopt an approximation that scattering by the molecule is identical to that by atomic hydrogen except that the cross section of Rayleigh scattering is enhanced by a factor of two per electron (Sunyaev et al. 1999). Since the second most abundant element, helium, does not have analytical expressions of atomic properties, we use numerical calculations of differential cross sections (Vainshtein et al. 1998). Appendix B provides details of the implementation of the scattering processes.

For verification of our simulation, we performed a benchmark calculation of an illuminated cloud. Results of this simulation can be directly compared with Figure 9 of Sunyaev & Churazov (1998) under identical simulation conditions. In the simulation setup, a point X-ray source was located at the center of a spherical uniform cloud which has a radial Thomson depth of 0.2 from the center to the surface. Figure 2 shows emerging spectra from the cloud in different time intervals when a short flare occurred at t = 0. The initial spectrum was a power law with a photon index of 1.8 in an energy range of 1–18 keV. Our results nicely agree with the former calculation on the wide-band continuum and the iron Kα line shape including the Compton shoulder for each time interval. In addition, the spectra from our calculation present iron Kβ at 7.1 keV and K-shell fluorescent lines of other metals, which were neglected in Sunyaev & Churazov (1998).

Zoom In Zoom Out Reset image size

We used the simulation framework described in Section 2.1 to model X-ray reflection from the Sgr B2 giant molecular cloud in the GC region. Sgr B2 is one of the brightest XRNe emitting strong 6.4 keV iron fluorescence, which also shows variability over a timescale of a few years. Moreover, since this cloud has been intensively studied in other wavelengths such as radio and infrared bands, its structure and physical properties such as velocity fields and density distribution of molecular gas are relatively well measured compared to other molecular clouds. Sgr B2 is, therefore, the best target for investigating X-ray reflection from molecular clouds.

Figure 3 shows the geometrical setup of the simulation. A photon generator or X-ray source simulating the Sgr A* supermassive black hole is located at the center. The projected distance between Sgr B2 and Sgr A* is fixed to 100 pc. Although a trigonometric parallax of Sgr B2 with the Very Long Base Array (Reid et al. 2009) and X-ray absorption measurements toward the GC region with the Suzaku satellite (Ryu et al. 2009) suggested that Sgr B2 is positioned in front of the GC, there still remains uncertainty of the position of the cloud along the line of sight. The line-of-sight position is important in studying the structure of the Galaxy and also makes a difference in the time of the Sgr A* outburst. We assumed three different line-of-sight positions of the cloud: y = +100 pc, 0 pc, and −100 pc, to investigate effects of the relative position of the cloud to the center. In the present paper, we assume that the distance to the GC from the Earth is 8 kpc.

Zoom In Zoom Out Reset image size

The Sgr B2 complex is the most massive giant molecular cloud in our Galaxy and has a complicated structure containing several active star-forming H ii regions. As intensive observations in radio bands have provided models of the cloud's structure (Lis & Goldsmith 1990; Hasegawa et al. 1994; de Vicente et al. 1997), it can be presumed that the cloud consists of three components: (1) very dense cores with diameters of ∼0.5 pc containing star-forming regions, (2) a dense envelope that seems to correspond with an XRN with a number density of hydrogen molecule of $n_\mathrm{H_2}=10^{4}\hbox{--}10^{5} \;\mathrm{cm^{-3}}$, and (3) a large, diffuse component with an almost constant density of $n_\mathrm{H_2}\sim 10^{3}\;\mathrm{cm^{-3}}$ extending to a diameter of 45 pc.

For simplicity, we constructed a spherically symmetric geometry model composed of a single dense core and an envelope that has a diameter of 10 pc and a realistic density distribution. We determined the diameter of the envelope according to Hasegawa et al. (1994), and the value well agrees with X-ray morphology obtained by Chandra observations (Murakami et al. 2001). In the model, the radial profile of the number density in the cloud is represented as

$$\begin{equation} n(r)=\left\lbrace \begin{array}{@{}lll}5\times 10^6\;{\rm cm}^{-3} &\!({\rm core{:} }\, r<0.25\;{\rm pc}) \\ n_0 \left(\frac{r}{r_0}\right)^{-\alpha } \;{\rm cm}^{-3} &\!({\rm envelope{:}}\,0.25\;{\rm pc}<r<5.0\;{\rm pc}). \end{array} \right. \end{equation} \tag{ 1 }$$

The mass of the core is 2 × 10⁴ M_☉, which is a very small fraction of the total mass of the cloud. Table 1 shows parameters of our cloud models in the simulation. We prepared four models with different n₀ values corresponding to total cloud masses from 2.5 × 10⁵ M_☉ to 2 × 10⁶ M_☉ (Models 1–4). In Models 1–4, the power-law index α = 1 was assumed, as radio observations have reported values of 2 (Lis & Goldsmith 1990) or 0.87 (de Vicente et al. 1997). To investigate effects of different density profiles, we built a model with α = 0 (Model 5) and a model with α = 2 (Model 6), fixing the total mass to 5 × 10⁵ M_☉. While we assumed a metal abundance of 1.5 protosolar value in Models 1–6 as a standard value in the GC region (Nobukawa et al. 2010), we additionally checked two models of different values of 1.0 protosolar (Model 7) and 2.0 protosolar (Model 8). The values of the protosolar (solar system) abundances are taken from Lodders (2003). We ignored the third component surrounding the dense envelope to reduce computation costs; instead, we introduced absorption of the initial spectrum by the third component as

$$\begin{equation} F(E)\propto \exp (-N_\mathrm{H}\sigma _{\rm abs}(E))E^{-\gamma }\quad (1\;\mathrm{keV}<E< 400\;\mathrm{keV}), \end{equation} \tag{ 2 }$$

where N_H = 6 × 10²² cm⁻² is an equivalent hydrogen column density of the surrounding diffuse component and σ_abs(E) is the photoelectric absorption cross section at energy E per hydrogen. The photon index γ of the initial spectrum is 1.8, which is consistent with flares of Sgr A* or X-ray emission from Seyfert galaxies.

We investigated morphologies of iron 6.4 keV line emissions and hard-X-ray emissions by extracting images from the simulation results. Since X-ray emissions extend over the whole cloud, a duration of the illuminating flare of Sgr A* must be longer than the light crossing time of the cloud. As the exact value of the duration does not significantly affect the emissions unless it is shorter than the light crossing time, we assumed that the flare began at t = −65 yr, kept a constant luminosity, and then ended at t = 0. The duration of 65 yr is twice as long as the light crossing time of the 10 pc large cloud. The X-ray reflection from the cloud becomes brightest just before the rear end of the direct light front reaches the cloud edge confronting the illuminating source. Time of the brightest moment depends on the relative position of the cloud to the flare source and the observer, as summarized in Table 2. The assumed luminosity of the flare is L_X = 1.3 × 10³⁹ (d/100 pc)² erg s⁻¹ in an energy range of 1–10 keV, where $d=\sqrt{x^2+y^2}$ is the distance between the cloud and the illuminating source. This is normalized so that the iron line flux of Model 2 (y = 0 pc) at the brightest moment agrees with an observational value of 1.0 × 10⁻⁴ photons s⁻¹ cm⁻² obtained with Chandra in 2000 (Inui et al. 2009).

Figure 4 depicts iron 6.4 keV line morphologies at the brightest moment for the eight different cloud models and three different line-of-sight positions. The X-ray illuminates the cloud from the right side of each image. The case of the middle position (y = 0 pc) is shown in the middle row of Figure 4. In this case, the cloud becomes brightest at t = 303 yr. As the cloud mass increases (from Model 1 to 4), the cloud becomes opaque, making the morphology semicircular. Corresponding hard-X-ray images extracted in an energy range between 20 keV and 60 keV are shown in the middle row of Figure 5. The hard-X-ray morphologies are obviously different from the iron fluorescence due to their large penetrating power. Hard X-rays are therefore useful for probing the structure of the dense molecular clouds. Images of Models 2, 5, and 6 show the difference due to the density profile of the cloud. Since the cloud with a steeper radial profile (Model 6) has denser inner layers and thinner outer layers, its morphology becomes more concentrated and concave. In Model 5, the small opacity of the uniform envelope exposes the central core in the hard-X-ray band. The metal abundance in the cloud also changes the iron line emission while the hard-X-ray images are similar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If the line-of-sight position of the cloud is different, the incident direction of the illuminating photons changes and significantly different morphologies appear. The top rows in Figures 4 and 5 show X-ray morphologies at the brightest moment, t = 756 yr, in the front-illuminated case or the position of y = +100 pc. The delay of the brightest moment is due to light propagation via the cloud behind the GC. In this case, the iron line emissions form into a round shape since a large proportion of the cloud front is illuminated by the flare. On the other hand, the bottom rows of Figures 4 and 5 show morphologies in the back-illuminated case or the position of y = −100 pc, and the reflection becomes the brightest at t = 122 yr. Morphologies of the iron line form into a concave shape; moreover, the shape becomes a crescent if the cloud is more massive than ∼1 × 10⁶ M_☉. In contrast to the iron emission, morphologies of scattered hard X-rays do not largely depend on the relative position of the cloud or the direction of the illuminating source unless the cloud is Compton thick.

Since Sgr B2 displays fading of the X-ray emission (Koyama et al. 2008; Inui et al. 2009; Terrier et al. 2010), we would observe evolution of the X-ray morphologies in the XRN model as shown in Figure 6. At 6.4 keV the central region of the cloud will become darker in the fading phase than outer regions if the cloud is massive enough to absorb most of the incident soft- X-ray photons in the central regions. In contrast, hard X-rays do not suffer from such significant absorption; thus, morphological disagreement between the iron line and the hard X-ray is notable in the decay phase. An apparent light front intersecting the cloud from the observer moves fast when the cloud is positioned in front of the flare source (Sunyaev & Churazov 1998). For the cloud at y = −100 pc morphological evolution is very fast even in the hard-X-ray band (Figure 6, bottom).

Zoom In Zoom Out Reset image size

In addition to an image of the reflected emission, its spectrum contains important information about the cloud and the flare. Figure 7 shows spectra extracted from the simulation results for the models with different masses (Models 1–4). While the assumed incident spectrum is a power law with a photon index of 1.8 absorbed by the outer diffuse region (N_H = 6 × 10²² cm⁻²), absorption due to the interstellar medium is not considered here. Since scattered hard X-rays above 20 keV do not significantly suffer from absorption, their flux increases with the cloud mass. In contrast to the straightforward behavior of hard X-rays, dependence of the iron line on the cloud mass is complicated due to heavy absorption in dense clouds. In addition to the iron lines, fluorescence lines of other heavy elements (e.g., Cr, Mn, Co, Ni) can be observed with high energy resolution. Spectral evolutions are considerably different in each cloud mass since absorption plays an important role in the fading phase of the reflected emission. A more massive cloud displays faster decay, enhancing the line equivalent width and the Compton shoulder.

Zoom In Zoom Out Reset image size

In order to evaluate the spectra quantitatively, we extracted the hard-X-ray flux (20–60 keV), the iron line flux, the equivalent width of the iron line, and the shoulder-to-peak ratio of the Compton shoulder as a function of time. The iron line was divided into a peak and a shoulder. The peak flux and shoulder flux are integrated over energy ranges between 6.0 keV and 6.38 keV, and in a range between 6.38 keV and 6.42 keV, respectively. A continuum is evaluated by fitting to a power law over a range between 5.5 keV and 6.8 keV without several fluorescent line energies, and is then subtracted from both fluxes. The equivalent width and the shoulder-to-peak ratio were calculated using these values.

Figure 8 shows the spectral parameters as a function of time elapsed since the end of the Sgr A* flare. After t = 303 yr, the hard-X-ray and iron line fluxes gradually decrease while the equivalent width and shoulder-to-peak ratio increase with time. Although the iron line flux is almost independent of the cloud mass in the brightest phase, its decay speed depends on the mass because of absorption. The equivalent width of the iron line keeps constant and does not depend on the cloud mass when the whole cloud is illuminated. In the fading phase, however, its variation significantly depends on the mass; it shows a rapid increase in the dense clouds. One reason for this is a difference in photoelectric absorption probability (Sunyaev & Churazov 1998). In a dense cloud, continuum photons around 6 keV are largely and rapidly absorbed, though hard X-rays above 10 keV, which are responsible for producing the iron line, are not so much. As a result, a large equivalent width of the iron line appears. Since the Compton shoulder is generated via scattering of the iron fluorescence, it is also a good probe of the Thomson depth of the cloud or the total mass of the cloud. The equivalent width and the shoulder-to-peak ratio also depend on the metal abundance in the cloud, as shown in Figure 9. The equivalent width is nearly proportional to the metal abundance. A large metal abundance increases photoabsorption probability, suppressing scatterings to generate the shoulder. Note that the metal abundance controls amounts of all metals, not only iron; changing only the iron abundance could bring different results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For different line-of-sight positions of the cloud (Figures 10 and 11), the dependences of the iron line flux on the cloud mass are different from the cloud at y = 0 pc, while other parameters behave in the same way. Since the flux dependence is a monotonic function of the cloud mass, these behaviors of the iron line are somewhat easier to interpret than those of the cloud at y = 0 pc. For the front-illuminated cloud at y = +100 pc, it is an increasing function because it depends simply on the total mass of matter that generates fluorescence. However, the line flux is saturated at M_cloud ∼ 1 × 10⁶ M_☉. For the back-illuminated cloud at y = −100 pc, on the other hand, the dependence is a decreasing function due to heavy absorption in a dense cloud. Moreover, the line-of-sight position makes a difference to the decay time of the hard-X-ray emission without any significant biases of the cloud structure. In contrast, the time evolution of the iron line flux depends not only on the position but also strongly on the amount of absorbing matter. The decay time of the hard X-rays to decrease by half from the brightest moment is 13, 23, or 28 yr for the different positions along the line of sight of y = −100, 0, or +100 pc, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If we assumed that radiation did not affect the surrounding medium, X-ray reprocessing (fluorescence or scattering) in the system would be regarded as a linear response of the system. The linear nature of the system greatly simplifies the problem because the final results are obtained by convolving the response function, which describes how the system changes the initial photon into the final state of the photon. In other words, the problem is reduced to obtaining a Green's function of a linear operator that describes the system. X-ray reflection by cold matter can be treated in this manner.

Thus, the purpose of the Monte Carlo simulation is to calculate a conditional probability density of emitting a photon with an energy E₁ and direction $\bm {\Omega }_1$ at time t₁ and position $\bm {x}_1$ under initial conditions of a photon emitting with an energy E₀ and direction $\bm {\Omega }_0$ at time t₀ and position $\bm {x}_0$. Figure 1 shows a schematic concept of the simulation. In each Monte Carlo trial, we track a photon from certain initial conditions until the photon escapes from the system, recording the initial condition and the conditions of photons at the last interaction before escaping. After M trials of Monte Carlo simulation, the conditional probability density can be written as

$$\begin{eqnarray} &&p(E_1, \bm {\Omega }_1, t_1, \bm {x}_1|E_0, \bm {\Omega }_0, t_0, \bm {x}_0)\nonumber\\ &&= \frac{g(E_1, \bm {\Omega }_1, t_1, \bm {x}_1|E_0, \bm {\Omega }_0, t_0, \bm {x}_0)}{\beta (E_0, \bm {\Omega }_0, t_0, \bm {x}_0)}, \end{eqnarray} \tag{ A1 }$$

where g denotes the number of events that have an initial condition of $X_0=(E_0, \bm {\Omega }_0, t_0, \bm {x}_0)$ and a final condition of $X_1=(E_1, \bm {\Omega }_1, t_1, \bm {x}_1)$, and β denotes a distribution function of initial conditions of the simulation, which is normalized as

$$\begin{equation} \int \beta (X_0) dX_0 = M. \end{equation} \tag{ A2 }$$

Then, we can write a distribution function of emitting photons in the system as

$$\begin{equation} f(X_1)=\int p(X_1|X_0) \alpha (X_0) dX_0=\int g(X_1|X_0)\frac{\alpha (X_0)}{\beta (X_0)}dX_0, \end{equation} \tag{ A3 }$$

where α(X₀) is a distribution function of the initial photons that are emitted by the source, or simply the source function. The simulation output is a list of the escaping photons with information on the initial and final conditions. For the list, f(X₁) is calculated by

$$\begin{equation} f(X_1)=\sum _{i=1}^{M} w_i(X_1^{\prime }, X_0), \end{equation} \tag{ A4 }$$

$$\begin{equation} w_i(X_1^{\prime }, X_0) = \left\lbrace \begin{array}{@{}ll}\frac{\alpha (X_0)}{\beta (X_0)} & (X_1^{\prime } = X_1) \\ \ms 0 & (X_1^{\prime }\ne X_1). \end{array} \right. \end{equation} \tag{ A5 }$$

Here, w_i is a weight of each event and the summation is performed over all Monte Carlo events. When the geometry of the system is static or has a symmetry of time translation, it is possible to fix the initial time t₀ in the simulation to zero; then Equation (A3) is reduced to

$$\begin{eqnarray} f(E_1, \bm {\Omega }_1, t_1, \bm {x}_1)&=&\int g(E_1, \bm {\Omega }_1, t_1-t_0, \bm {x}_1|E_0, \bm {\Omega }_0, 0, \bm {x}_0)\nonumber\\ &&\times\frac{\alpha (E_0, \bm {\Omega }_0, t_0, \bm {x}_0)}{\beta (E_0, \bm {\Omega }_0, 0, \bm {x}_0)}dE_0 d\bm {\Omega }_0 dt_0 d\bm {x}_0.\nonumber\\ \end{eqnarray} \tag{ A6 }$$

To simulate observations, we must consider time differences of light propagation to an observer from different positions. A distant observer directed toward a unit vector $\bm {r}$ observes the emission at $t_{\rm obs}= t_1-\bm {x}_1\cdot \bm {r}/c$, where c is the speed of light and the origin of t_obs is determined so that an emission at t₁ = 0 and $\bm {x}_1=\bm {0}$ is observed at t_obs = 0 through direct propagation (without any interactions). Events are selected to obtain observational spectra or images if the direction agrees with $\bm {r}$ within a certain small tolerance. The weight of the selected event is calculated by Equation (A5) in which $t_1=t_{\rm obs}+\bm {x}_1\cdot \bm {r}/c$ is substituted. If the system has spatial symmetry such as spherical or axial symmetry, we can extend the range of the direction selection to reduce statistical uncertainties unless the symmetry breaks.

Our simulation framework, MONACO, provides methods for the Monte Carlo simulation and for the convolution to generate observed spectra or images (Appendix A.1). An outline of the calculation to obtain observational data is shown in Figure 13. The first part uses Monte Carlo simulations to calculate the probability of obtaining certain final conditions, given certain initial conditions. This part requires information about the initial conditions of the photon and the geometry including physical conditions of matter as input. The output here is a list entry for an event that has information about the initial state and the final state (emission). If a photon does not experience an interaction with matter, the final state is identical to the initial state.

Zoom In Zoom Out Reset image size

Using the event list generated by the simulation part as an input, the second part, referred to "convolution" in Figure 13, generates simulated observational data such as spectra or images by performing the convolution with the source function described in Appendix A.1. The output here is the final results, which can be compared with real observations. The list input can be a combination of different event lists to reduce statistical uncertainties or to extend the ranges of initial conditions. To perform the convolution by Equation (A5), this part requires as inputs the initial conditions of photons both for the simulation and for the real source (source function). It also requires information on the observer: direction, time, distance, and miscellaneous parameters.

The first part (simulation part) performs the photon-tracking simulation including treatment of X-ray reprocessing. This part is an extended and improved version of original Monte Carlo codes developed by Watanabe et al. (2003, 2006). We use the Geant4 simulation toolkit library (Agostinelli et al. 2003; Allison et al. 2006) only for treatment of complicated geometry, while physical processes are our original implementations. The simulation framework consists of a manager and four user-defined modules. The manager, which is provided by the Geant4 library, handles communications between the modules and treats particle tracking based on Monte Carlo simulations and geometry navigation. The user-defined modules are for geometry, physical processes, primary generation, and data acquisition. The geometry module describes the geometry of the astrophysical object to be simulated as well as physical conditions of the matter. The module for physical processes organizes X-ray interactions with matter. The primary generator describes how to create initial photons to be tracked. For example, a point-like source is specified by a position vector of the source and an energy spectrum such as a power law with a photon index and an energy range; in each Monte Carlo trial, the direction and energy of the primary photon are sampled from an isotropic distribution and the power-law distribution, respectively. The data acquisition module forms an event list of the simulation data that contains information on the initial and final states of the photon and saves it to a file.