Validation of Radiative Transfer Computation with Monte Carlo Method for Ultra-Relativistic Background Flow
We developed a three-dimensional radiative transfer code for an ultra-relativistic background flow-field by using the Monte Carlo (MC) method in the context of gamma-ray burst (GRB) emission. For obtaining reliable simulation results in the coupled computation of MC radiation transport with relativistic hydrodynamics which can reproduce GRB emission, we validated radiative transfer computation in the ultra-relativistic regime and assessed the appropriate simulation conditions. The radiative transfer code was validated through two test calculations: (1) computing in different inertial frames and (2) computing in flow-fields with discontinuous and smeared shock fronts. The simulation results of the angular distribution and spectrum were compared among three different inertial frames and in good agreement with each other. If the time duration for updating the flow-field was sufficiently small to resolve a mean free path of a photon into ten steps, the results were thoroughly converged. The spectrum computed in the flow-field with a discontinuous shock front obeyed a power-law in frequency whose index was positive in the range from 1 to 10 MeV. The number of photons in the high-energy side decreased with the smeared shock front because the photons were less scattered immediately behind the shock wave due to the small electron number density. The large optical depth near the shock front was needed for obtaining high-energy photons through bulk Compton scattering. Even one-dimensional structure of the shock wave could affect the results of radiation transport computation. Although we examined the effect of the shock structure on the emitted spectrum with a large number of cells, it is hard to employ so many computational cells per dimension in multi-dimensional simulations. Therefore, a further investigation with a smaller number of cells is required for obtaining realistic high-energy photons with multi-dimensional computations.
keywords:Gamma-ray burst, Relativistic jet, Radiative transfer, Monte Carlo method
Relativistic radiation hydrodynamics computation is widely used in the field of high-energy astrophysics, e.g., gamma-ray bursts (GRBs), supernovae, accretion discs, and active galactic nuclei. It is difficult to analytically solve the radiative transfer equation since it is a seven-dimensional equation; therefore, some approximated methods have been proposed to obtain physically reasonable solutions (1); (2); (3).
The moment method, which wraps up angular dimensions of photons, is frequently used in radiation hydrodynamics computation (1). However, it is valid only for a closely isotropic radiation field (that is, an optically thick regime), unless an appropriate closure is given. For a highly relativistic regime, it is inadequate because the radiation field holds strong anisotropy as a result of the beaming effect. The moment method was also formalized in the relativistic regime, but it is applicable for simple limiting cases (4). The discrete-ordinate method, which solves the radiative transfer equation with the finite differencing of direction components, is suitable for an anisotropic regime (5). However, it generally requires huge computational costs for resolving six-dimensional phase space. The Monte Carlo (MC) method is useful to statistically obtain a solution of a multi-dimensional integro-differential equation and then to solve the radiative transfer equation including scattering process. It is appropriate in an optically thin and a mildly scattering regime such as a radiation field accompanying a relativistically expanding jet, resulting in GRB emission.
GRBs are highly energetic explosion phenomena in which extremely large amounts of energy are emitted in a few seconds to minutes. In particular, long GRBs of duration greater than 2 s are thought to occur in association with relativistic jets formed around collapsing massive stars and have a collimated configuration, while short GRBs of duration less than 2 s are considered to result from a binary merger of compact objects, e.g., neutron star–neutron star or neutron star–black hole. The prompt gamma-ray emission produced in internal shocks of an ultra-relativistic jet is typically interpreted as a result of synchrotron radiation from shocked-accelerated electrons (6); (7). However, the internal shock model could provide insufficient radiation efficiency for explaining GRB emission (8); (9); (10). On the other hand, the effect of a photosphere position at which thermal photons emerge has been discussed for the thermal component of GRB spectra (11); (12). Such a photospheric emission model, in which thermal photons are Comptonized at the shock front, has high radiation efficiency for GRBs. The structure of relativistic jets, which develops with time, has also been studied through multi-dimensional relativistic hydrodynamics simulations in the context of GRBs (13); (14); (15); (16); (17); (18), and light curves and spectra have been estimated based on such simulations (17); (19); (20); (21); (22). Although the observed spectra can be characterized by a broken power-law shape (23), their spectral properties have not been accurately reproduced in the numerical works.
The MC radiative transfer computations have been implemented in the relativistic regime relevant to GRBs (31); (32); (33); (34); (35). Numerical studies with the MC technique have been also reported for explaining radiation and neutrino transport in supernova explosions (24); (25); (26); (27); (28); (29); (30). Some observations indicate the GRB spectra includes a thermal component (36); (37); therefore, the radiation transport of thermal photons produced at the photosphere has been paid attention, and the non-thermal feature of the spectra was obtained by overlapping thermal spectra with various escaped angles at different time (38). The non-thermal spectra were also obtained by taking into account the gradual energy dissipation by magnetic reconnection (39); (40). On the other hand, the structure of an ultra-relativistic shock wave in the self-similarly expanding fluid has been numerically investigated (41), and the high-energy component of GRB spectra could be explained by bulk Compton scattering in such a shock wave by using MC computations (42); (43). The bulk Compton scattering occurs when photons traveling across the shock wave collide with relativistic electrons. The non-thermal spectra of GRBs have been also explained through bulk Comptonization in an ultra-relativistic jet with some shells of various flow velocities (44); (45).
Past works with multi-dimensional relativistic hydrodynamics simulations showed that structure of a relativistic jet exhibits highly inhomogeneous developing, which can affect the observed spectra (13); (14); (16); (17). Although some radiative transfer simulations were conducted with an ultra-relativistic steady flow-field, those with a time-dependent flow-field may have a significant impact on the detailed analysis of GRB spectra. Relativistic radiation hydrodynamics simulations were implemented with the MC radiative transfer; however, some of them do not appropriately take into account the feedback from interaction of radiation with flow-field matter (46); (47), while others do not sufficiently perform the test calculations for the ultra-relativistic flow-field (48). On the other hand, radiation hydrodynamics calculations for predicting the emission of internal magnetized shocks including the feedback of the radiation on the dynamics were examined in the relativistic regime without the MC technique (49). The coupling of the MC radiative transfer with relativistic hydrodynamics has not been sufficiently performed because of computational difficulties as introduced in Sec. 2. In this paper, we show the validation of the formalism to include the MC method in the ultra-relativistic flow-field and examine the appropriate simulation conditions such as time interval for developing the flow-field and spatial resolution for obtaining converged results.
The remainder of this paper is organized as follows. We introduce the difficulty in performing MC radiative transfer computation in the ultra-relativistic flow-field in Sec. 2, and present the numerical method in Sec. 3. The validation of the radiative transfer computation in different inertial frames is presented in Sec. 4, and the effect of the flow-field resolution on radiation transport is discussed in Sec. 5. We summarize this paper in Sec. 6.
2 Difficulties in performing Monte Carlo radiative transfer computation coupled with ultra-relativistic hydrodynamics
Some MC methods for radiation transport have been developed for the coupling with non-relativistic hydrodynamics computation in previous works as follows. Implicit MC (IMC) schemes have been employed for efficient computations in optically thick systems, and they allow employing larger time steps, , in the numerical simulation of the flow-field (50); (51). The diffusion approximation which can reduce a computational cost is reasonable in opaque regions, so a hybrid method consisting of IMC and MC diffusion method has been developed (52); (53); (54). Since the MC method includes statistical errors with a small number of samples, some techniques for reduction of statistical noise have been investigated with a moderate computational cost (55); (56). These techniques are actually effective in the non-relativistic regime.
In the ultra-relativistic regime (Lorentz factor ), however, since the velocity of matter is almost the same as the speed of light, computation with an excessively large for developing the flow-field leads to false judgment on whether a photon crosses the shock front (Fig. 1). Photons in the upstream side of the shock wave are scattered with a different probability from those in the downstream side. If the Compton scattering is considered, scattered photons not only shift their directional angles but also undergo energy exchange with matter. Therefore, the simulation of radiation transport could not produce accurate results in directional-angle distribution and an energy spectrum with a mistake in judgment on whether a photon crosses the shock front due to the large . A small is essentially necessary in the ultra-relativistic regime, so the IMC scheme for a large might be meaningless, and it is controvertible to develop the IMC scheme for highly relativistic situations. Since a strongly anisotropic flow should be taken into account in the ultra-relativistic regime because of the beaming effect, the diffusion approximation cannot be applicable.
In computations including feedback from interaction between photons and flow-field matter, a huge computational cost is needed. Since the flow-field profile is affected by the radiation transport, the post-process computation of the radiation transport in the steady background flow-field computed in advance is not adequate. Therefore, the computations of radiative transfer and hydrodynamics should be alternately performed. Furthermore, a number of photons should be put in each cell of the hydrodynamics computation for obtaining converged simulation results, however, it leads to a huge-computational-cost problem. Thus, some techniques for reduction of statistical noise with a small number of samples as stated above are required, and the coupled computation may be feasible if restricted to some regions in the flow-field.
Moreover, numerical diffusion is inevitable in hydrodynamics simulation, which especially affects shock structure. Photons near the shock front are scattered with the bulk-Compton-scattering process and obtain energy from relativistic electrons in the flow-field according to the local flow velocity and density. The local flow velocity affects the obtained energy from the scattering, and density is related to the probability of the scattering. The shock structure, therefore, affects the spectrum resulting from the radiative transfer computation. Unphysical finite width of a smeared shock front due to numerical diffusion depends on spatial resolution of the hydrodynamics computation. Appropriate simulation conditions should be assessed such as the time interval for updating the background flow-field and the spatial resolution for feasible computations of radiation transport before the coupled computation of the MC radiative transfer with ultra-relativistic hydrodynamics can be performed.
In the present study, we construct radiative transfer algorithm for an ultra-relativistic background flow-field with consistent transformation between a comoving frame (CMF) and an observer frame (OBF), and assess the interval of time-step values which yield accurate computations. The effect of the spatial resolution of the background flow-field on the radiative transfer computation is also discussed.
3 Numerical method
We developed a simulation code of radiation transport for a highly relativistic background with the MC technique. The radiative transfer equation in a certain inertial frame that takes scattering into account is expressed as follows (1):
where is the specific intensity, which is a function of the position vector , the traveling direction vector , photon frequency , and time . The symbols and denote the speed of light and density, respectively. The absorption and scattering cross-sections are denoted by and , respectively. The emissivity depends on the matter temperature . In the integrand, , , , and denote the scattered direction, incident direction, scattered frequency, and incident frequency, respectively. The fluid velocity varies depending on the position; therefore, the cross-sections are dependent on in an arbitrary inertial frame.
Photons are thermally emitted in the flow-field and travel in a medium in which scattering opacity is assumed to be greater than absorption opacity. Thomson and Compton scattering processes are taken into account. We adopt the MC technique to solve the radiative transfer equation as in the previous works (33); (35).
By tracking a large number of particles, we can obtain an approximate solution of the equation with reasonable computational costs. A cluster of photons (called as ‘packet’) having single frequency is considered as a sample particle. The packet has energy of , where and are Planck constant and the number of photons in the packet, respectively. Emission power at any position is described as the spectral integration of the emission coefficient in the CMF, , given by
where is the photon frequency in the CMF; the initial value of is set randomly as weighted for Planck energy distribution in the CMF. The holding energy for single packet in the CMF, , can be obtained as product of , the total cell volume, and a certain time duration divided by the total number of emitted packets during the time interval, . Since is the energy in the CMF, we need to Lorentz-transform it to that in the OBF. The frequency in the CMF, , is transformed to that in the OBF as well. The initial direction of the photon is also set by two random numbers with the assumption of isotropic emission in the CMF; the initial direction vector can be expressed by
where and are random numbers for the direction. The traveling direction in the CMF, , is transformed to obtain the one in the OBF, , as follows:
where is the Lorentz factor of the flow velocity. Using the direction and the flow velocity , the frequency and packet energy in the OBF are given by
Photons are categorized into two energy groups because they have different cross-sections depending on their energy regime. The lower-energy photons are called as ‘optical ray’, and the higher-energy ones are expressed by ‘gamma ray’. The threshold of two groups is set to 10 keV. The scattering cross-section for optical ray is Thomson scattering one, . The electron number density, , is estimated by assuming that the species involved in the background is only fully-ionized helium gas in this paper, so the mass scattering cross-section in the CMF, , for optical ray can be calculated as . In the Thomson scattering, the energy of the incident photon in the CMF is fully transferred to the scattered one. The scattered angle and direction are determined by assuming unpolarized light and using a phase function , where is the scattering angle.
Gamma-ray photons are scattered in the Compton regime. Thus, the corresponding Compton cross-section is employed for them. The energy of the incident gamma ray is transferred to the scattered light and matter, and the photon loses its own energy to become an optical ray through multiple scatterings. The scattering cross-section in the CMF can be estimated by integrating the Klein-Nishina (K–N) formula. The scattering angle, , can be determined by a random number, , employing the K–N formula;
where and are the energy of the incident photon in the CMF and K–N scattering cross-section, respectively. Here, the differential K–N cross-section is expressed as follows:
where and . Here, and are the electron rest mass and the elementary charge, respectively. The total K–N cross-section is obtained by integrating Eq. (9) over the solid angle (57). The angle in a plane perpendicular to is randomly chosen by , where is a random number. The energy of the incident packet is divided into the scattered light and the scattered electron, and the fraction of the energy of the scattered light is denoted as . So, the energy of the scattered light in the CMF, , is . Here, electron energy spectrum is not taken into account for simplicity.
A free path of a photon in the CMF, , between the collisional events depends on the cross-section in the CMF. The optical depth of the matter corresponding to is given by
The free path of the photon is determined by probabilistic manner using . The probability that the photon freely travels can be represented by
The uniform random number for the free path, , equals to the integrated probability in the following form:
so that the optical depth is readily expressed with it;
Therefore, the flying distance of the photon, , is given by
Since photon transport is treated in the OBF, the free path in the CMF, , should be transformed to that in the OBF, , by
The time it takes for the packet to travel a distance at light speed is . The time for the packet should be incremented to between the events.
We validate our simulation code of radiative transfer through two test problems: (1) comparing the simulation results calculated in different inertial frames, and (2) evaluating the effect of the smeared shock front mimicking numerical diffusion on the radiative transfer simulation as a preliminary result for the coupled computation of radiative transfer with relativistic hydrodynamics.
4 Comparison among different inertial frames
We performed radiative transfer simulations in a shock rest frame and two shock moving frames. In the MC radiation transport computation, the part of collision with electrons is treated in the CMF, and that of transport of each photon is treated in the OBF. Thus, the Lorentz transform between the CMF and OBF is performed. By comparing the results computed in the different frames, we verify that our simulation code can properly handle Lorentz transformation and determine the appropriate computational conditions for obtaining converged solutions.
4.1 Simulation condition
A relativistic flow-field with a steady shock wave was set up, and a cylindrical coordinate system with one cell in the -direction and two cells in the -direction was adopted. The density , pressure , and flow velocity are the physical parameters at the upstream side of the shock, as shown in Fig. 2. Similarly, the physical parameters , , and are the values at the downstream side. These quantities satisfy the Rankine–Hugoniot (R-H) relations between the upstream and downstream sides of the shock wave; therefore, the cell interface in the -direction represents the shock front. The width of the cell in the -direction, , is determined by the optical depth measured from the boundary at the upstream side along the -axis as
such that and at the upstream and downstream sides, respectively, as shown Fig. 2. Here, the density and the scattering cross-section are the values in the CMF, and is the Lorentz factor of the flow velocity in the shock rest frame. The widths of the computational cells in the shock moving frames can be fixed through the Lorentz contraction. To make the computational cells equivalent among different inertial frames, the length is measured by means of . The width of the computational cell in the -direction is set to be large enough not to affect simulation results. All photons are emitted once at a pre-established initial time and at one point located at an optical depth in the -direction. The location and timing of emitted photons are coincident on a space-time diagram among each frame. A lot of photons are tracked until the photons reach the boundary of the computational domain. We record directional angles and energy of the photons at the boundary and sum them for all the photons to obtain directional angular distribution and a spectrum.
The relativistic shock wave in the background can be described by solving the following relativistic R-H relations:
where , , and are the four-velocity, specific enthalpy, and pressure, respectively. Here, , where is the Lorentz factor and is the three-velocity of the th-direction component. The specific enthalpy is defined by , where is the specific internal energy given by the equation of state for an ideal gas, . A constant specific heat ratio is assumed. The upstream quantities , , and and the shock speed are given; then, the downstream quantities , , and can be obtained by solving these relations. The physical quantities in the shock rest frame are adopted as summarized in Table 1. These quantities are employed by reference to the previous work on the relativistic hydrodynamics computation of a jet (17), in which the flow-field transits from the optically thick to the optically thin regime. Here, the value of the flow velocity corresponds to , where is the Lorentz factor of the upstream flow velocity in the shock rest frame, and is equal to the Lorentz factor of the shock speed , , in the rest frame of the upstream flow velocity.
The temperature of the matter is calculated from the equation of blackbody radiation by assuming the radiation equilibrium gas
where is the Stefan–Boltzmann constant. The emission coefficient and initial frequency of the packet are determined according to , and electrons are assumed to be at rest in the CMF. Here, in the downstream side of the shock wave corresponds to . The radiation pressure is assumed to be equal to the total pressure for simplicity.
4.2 Transformation of inertial frames
We transform the flow velocities from the shock rest frame to the shock moving frames in both the upstream and downstream sides of the shock wave to prepare different frames for the same shock-wave strength. The transformation is implemented according to the following equation (58):
where is the flow velocity in the shock rest frame, is the flow velocity in the shock moving frame, and is the relative velocity between two frames (i.e., the apparent shock speed in the arbitrary frame). The flow velocities in each observer frame are shown in Fig. 3.
We compare the simulation results for three different inertial frames in this paper. One is the shock rest frame, and the other two are shock moving frames corresponding to and , where is the Lorentz factor of the apparent shock speed, . Here, the situation with corresponds to the rest frame of the upstream flow velocity.
To ensure that the computational grids on the space-time diagram are consistent between the shock rest frame and the shock moving frames, the computational grids move at the shock speed in the shock moving frames, and the shock front is located at the cell interface in each time step so that it has an exact discontinuity at all times. The shock front is assumed to be static during a time step . The radiation transport is computed on the flow-field with the static shock front during the current . Then, the position of the shock front is updated by the next , and the radiation transport is computed on the updated flow-field with the static shock front again. The shock front movement and the radiation transport are repeated alternately.
4.3 Constraint of time interval
We consider the required constraint of the computational time interval updating the flow-field, , for obtaining converged solutions from the angular distribution of the photons escaping from the computational domain. The angular distribution between the photon traveling direction and the -axis with various values, which are calculated in the frame with sample particles and transformed to the shock rest frame, are shown in Fig. 4 (a). For simplicity, we ignored Compton scattering and only took into account Thomson scattering with the low-frequency approximation of the K-N cross-section. In the shock rest frame, the direction of the flow velocity is negative along the -direction as shown in Fig. 3; most of photons are then scattered to the downstream, and the peak of the directional angle appears in the backward along the -direction.
Now, the relation between the mean free path and photon traveling distance during is expressed as follows:
where is the mean free path of a photon traveling opposite to the flow-velocity direction. Here, the free path of the photon in the CMF is transformed to that in the OBF by the same transformation as Eq. (15), and is the shortest mean-free-path when . We employed for obtaining the strictest constraint for . This prescription is set from the idea that the optical depth is invariant in any inertial frame. The free path decreases as the shock speed increases due to the Lorentz contraction, and the required is shortened.
In Fig. 4 (a), the result with is significantly different from that with , and the result with is slightly different from that with . On the other hand, the result with is in good agreement with that with . Thus, the result of the angular distribution in the frame is converged at . The relative differences between two distributions with different values in the frame are shown in Fig. 5 (a). Here, the angular distribution with is regarded as a sufficiently converged one. The relative difference between the result with and that with is significant and the relative difference between the result with and that with is smaller, while the relative difference between the result with and that with are almost negligible. The result is significantly different between and because the condition of implies that the mean free path is almost the same as the photon traveling distance during . Hence, photons may not be scattered before updating the flow-field due to insufficiently small . That is, can be regarded as a threshold condition on whether photons are scattered or not scattered before the flow-field update. Indeed, with smaller than unity, the result approaches to the converged one as shown in Fig. 4. Here, corresponding to is s in the frame. The time step for photon transport is comparable to or smaller than for updating the flow-field.
Similarly, in the frame, the result is converged at , as shown in Fig. 4 (b). Here, corresponding to is s in the frame. This implies that the applicable in the frame is smaller by one order of magnitude than that in the frame because in the frame is smaller by one order of magnitude compared to that in the frame as a result of the Lorentz contraction. The angular distribution results in Fig. 4 were computed in two different inertial frames; however, these were transformed to the shock rest frame, so that similar features were obtained. The differences between two distributions with different values in the frame are also shown in Fig. 5 (b), and the difference between the result with and that with is smaller than that of the other two cases and is almost negligible.
These differences with different values are interpreted as follows. The simulation result is never dependent on in the shock rest frame, as shown in Fig. 6. In the shock moving frames, however, photons can escape from the computational domain or overtake the shock wave front before scattering at incorrect timing when is not sufficiently small because the downstream-side and upstream-side computational boundaries and the cell interface (i.e., the shock front) move with the shock speed. As the difference between the boundary moving velocity (that is, shock front moving velocity) and the speed of light along the -axis is small because the optical depth is not so large, the large fails to produce accurate simulation results. The computational grid catches up with the photons or the photons catch up with the grid before the photons complete their free path due to a large ; consequently, the angular distribution varies according to the value of .
In both the and frames, with a suitable time interval, the value of is less than . Therefore, we should adopt that resolves the mean free path to almost ten steps. Now, the value of is not a critical value, but rather a sufficient value, because we change by one order of magnitude in the range of – for examining whether the angular distribution is converged. Examination with a smaller change of in the range from to is required for determining the critical value of .
4.4 Comparison of spectra in different inertial frames
The energy spectra of photons escaped from the computational domain are shown in Fig. 7 (a) in each inertial frame with sample particles. Here, ’s in the shock moving frames are set as , as in the earlier section. The flow velocity at the downstream side of the shock wave increases along the -axis as the apparent shock-Lorentz-factor increases, as shown in Fig. 3. Consequently, the peak energies of the emission spectra are shifted to the high-energy side by the Doppler effect. The energy spectra transformed from each frame to the shock rest frame are shown in Fig. 7 (b). The spectra are identical among the three frames after transformation, and they are converged with satisfying . The second peak in the high-energy side is formed through bulk Compton scattering, in which the photon energy is boosted by collisions with relativistic electrons. The energy of the scattered photon is higher by a square of the flow Lorentz factor compared to that of the incident photon (59). Here, the Lorentz factor of the flow velocity in the upstream side of the shock wave in the shock rest frame, , is . Since the first peak energy is at a few times of keV in Fig. 7 (b), the second peak energy must be located at a few times of , that is, a few times of keV, which is in good agreement with the second peak position in Fig. 7 (b). GeV-order photons are also found in the spectrum because the K-N cross-section and the energy loss by the Compton scattering is not included here and energy cut-off does not appear.
Note that we pick photons that escaped from all the computational boundaries, not only the upstream side of the shock wave but also the downstream side and the boundary parallel to the cylinder axis. As we identify whether the simulation results are in good agreement among all the different frames, all photons treated in the computation should be sampled. In fact, if the observer is at a certain position, photons traveling in the direction opposite to the observer can never be observed. Only photons that enter the view angle of the stationary observer should be sampled to obtain the spectra measured in the observations; however, the number of sampled photons decreases, and large statistical errors may appear in the MC computation. We discuss about this issues in A.
4.5 Effect of shock-wave position
So far, we have considered that the shock wave does not move during . The shock speed is assumed to be constant. Thus, in this section, when photons travel across the shock front during , the true position of the shock front can be calculated analytically with . Therefore, we can correctly determine whether photons can overtake the shock wave since we can compute the true shock position. Similarly, we can correctly determine whether photons escape from the computational boundary because we can compute the true boundary position (which moves at the same speed as the shock wave). This additional technique is required for calculating the scattering process with the true shock position; thus, a more accurate energy spectrum can be obtained with this technique.
Let us assume that a photon, initially located at a position , catches up with the grid (which moves at a uniform speed ) after a time duration . The relation between the positions of the photon and of the moving grid are shown in Fig. 8, where denotes the position of the grid (cell interface) at an initial time. After a time interval , the interface moves to a new location . The distance that the photon travels in the -direction during can be computed as , where is the -component of the photon directional unit vector.
We split a time interval as , where . After a photon travels during a time interval , it stops before the cell interface; otherwise, the mean free path cannot be accurately computed. This is because the local conditions along the photon path (e.g., the density or the flow velocity) might not be known if the photon overtakes the cell interface. This procedure prevents the artificial escape of photons at incorrect time from the computational domain and avoids that photons overtake the shock front at false locations.
The spectra obtained after transforming to the shock rest frame are shown in Fig. 9. We used sample particles, and values in the shock moving frames are set to satisfy as in the previous section. The energy distributions in the second peak of the high-energy side () computed in the three different frames are in much better agreement than those displayed in Fig. 7 (b). In a coupled computation of radiative transfer with hydrodynamics, it is not possible to accurately calculate the position of the shock front since the shock speed may develop with time. Thus, the time-step size in coupled computation should be set sufficiently small to minimize the errors in determining the location of the photon when it crosses the shock front.
Here, since the constraint in the time step is caused by the shock displacement, it should not be considered in the shock rest frame because the position of the shock front is always fixed. Indeed, the comparison among the different ’s in the shock rest frame is shown in Fig. 6, and the angular distributions with any ’s are in good agreement each other.
4.6 Effect of Compton scattering
Here, we also examined the computation with Compton scattering by employing the K-N cross-section and considering energy loss for the photons with energy greater than 10 keV, and the cross-section for Thomson scattering was employed for the other photons. As in the earlier sections, we used sample particles. The time interval in the shock moving frames are set as satisfying , and the true position of the shock front is calculated when a photon overtakes the shock wave, as in the Sec. 4.5. The spectrum after transforming to the shock rest frame with Compton scattering is shown in Fig. 11. The simulation results of all the frames agree each other even with Compton scattering. The cut-off appears in the high-energy side, in contrast to the case in which only Thomson scattering is considered, because the high-energy photon loses its energy through the Compton scattering process. An increase in the number of photons in the range of a few times of results from the shift of the down-scattered photon losing its energy down to the electron rest-mass energy. The readers can refer to B for the details of the peak found at the electron rest-mass energy.
The angular distribution with Compton scattering transformed to the shock rest frame is shown in Fig. 11. The property that a peak appears in the backward direction at along the -axis is similar to that in the case in which only Thomson scattering is considered. The angular distributions after transforming to the shock rest frame with various values are shown in Fig. 13. The result with is significantly different from the others, and the results with and are in good agreement each other. This means that the angular distribution is converged for , and the value of to obtain a converged result is larger by about one order of magnitude than in the case in which only Thomson scattering is included. This happens because the effect of the K-N cross-section acting in the case in which Compton scattering is considered. The variation in the ratio of the total K-N cross-section, , to the total cross-section for Thomson scattering, , depending on the photon energy is shown in Fig. 13. The cross-section decreases as the photon energy increases. Therefore, the high-energy photons are not scattered as frequent as their low-energy counterparts. Thus, the mean free path in the case in which Compton scattering is included is larger than that with only Thomson scattering. Consistently, the value of needed for obtaining a converged angular distribution is larger than that with only Thomson scattering (see Eq. (22)).
Moreover, the distribution at the direction of with respect to the -axis increases rapidly in Fig. 13. The flow velocity in the shock rest frame is along the negative direction of the -axis, and the photons traveling along the direction of the flow (in this situation, the angle between the flow velocity and the direction of photon travel is small) have high energy, as indicated in Eq. (6) for . Therefore, the high-energy photons in the direction of with respect to the -axis are not greatly scattered because of the energy dependence in the K-N cross-section, and they can escape from the computational domain while maintaining the direction. This is why the angular distribution increases rapidly near the direction of . In order to explore the relation of the direction with the energy of photons, the spectra are plotted for photons in the directions of and with respect to the -axis, as shown in Fig. 15, where a bin of the direction angle is . Photons with energy greater than a few MeV do not appear in the cases of and , while many high-energy photons appear at MeV in the case of .
We also explored the above results through an analytical approach to derive Figs. 13 and 17 from the analytical formula of the K-N cross-section and Eq. (24). The relation between the direction angle of the incident photon, , and scattered photon, , is illustrated in Fig. 15. Moreover, the variation of the scattered photon energy in the CMF, , with is shown in Fig. 17. The energy of the incident photon, , is set at the first peak energy in Fig. 11, i.e., 3.6 keV. The critical value of , for which the scattered photon energy in the CMF steeply jumps, corresponds to , and energy after the jump exceeds keV. Photons with direction angles greater than this value are rarely scattered since the cross-section of the photon that has energy exceeding keV is close to zero, as shown in Fig. 13. The critical angle is almost in agreement with the angle at which the angular distribution starts to rapidly increase in Fig. 11, which is found as .
5 Effect of smeared shock front
The effect of the spatial resolution of background flow-field on the radiative transfer computation was discussed in this section. The shock front is inevitably smeared in hydrodynamics computation due to numerical diffusion; then, we examined the significance of the effect on the emitted spectrum. Radiative transfer calculations were implemented for a flow-field with a discontinuous shock wave and for one with an artificially smeared shock wave for comparison.
5.1 Numerical modeling
The setting of the computational domain for the flow-field with the smeared shock front is shown in Fig. 17, where value indicates the initial optical depth. Computations were executed on a one-dimensional system with uniform computational cells in the -direction. The size of the computational domain in the -direction is , which corresponds to . All photons are initially placed at a single point immediately behind the shock wave. The initial position of the shock front corresponds to . Only photons initially emitted in the forward direction along the -axis in the shock rest frame are employed, and the others are omitted. Radiative transfer computations were implemented on the flow-field with both the discontinuous and smeared shock front. In this computation, the inertial frame of interest is the rest frame for the upstream flow. Only photons that escape from the forward boundary, that is, the photons that can overtake the shock front, were sampled. The smeared shock front is set in the shock rest frame by determining the density distribution as follows:
where and are the position of the shock wave at a certain time, , and shock width, respectively. The maximum value of the density, , and the minimum value, , are associated by the Rankine–Hugoniot relations and are set as the downstream and upstream values shown in Table 1, respectively. The distribution of the flow velocity is determined as satisfying the equation of continuity, , where is the density in the CMF and is the four-velocity in the shock rest frame. The equation of continuity can be transformed as follows:
where and are the Lorentz factor of the shock velocity and the three-velocity, respectively. Therefore, the equation is simply expressed with the density in the shock rest frame, , and the three-velocity, , without the density transformed to the one in the CMF. The shock speed corresponds to . The initial distribution of density and flow velocity in the shock rest frame with a discontinuous shock wave are shown in Fig. 18 (a), and those with a smeared shock wave are shown in Fig. 18 (b) with the shock width corresponding to 4-computational-cell length.
5.2 Result of spectrum
We employed sample particles and set as satisfying in the situation with of the earlier section; that is, s. The spectrum sampling the initially emitted photons and the one sampling the photons escaped from the computational domain with the discontinuous shock wave are shown in Fig. 19 (a), and those with the smeared shock wave are shown in Fig. 19 (b). The smearing width of the shock wave is set to a 4-cell length. Both the spectra sampling the escaped photons are different from the initial spectra. The high-energy photons of 1–10 MeV order form a power-law curve whose index is positive in the case of the discontinuous shock wave, while a power-law curve whose index is negative is formed in the case of the smeared shock wave. High-energy photons over 10 MeV disappear in both the cases since only photons that escape from the forward boundary are sampled. Photons obtaining high energy via the bulk Compton scattering are scattered to the backward boundary along the -axis, and the photon with energy greater than 10 MeV is rarely scattered any more since the scattering cross-section is small, as shown in Fig. 13. Therefore, high-energy photons escape from the backward boundary and cannot contribute to the spectra in these cases.
The spectra with the discontinuous shock wave and smeared shock waves with shock widths of 2-cell, 4-cell, and 6-cell length are shown in Fig. 20. The spectrum in the case with the smeared shock wave has the negative-index power-law in the high-energy side and becomes slightly steeper as the shock width increases because energy obtained via bulk Compton scattering decreases. We discuss this effect in the next section.
Actually, it is difficult to use so many computational cells ( cells in the computational domain of cm) in practical multi-dimensional hydrodynamics simulations. We used such a large number of computational cells to show that the shock structure affects the emitted spectrum even when such a fine computational mesh is used. If we have 10 or 100 times worse resolution, the spectral changes depending on the shock width in the high-energy side may be more remarkable.
5.3 Effect of shock width
We explored the relation between the shape of the spectra and the shock width. The initial distributions of optical depth per cell, , and the Lorentz factor of the flow velocity in the OBF with the discontinuous shock wave are shown in Fig. 21 (a), and those with the smeared shock wave are shown in Fig. 21 (b). The optical depth is an indicator of the probability of scattering, and the Lorentz factor is an indicator of energy obtained via bulk Compton scattering. For obtaining energy via bulk Compton scattering and traveling toward the observer, photons should be emitted at the downstream side of the shock wave, then scattered at the upstream side to the backward direction, and scattered again at the downstream side to the forward direction. Therefore, the optical depth is required to be rather large around the flow velocity jump. The optical depth in the case with the discontinuous shock wave is greater than that in the case with the smeared shock wave immediately behind the velocity jump. Photons are rarely scattered at the downstream side in the case with the smeared shock wave, in contrast to the case with the discontinuous shock wave; consequently, high-energy photons traveling to the forward boundary decrease.
Since the structure of the shock wave can considerably affect the radiative transfer computation in even the one-dimensional background flow-field, we should also consider the effect in the radiation-hydrodynamics coupled computation. We implemented the coupled computation neglecting the radiation back reaction and discussed in C about the difference of the spectral feature from that in the model computations shown in this section.
Radiative transfer computation in the ultra-relativistic flow-field was validated as the preliminary study for the coupled computation of the radiative transfer and relativistic hydrodynamics with the appropriate simulation conditions in the present paper. We have developed a three-dimensional radiative transfer code for the ultra-relativistic background flow-field by using the MC technique, and some test calculations were implemented.
We performed the radiative transfer computations in three different inertial frames in which the apparent Lorentz factor of the shock speed, , is 1, 10, and 100 and compared the simulation results in the same frame. The angular distributions computed with different values showed that the simulation result was converged with sufficiently small . The value of that resolves the mean free path into ten steps was the sufficient condition for obtaining thoroughly converged results. The value of tested in this paper is appropriate for only the case without including the radiation back-reaction on the flow-field. In the case with radiation back-reaction, the time step for updating the flow-field will be smaller. In a fully coupled radiation hydrodynamics simulation, energy-momentum is extracted out of hydrodynamic cells and carried away by photons. Thus, the energy (or pressure) of the cells may become negative unless a sufficiently small time step is employed. In optically thick regions of the flow-field, this restriction can be so tight and an implicit or partly-implicit time integration algorithm would be needed.
The spectra computed in each inertial frame showed that the peak energy was shifted to the high-energy side as the flow velocity increases because of the Doppler effect. The spectra transformed from each frame to the shock rest frame were in good agreement, validating the transformation among the different frames in the ultra-relativistic regime. The spectra were in better agreement on considering the true position of the shock front when photons crossed the shock wave. The spectra with Compton scattering had peak energy at a few times of because photons were scattered and lost their energy down to the electron rest mass energy. The angular distribution with Compton scattering showed a rapid increase along the direction of with respect to the -axis. The K-N cross-section decreases with the increase in the photon energy, and the photons traveling toward the direction of had high energy; then, these photons escaped from the computational domain without scattering any more. With Compton scattering, for obtaining the converged angular distribution was larger than that in the situation with only Thomson scattering because of the energy dependence of the K-N cross-section. Therefore, the constraint of may be relaxed in the case that both of Thomson and Compton scatterings are considered.
We assessed the effect of the spatial resolution of the background flow-field on the radiative transfer computation. Radiative transfer was computed in scenarios with a discontinuous and an artificially smeared shock fronts, and the simulation results were compared. High-energy photons with 1–10 MeV order decreased in the scenario with the smeared shock front compared to that with the discontinuous shock front. The optical depth immediately behind the shock wave was smaller with the smeared shock front compared to that with the discontinuous shock front, and photons that obtain high energy through bulk Compton scattering and travel to the observer decreased. The structure of the shock wave affected the radiation transport even in the one-dimensional computation. In the one-dimensional models presented here, we have employed hydrodynamical cells to accurately compute the high-energy tail of the spectrum; however, it is difficult to employ so many cells per dimension in multi-dimensional simulation. Extending the method proposed here to more than one dimension is handicapped by the huge numerical resolution. This is a future work toward the multi-dimensional computation.
Some test calculations were performed for validating the radiative transfer computation on the ultra-relativistic background flow-field as the preliminary study for the coupling computation of MC radiation transport with relativistic hydrodynamics. In the future, we try to combine the developed code in this paper with relativistic hydrodynamics simulation to achieve radiative transfer computation on the time-dependent ultra-relativistic flow-field involving the radiation feedback to the flow-field and firstly perform some test computations restricted to some regions in the flow-field.
We are grateful to the anonymous referees for their so fruitful comments on this manuscript.
This work was supported by JSPS KAKENHI Grant Number 52638590,
a Grant-in-Aid for Scientific Research from the Ministry of
Education, Culture, Sports, Science, and Technology (MEXT) of Japan (24103006, 24740165, and 24244036),
and the HPCI Strategic Program of the Japanese MEXT.
This work is partly supported by the Grant-in-Aid for Scientific Research (S:16H06341) and
the Grant-in-Aid for Young Scientists (B:16K21630) from the MEXT of Japan.
H. N. is supported in part by JSPS Postdoctoral Fellowships for Research Abroad No. 27–348,
and also supported at Caltech through NSF award No. TCAN AST–1333520.
H. I. is supported in part by the RIKEN pioneering project âInterdisciplinary Theoretical Science (iTHES)â.
Appendix A Convergence test for only photons traveling to the observer
We computed the spectrum equivalent to Fig. 7 but only for photons that would reach an observer. The observer was assumed to be at the upstream side of the shock wave and far enough away from the shock front. Therefore, we sampled only photons escaped from the forward boundary of the computational domain. That is, photons moving at an angle with respect to the shock normal propagation direction were sampled. The spectrum result has larger statistical errors than that for all photons. We compared the result with that for 10 times larger number of photons ( particles) and smaller number of photons ( particles) as shown in Fig. 22. The spectrum for particles has large statistical errors in the high-energy side, and the convex feature is not clear in the range of keV. Although the spectrum for particles has somewhat statistical errors, the convex feature can be found in the high-energy side. The spectrum for particles is better converged. Thus, more than particles may be needed for obtaining a realistic observed spectrum.
Appendix B Energy peak formed by Compton scattering
The peak in the range of in Fig. 11 is formed due to the down-scattered photons from the high-energy side. The separated spectra depending on the photon traveling path across the shock wave are shown in Fig. 23. The word of ‘stay’ in the legend denotes the spectrum for the photons staying in the downstream side from emission to escape out from the computational boundary. Similarly, the word of ‘one way’ denotes the spectrum for the photons traveling from the downstream side to the upstream side crossing the shock wave once, ‘round trip’ denotes that the photons travel from the downstream side to the upstream side and subsequently from the upstream side to the downstream side crossing the shock wave twice, and ‘others’ denotes the other photons. The spectrum of ‘stay’ is not significantly different from the initial emitted spectrum; however, that of ‘one way’ is shifted to the high-energy side due to bulk Compton scattering. Moreover, the spectrum of ‘round trip’ has the peak at since the high-energy photons produced in the upstream side are carried back to the downstream side and Compton scattered several times losing their energy down to the electron rest mass energy.
The photon energy after scattering in the CMF, , is calculated as follows:
The second term in the denominator in the right-hand side is predominant when the photon energy is larger than the electron rest mass energy, ; consequently, the photon energy after scattering decreases. In contrast, when the photon energy is smaller than , this term becomes negligible and the photon energy after scattering is not noticeably changed from the incident energy. Therefore, the photon with high energy loses its own energy through several scattering processes down to , and if its energy falls less than , the photon energy is unchanged any more.
The flow velocity in the downstream side in the shock rest frame is not highly relativistic, so the photon energy corresponding to the electron rest mass in the CMF is not so different from that in the shock rest frame. Therefore, most of the photons that have high energy after experiencing bulk Compton scattering in the upstream side and carried to the downstream side hold energy of .
In order to show that the peak energy is located at under any condition, the spectra with different temperatures are shown in Fig. 24. Here, denotes criterial temperature. The position of the first peak energy in the left hand is shifted to the higher energy side by one order of magnitude as the temperature increases by one order of magnitude (from (a) to (c)). However, the peak of remains since the peak position is determined by the photon energy corresponding to the electron rest mass, and the photon energy is unchanged as long as being calculated in the same inertial frame.
Appendix C Post-processed radiative transfer computation of relativistic hydrodynamics simulation
It is challenging task to reproduce the shock wave with in numerical simulation with any numerical schemes. So, for developing the relativistic radiation hydrodynamics code, it is necessary to develop not only the MC radiative transfer code but also hydrodynamics simulation code which is appropriate for high ; however, it is out of the scope for this article, and we show the result of radiative transfer computation on the relativistic hydrodynamical flow-field with only as a guide.
We performed post-processed computation of the MC radiative transfer on the time-dependent flow-field with relativistic hydrodynamics computation but no feedback from interaction of photon with flow-field matter. Relativistic hydrodynamics simulation code was the same as that in the previous work (17). Hydrodynamics simulation was performed in advance, and radiative transfer was computed as a post-process. The computation was conducted on a one-dimensional system with uniform computational cells in the -direction using sample particles for comparison with the model computations presented in Sec. 5. The initial conditions were the same as shown in Table 1 obtained by analytically solving R-H relations, and the velocities were transformed from the shock rest frame to the shock moving frame with . The shock front was initially located near the left boundary. The computation was terminated when the shock front reaches at the right boundary. The distributions of optical depth and Lorentz factor of the flow velocity in the hydrodynamical flow-field at a certain time of s are shown in Fig. 25. The smeared shock front is similar to that in the modeled flow-field as shown in Fig. 18 (b).
Fig. 26 shows the photon spectrum produced for two types of models. Models labeled with “model with discontinuous shock” and “model with smeared shock” refer to calculations in which the modeled background flow-field is set up as in Sec. 5. The time interval for updating the flow-field in the model computations are set as satisfying in the situation with in Sec. 4.3; that is, s. The other spectrum of Fig. 26 corresponds to post-processed radiative transfer computation of relativistic hydrodynamics simulation in which the background flow-field is evolved employing the hydrodynamical time step ; namely, s. We note that the value of the hydrodynamical time step is larger by one order of magnitude than in the model computations shown in Sec. 5. The reason to set such large value of is to limit the computational cost. The spectrum in the high-energy side has the positive-index power-law in the model computation with the discontinuous shock front, while the negative-index power-law appears in the model computation with the smeared shock front. The negative power-law slope becomes steeper as the shock width increases because bulk Compton scattering rarely occurs as mentioned in Sec. 5.3. In the hydrodynamics computation, since the shock width becomes wider due to numerical diffusion as time advances, the high-energy photons further decrease. The spectrum obtained from the post-processed simulation has the small number of photons in the high-energy side, and this feature is consistent with the model computations.
In the numerical simulation, the flow-field cannot keep a single shock structure even when the initial conditions obtained by analytically solving R–H relations are employed. The R–H equations are the relation that links between the upstream and downstream sides of the shock wave. So, if the flow-field structure is calculated analytically with the same initial conditions, it can keep a single shock structure. This difference between analytical solution and numerical one is caused by using the approximate Riemann solver in the numerical simulation. Therefore, the flow-field evolves to the Riemann problem as shown in Fig. 25. The realistic flow-field developing with time cannot be represented by the simple model flow-field. So, the shock structure in our model flow-field corresponds to a snapshot of the flow-field computed by hydrodynamics computation at an early time. For the above reason, there is the lag between two jumps of the optical depth and the velocity as shown in Fig. 25. The optical depth goes up in the region which already reaches the terminal (downstream) velocity where the velocity does not significantly change. So, a large number of photons cannot experience the difference in the bulk Lorentz factor in the two scattering positions. This should kill lots of the high energy photons produced by the bulk Compton scattering which can be produced in the model computation with the discontinuous shock front.
Therefore, radiation hydrodynamics computation with coarse computational grids resulting in the expanded shock width and with the lag between the optical depth and the velocity jumps cannot produce a lot of high-energy photons due to less bulk Compton scattering. The effect of the simulation condition on the spectrum is sustained in the result of the test calculation performed in this paper.
- journal: Journal of Computational Physics
- G. C. Pomraning, The Equations of Radiation Hydrodynamics, Pergamon Press, California (1973).
- G. M. Webb, Astrophys. J. 296 (1985) 319.
- J. G. Kirk, D. B. Melrose, and E. R. Priest, Plasma Astrophysics, Springer-Verlag, Berlin (1994).
- J. L. Anderson and E. A. Spiegel, Astrophys. J. 171 (1972) 127.
- N. Tominaga, S. Shibata, and S. I. Blinnikov, Astrophys. J. Suppl. Ser. 219 (2015) 38.
- M. J. Rees and P. Meszaros, Astrophys. Lett. 430 (1994) L93.
- R. Sari and T. Piran, Astrophys. J. 485 (1997) 270.
- S. Kobayashi, T. Piran, and R. Sari, Astrophys. J. 490 (1997) 92.
- P. Mimica and M. A. Aloy, Mon. Not. R. Astron. Soc. 401 (2010) 525.
- P. Mimica and M. A. Aloy, Mon. Not. R. Astron. Soc. 421 (2012) 2635.
- D. Eichler and A. Levinson, Astrophys. J. 529 (2000) 146.
- P. Meszaros and M. J. Rees, Astrophys. J. 530 (2000) 292.
- M. A. Aloy et al., Astrophys. J. Lett. 531 (2000) L119.
- M. A. Aloy et al., Astron. Astrophys. 396 (2002) 693.
- W. Zhang, S. E. Woosley, and A. I. Macfadyen, Astrophys. J. 586 (2003) 356.
- A. Mizuta et al., Astrophys. J. 651 (2006) 960.
- H. Nagakura et al., Astrophys. J. 731 (2011) 80.
- J. Matsumoto and Y. Masada, Astrophys. J. Lett. 772 (2013) L1.
- D. Lazzati, B. J. Morsony, and M. C. Begelman, Astrophys. J. 700 (2009) L47.
- A. Mizuta, S. Nagataki, and J. Aoi, Astrophys. J. 732 (2011) 26.
- C. Cuesta-Martinez, M. A. Aloy, and P. Mimica, Mon. Not. R. Astron. Soc. 446 (2015) 1716.
- C. Cuesta-Martinez et al., Mon. Not. R. Astron. Soc. 446 (2015) 1737.
- M. S. Briggs et al., Astrophys. J. 524 (1999) 82.
- H. T. Janka, Nuclear Astrophysics: Proceedings of the Workshop 287 (2005) 319.
- H. T. Janka and W. Hillebrandt, Astron. Astrophys. Suppl. Ser. 78 (1989) 375.
- H. Th. Janka, Astron. Astrophys. 256 (1992) 452.
- M. Th. Keil and G. G. Raffelt, Astrophys. J. 590 (2003) 971.
- E. Abdikamalov et al., Astrophys. J. 755 (2012) 111.
- K. Maeda, Astrophys. J. 644 (2006) 385.
- A. Suzuki and T. Shigeyama, Astrophys. J. 719 (2010) 881.
- L. B. Lucy, Astronomy & Astrophys. 429 (2005) 19.
- A. M. Beloborodov, Astrophys. J. 737 (2011) 68.
- A. Ishii et al., High Energ. Dens. Phys. 9 (2013) 280.
- S. Shibata, N. Tominaga, and M. Tanaka, Astrophys. J. Lett. 787 (2014) L4.
- A. Ishii et al., High Energ. Dens. Phys., (2014).
- F. Ryde and A. Pe’er, Astrophys. J. 702 (2009) 1211.
- C. C. Thöne et al., Nature 480 (2011) 72.
- A. Pe’er and F. Ryde, Astrophys. J. 732 (2011) 49.
- D. Giannios and H. C. Spruit, Astron. Astrophys. 469 (2007) 1.
- D. Giannios, Astron. Astrophys. 480 (2008) 305.
- K. Nakayama and T. Shigeyama, Astrophys. J. 627 (2005) 310.
- R. Takagi and S. Kobayashi, Astrophys. J. 622 (2005) L25.
- Y. Ohtani, A. Suzuki, and T. Shigeyama, Astrophys. J. 777 (2013) 113.
- H. Ito et al., Astrophys. J. 777 (2013) 62.
- H. Ito et al., Astrophys. J. 789 (2014) 159.
- U. M. Noebauer et al., Mon. Not. R. Astron. Soc. 425 (2012) 1430.
- H. Ito et al., Astrophys. J. Lett. 814 (2015) L29.
- N. Roth and D. Kasen, Astrophys. J. Suppl. Ser. 217 (2015) 9.
- P. Mimica et al., Astron. Astrophys. 441 (2005) 103.
- J. A. Fleck JR. and J. D. Cummings, J. Comput. Phys. 8 (1971) 313.
- R. G. McClarren and T. J. Urbatsch, J. Comput. Phys. 228 (2009) 5669.
- N. A. Gentile, J. Comput. Phys. 172 (2001) 543.
- R. T. Wollaeger et al., Astrophys. J. Suppl. Ser. 209 (2013) 36.
- M. A. Cleveland and N. Gentile, J. Comput. Phys. 291 (2015) 1.
- J. D. Densmore and E. W. Larsen, J. Comput. Phys. 291 (2015) 1.
- J. D. Densmore et al., J. Comput. Phys. 284 (2015) 40.
- W. Heitler, The Quantum Theory of Radiation, Dover Publications, New York (2010).
- B. Schutz, A First Course in General Relativity, Cambridge University Press, New York (2009).
- G. B. Rybicki and A. P. Lightman, Radiative Processes in Astrophysics, Wiley-VCH, Mörlenbach (1985).