A novel, numerically stable co-moving method to speed up line radiative transfer.

Thomas Ceulemans 1

  • 1 Ku Leuven, Leuven

Abstract

In Asymptotic Giant Branch (AGB) stars, complex morphologies of the stellar winds have been observed. These can be modelled using hydrodynamical simulations. In order to bridge the gap from hydrodynamical simulations towards observations, radiative transfer modelling is required. As light interacts with the circumstellar matter, observations will contain information about the chemical composition, geometry and hydrodynamical properties of the stellar wind.

It is well known that Doppler shifts will shift the frequencies at which spectral lines appear. However, from a computational point of view, this provides a challenge for non-local thermodynamical equilibrium (NLTE) line radiative transfer as this requires the intensities around the line centers to be computed self-consistently with the properties of the medium.

The line radiative transfer code MAGRITTE computes intensities using rays traced throughout a 3D model. Traditionally, the radiative transfer equation is solved using a fixed frequency discretisation around the line centers. Due to Doppler shifts, this discretisation changes for each point on the ray. This weighs heavily on the computation time, as the frequencies for which one needs to compute the intensity are not aligned. As 3D radiative transfer calculations are computationally challenging even today, reducing the amount of computation time is paramount.

By using a numerically stable modification of the radiative transfer equation, we have rewritten the radiative transfer equation in order to cope with Doppler shifts. In this manner, we can shift the frequency for which we compute the intensity during the computation. As a result, we are able to reduce the number of frequencies at each point for which one needs to compute the intensity by following the Doppler shift on the ray. To stabilise the numerical algorithm, we introduce a novel and intuitive insight into numerical mathematics.

In order to show the effectiveness of the new algorithm, we apply it to 3D hydrodynamical models, comparing it to established radiative transfer algorithms. In particular, we test the solver using PHANTOM hydrodynamical simulations of AGB binary systems. These AGB models exhibit a complex velocity field, which is handled by an accurate implementation of boundary conditions in frequency space, enabling it to handle non-monotonic velocity gradients on a ray. The speedup obtained from the new solver brings us closer to being able to include consistently radiative line cooling in hydrodynamical simulations, instead of radiative transfer only being used to post-process these hydrodynamical models.

Co-moving method

In line radiative transfer, the intensity at a specific frequency changes when that frequency is in the neighborhoud of an atomic/molecular line. Due to Doppler shifts, this region moves in frequency space if the velocity changes.

Using the default radiative transfer equation on a ray,

\frac{\partial I(x,\nu)}{\partial x}=\eta(x,\nu)-\chi(x,\nu) I(x,\nu)
\text{In which } I \text{ is the intensity, } \eta \text{ is the emissivity, } \chi \text{ is the opacity, } x \text{ is the position and } \nu \text{ is the frequency.}

does not allow for a change in the frequency 𝜈. Thus computation time will be spent computing intensities in frequency regions where the intensity does not significantly change.

For computational efficiency, we propose a slight alteration of the radiative transfer equation.

\frac{d I(x,\nu)}{d x}=\eta(x,\nu)-\chi(x,\nu) I(x,\nu) \boldsymbol{+ \frac{d\nu}{dx}\frac{\partial I}{\partial \nu}(x,\nu)}
\text{In which } I \text{ is the intensity, } \eta \text{ is the emissivity, } \chi \text{ is the opacity, } x \text{ is the position and } \nu \text{ is the frequency.}

This formulation allows for a change in frequency when solving this equation.

I exploit this feature in two applications:

  • Computing NLTE intensities for line radiative transfer: The radiative properties of the medium depend on the radiation field at frequencies near the (Doppler shifted) line frequencies.
  • Imaging models with non-monotonic velocity gradients: Both a wide frequency range and a fine frequency distribution is required for the line profile.

This comoving description of the radiative transfer equation similar to Baron & Hauschildt 2004 (and others in the same series). This poster adds the following improvements:

  • A novel way to stabilize the equations
  • Implementation of accurate boundary conditions which support non-monotonic velocity gradients

