Introduction

TODO

Source term

TODO

Helmholtz approach

TODO

Monte Carlo approach

In our Monte Carlo method for photoionization, we assume that photon scattering can be neglected, and that the direction of ionizing photons is isotropic. The method can be divided in three parts:

• Determine the coordinates at which photons are produced, and store these coordinates in a list \(L_{src}\).
• Determine the coordinates at which photons are absorbed, and store these coordinates in a list \(L_{dst}\).
• Compute the resulting photoionization profile on a mesh.

The implementation of these steps is described below.

The source of ionizing photons

In a plasma fluid simulation, the computational domain is divided into cells. We assume that for each cell the production of ionizing photons \(I\) is known within a given time step \(\Delta t\). In our Monte Carlo method, this information is converted to a list \(L_{src}\) of approximately \(N\) discrete photons in the following way.

1. Determine the total photon production \(I_{\Sigma}\) within the time step \(\Delta t\), by summing over all the cells.
2. For each cell \(n_{\gamma} = N I / I_{\Sigma}\) photons should be produced. To convert \(n_{\gamma}\) to an integer, first draw a uniform random number \(U (0 , 1)\). If \(n_{\gamma} - \lfloor{ n_{\gamma}} \rfloor > U\) round up, else round down.
3. For each produced photon, add the coordinate of the cell center to the list \(L_{src}\) of photons.

Some additional remarks: In principle the number of photons \(N\) can be chosen adaptively, for example to create discrete ‘super-photons’ with a given weight. It is also possible to assign different weights to different photons by modifying the above procedure, see section 11.3.4. Instead of rounding \(n_{\gamma}\) to an integer as described above, it can be more realistic to sample from the Poisson distribution with parameter \(n_{\gamma}\). For \(n_{\gamma} \ll 1\) the result would be almost identical, but for larger values there are differences. However, most of the random fluctuations in our method are due to the stochastic photon direction and absorption, as discussed in the next section. Therefore our rounding with uniform random numbers should be fine for most applications. When grid cells are large compared to typical photoionization length scales, one could determine a ‘subgrid’ production location in the cell, instead of using the cell center.

Absorption of ionizing photons

Now that we know where photons are produced, we need to determine where they are absorbed. This is done in the following way. First, we determine the absorption distance \(r\) for a photon. Given an absorption function \(f(r)\), the cumulative absorption function \(F(r) =\int_0^r f(r')dr'\) can be computed, either numerically or analytically. Then a so-called lookup table is constructed with columns \(F(r)\) and \(r\). Now we can do inversion sampling: given a random number \(U(0, 1)\), the corresponding distance \(r\) is obtained by linear interpolation from the table. A random orientation \(\hat{r}\) for the photon is determined, using a procedure for picking a point on a sphere proposed by Marsaglia [205]:

1. Get two random numbers \(U_1(-1,1)\) and \(U_2(-1,1)\). If \(U_1^2 + U_2^2 \leq 1\) accept them, otherwise draw new numbers.
2. Let \(a = U_1^2 + U_2^2\) and \(b = 2 \sqrt{1-a}\). The isotropic random direction is \(\hat{r} = (b U_1, b U_2, 1 - 2a)\) in Cartesian coordinates. In the special case of a 2D Cartesian \((x, y)\) coordinate system, we should not pick points on a sphere but points on a circle. Then the second step is replaced by

   \[ \hat{r} = ( \frac{U_1^2-U_2^2}{U_1^2+U_2^2},\frac{2 U_1 U_2}{U_1^2+U_2^2}) . \]

   Given the direction and the distance, the location of absorption \(\mathbf{r} = r \hat{\mathbf{r}}\) is known, which is added to the list \(L_{dst}\). Sometimes, the typical absorption length of photons is much larger than the domain size. Since most photons will not be absorbed within the domain, the Monte Carlo approximation becomes less accurate. This can be resolved by changing the lookup table. Suppose that all photons with absorption distances \(r > r_{max}\) can be ignored, then one can construct a lookup table up to \(r_{max}\) and scale the corresponding \(F(r)\) values to the range (0; 1). Each photon that is produced now gets a smaller weight \(F(r_{max})\).

A Lookup table with the absorption locations of the photons has been created. These coordinates need to be mapped to a mesh to get the photoionization profile Two options can be considered: the photons can be absorbed on a single mesh with a constant grid spacing, or they can be absorbed at different levels in a refined mesh. In both cases the nearest grid point (NGP) is used to map the photons to densities, which means that photons are absorbed within a single cell. With linear (and higher order) interpolation the resulting density profiles are smoother, but it becomes harder to handle refinement boundaries.