High Accuracy Photodissociation Rates in VULCAN with MCRT

Monte Carlo radiative-transfer (MCRT) techniques have long been used as when high accuracy is required for complex radiative-transfer calculations such as 3D multiple scattering simulations (e.g., Noebauer & Sim 2019). We develop a 1D MCRT plane-parallel method python module gCMCRT.py accelerated using numba (Lam et al. 2015) as part of the VULCAN photochemical model framework (Tsai et al. 2021).

For 1D plane-parallel simulations using the unweighted Lucy (1999) path length method, the zeroth moment of the intensity, the mean intensity, J_λ (erg s⁻¹ cm⁻² cm⁻¹), for a wavelength, λ (cm), in a layer is given by

$$\begin{eqnarray}&&{J}_{\lambda }=\displaystyle \frac{1}{4\pi {\rm{\Delta }}z{\rm{\Delta }}t}\sum {\epsilon }_{\lambda }l,\end{eqnarray} \tag{ 1 }$$

where Δz (cm) is the layer thickness, Δt (s) the time of simulation (typically 1s), _λ (erg cm⁻² cm⁻¹) the packet energy density and l (cm) the path length of the packet through the layer and μ the cosine angle of the packet.

One of the key quantities in photochemical modeling is the actinic flux, which is the effective number of photons per unit area per unit time per unit wavelength regardless of direction. The original VULCAN model uses the two-stream approximation for the UV radiative transfer calculation to obtain the diffuse component of the actinic flux. The conversion from the total diffuse flux to the required mean intensity is described in Tsai et al. (2021) (Section 2.4). For the MCRT model, the actinic flux is simply related to the mean intensity, F_actinic,λ  = 4πJ_λ . For shortwave radiation, each packet carries a fraction of the incident flux, F_inc,λ (erg s⁻¹ cm⁻² cm⁻¹),

$$\begin{eqnarray}&&\displaystyle \frac{{\epsilon }_{\lambda }}{{\rm{\Delta }}t}=\displaystyle \frac{{\mu }_{\star }{F}_{\mathrm{inc},\lambda }}{{N}_{\mathrm{ph}}},\end{eqnarray} \tag{ 2 }$$

where μ_⋆ is the zenith cosine angle of the stellar irradiation and N_ph the number of packets.

Due to this divisibility, energy is guaranteed to be conserved, since each photon packet carries a proportional fraction of the total energy (Lucy 1999). Packets are initialized at the top of the atmosphere with a direction of −μ_⋆, after which they are free to scatter or be absorbed by the atmosphere. Packets are terminated after being absorbed in the atmosphere or escape through the top of the atmosphere. We include surface Lambertian scattering should surface effects be required. We find only N_ph ∼ 1000+ packets per wavelength band are required to get a good signal of the mean intensity, allowing a quick computational time only a few times slower than two-stream methods.

The advantage of MCRT schemes is that no assumptions are made on the radiation propagation direction and scattering source function, as scattering phase functions are sampled randomly for each scattering event (e.g., Noebauer & Sim 2019). So far we include sampling isotropic, Rayleigh and Henyey–Greenstein phase functions in the model, but more complex phase functions can be easily added in the future.

We foresee the MCRT module being highly useful for complex multiple scattering scenarios such as high altitude clouds/hazes where strong albedos and forward scattering are expected, allowing accurate computation for the photolysis rates (e.g., Michelangeli et al. 1992). Figure 1 shows a possible haze scenario for the default VULCAN HD 189733b hot Jupiter simulation and how it affects the mean intensity profile. This method is generally applicable, and can be adapted for other modeling efforts which also require complex shortwave multiple scattering calculations.

Zoom In Zoom Out Reset image size

E.K.H.L. is supported by the SNSF Ambizione Fellowship grant (#193448). S-M. Tsai is supported by the University of California at Riverside.

Software: VULCAN (Tsai et al. 2021).