Numerical stability

The core idea of this section can be summarized in the following statement: extrapolation causes numerical issues.

\text{To practically use the comoving solver, we need to define the change in frequency } \nu(x) \text{.}

As the interesting frequency regions for line radiative transfer lie around the Doppler shifted line centers, an initial idea is for the frequencies to follow the Doppler shift, as illustrated in the figure on the right.

For line radiative transfer, we want to resolve details around the line center, thus a fine frequency binning is required. This configuration causes numerical issues when encountering velocity gradients:

\text{We need to compute a frequency derivative term }\frac{dI}{d\nu} \text{. If the frequency binning is too fine, or the velocity gradient too large, numerical issues appear.}

A 1D example of the differences between the default and the stabilized comoving solver. The model has an equidistant frequency discretization, a constant velocity gradient, a constant source function S = 1, and constant opacity χ = 1.

Instead of letting the frequency change match the Doppler shift, we minimize the frequency difference.This stabilizes the solver, as the frequency derivative is no longer extrapolated.

Application: NLTE radiative transfer

In 3D NLTE line radiative transfer, the radiation field itself impacts the properties of the medium. This requires an iterative procedure computing mean line intensities J(x).

J_{ij}(\boldsymbol{x})=\oint_{\Omega}\int_0^{+\infty} I(\hat{\boldsymbol{n}},\boldsymbol{x},\nu)\phi_{ij}(\nu)\text{d}\Omega \text{d}\nu

We compare the required computation time for a single NLTE iteration by applying it to a Phantom model (Malfait et al. 2021).

  • Feautrier: standard feautrier solver (Feautrier 1964), as implemented in Magritte
  • Comoving (non-local bdy): Comoving method, using an accurate implementation of the boundary conditions in frequency space
  • Comoving (local bdy): Comoving method, using blackbody intensities everywhere for the boundary conditions in frequency space (simplifies implementation)

Slice plot of the Phantom binary wind model by Malfait et al. 2021

Application: Imaging

In regions with large velocity gradients, the resulting line profile is wide in frequency space. For detailed line radiative transfer, this implies that a fine frequency discretization is required to properly resolve the line profile.

When using a static solver (for which the frequency does not change), the required computation time will scale linearly with the amount of requested frequencies in the synthetic observation. This wastes some computation power, as every static frequency only lies close to a line center for a close part of the model, due to large velocity gradients resulting in significant Doppler shift.

Using the comoving solver, we first compute intensities with a fine frequency grid around the Doppler shifted line centers. Afterwards, the computed intensities are interpolated, correctly taking into account the boundary conditions. and correctly interpolated. In this way, the amount of frequencies requested for the line profile has no impact on the computation time.

We apply this to a 2D AMRVAC model (Moens et al. 2022, Debnat et al. in prep) of a Wolf-Rayet wind with highly non-monotonic velocity gradients.

Synthetic line profile of the CO 1-0 line of the WR model, comparing the new comoving method with the established Feautrier method. The comoving method agrees with the Feautrier method, albeit the line profile is smoothed out.

The highly non-monotonic velocity profile of the 2D WR AMRVAC model.

Density profile of the 2D WR AMRVAC model.

Temperature profile of the 2D WR AMRVAC model.

Boundary conditions

When using the comoving solver for line radiative transfer, one requires boundary conditions not only in position space (at the start of a ray), but also in frequency space. This is because Doppler shifts change the relevant frequency range for which the intensity should be computed.

If one encounters a specific frequency range for the first time, blackbody intensities are sufficient as boundary condition. However, any subsequent time that frequency range is needed, we have to use the previously computed intensities at that frequency.

For more details, come talk to me at the conference.

Choosing 𝜈(x)

Solid line: static; dashed line: following Doppler shift; dotted line: minimal frequency differences

Check out the radiative transfer code Magritte

Personal website