compton.cu

Documentation for the main functions that handle Compton scattering.

Functions

Functions

__device__ void sample_scattered_photon(double k[4], double p[4], double kp[4], curandState *localState)

Computes the scattered photon 4-momentum after a collision with an electron. This high-level routine simulates the physical interaction between a photon and a specific electron. It handles the transition between classical Thomson scattering and quantum Klein-Nishina scattering based on the photon’s energy in the electron rest frame.

Routine:

1. Frame Shift: The incoming photon 4-momentum \( k^\mu \) is boosted into the electron’s rest frame.

2. Regime Selection:

• If \( E'_{photon} > 10^{-4} m_e c^2 \): Uses Klein-Nishina sampling for energy loss and angular deflection.

• If \( E'_{photon} \le 10^{-4} m_e c^2 \): Uses Thomson sampling (elastic scattering).

1. Angular Sampling: Constructs a local coordinate system aligned with the incident photon to sample the scattering angle \( \theta \) and azimuthal angle \( \phi \).

2. Inverse Transformation: The resulting scattered momentum is boosted back from the electron rest frame to the laboratory/fluid frame.

Parameters:

• k – Incoming photon 4-momentum in the laboratory frame.

• p – 4-momentum of the scattering electron.

• kp – [out] Scattered photon 4-momentum in the laboratory frame.

• localState – Pointer to the curand RNG state for the current thread.

Returns:

void

__device__ void boost(double v[4], double u[4], double vp[4])

Performs a general Lorentz boost into the frame of a given 4-velocity.

Transforms a 4-vector \( v \) from the laboratory frame into the frame moving with 4-velocity \( u \). The implementation uses the general Lorentz transformation matrix to handle boosts in arbitrary directions.

Parameters:

• v – Input 4-vector in the laboratory frame.

• u – 4-velocity of the frame to boost into (where \( u^0 = \gamma \)).

• vp – Resulting boosted 4-vector in the co-moving frame.

Returns:

void

__device__ double sample_thomson(curandState *localState)

Samples the scattering angle cosine for Thomson scattering.

Uses rejection sampling to pick the cosine of the scattering angle \( \mu = \cos\theta \) according to the Thomson differential cross-section distribution:

\[ \frac{dP}{d\mu} = \frac{3}{8}(1 + \mu^2) \]

Parameters:

localState – Pointer to the curand RNG state for the current thread.

Returns:

The sampled cosine of the scattering angle \( \mu \).

__device__ double sample_klein_nishina(double k0, curandState *localState)

Samples the post-scattering photon energy using the Klein-Nishina differential cross-section.

This function uses a rejection sampling algorithm to determine the scattered photon frequency \( k' \) in the electron rest frame. The sampling occurs within the kinematically allowed range:

\[ \frac{k}{1+2k} \le k' \le k \]

where the lower bound corresponds to a head-on collision ( \( \theta = \pi \)) and the upper bound corresponds to no deflection ( \( \theta = 0 \)).

Parameters:

• k0 – Dimensionless photon frequency in the electron rest frame before scattering.

• localState – Pointer to the curand RNG state for the current thread.

Returns:

The sampled dimensionless photon frequency after scattering \( k' \).

__device__ double klein_nishina(const double a, const double ap)

Computes the differential Klein-Nishina cross-section for a given frequency shift.

This function evaluates the probability density of a photon with dimensionless energy \( a \) scattering to energy \( a' \). It is used as the target distribution for the rejection sampling in the Klein-Nishina scattering routine.

The implementation calculates:

\[ \frac{1}{a^2} \left( \frac{a}{a'} + \frac{a'}{a} - \sin^2 \theta \right) \]

where the scattering angle \( \theta \) is implicitly determined by the Compton formula: \( \cos \theta = 1 + \frac{1}{a} - \frac{1}{a'} \).

Parameters:

• a – Initial dimensionless photon energy ( \( h\nu / m_e c^2 \)).

• ap – Scattered dimensionless photon energy ( \( h\nu' / m_e c^2 \)).

Returns:

The value of the differential cross-section (unnormalized).

__device__ void sample_electron_distr_p(double k[4], double p[4], double Thetae, curandState *localState)

Samples the 4-momentum of a scattering electron from a thermal distribution.

Implements the rejection-sampling scheme described by Canfield (1987) and Wong (2022). The routine picks an electron velocity and orientation relative to the incoming photon, weighted by the frame-dependent Klein-Nishina cross-section.

Routine:

1. Sample \( \gamma_e \) and \( \beta_e \) from a Maxwell-Jüttner distribution.

2. Sample the relative angle cosine \( \mu \) between the photon and electron.

3. Calculate the rest-frame photon energy \( K \).

4. Accept the candidate electron if a random variable is less than \( \sigma_{KN}(K)/\sigma_T \).

5. Construct a local orthonormal basis aligned with the photon’s direction to assign the electron’s spatial momentum components.

Parameters:

• k – Incoming photon 4-momentum in the local tetrad frame.

• p – [out] Sampled electron 4-momentum in the local tetrad frame.

• Thetae – Dimensionless electron temperature \( \Theta_e = k_B T_e / m_e c^2 \).

• localState – Pointer to the curand RNG state.

Returns:

void

__device__ void sample_beta_distr(double Thetae, double *gamma_e, double *beta_e, curandState *localState)

Samples the electron Lorentz factor and velocity from a thermal distribution.

This function determines the energy state of a scattering electron by sampling a transformed energy variable \( y \) and converting it into the Lorentz factor \( \gamma_e \) and the dimensionless velocity \( \beta_e \).

The sampling is based on the Maxwell-Jüttner distribution. The transformation used is:

\[ \gamma_e = y^2 \Theta_e + 1 \]

.

Parameters:

• Thetae – Dimensionless electron temperature \( \Theta_e = k_B T_e / m_e c^2 \).

• gamma_e – [out] Pointer to store the sampled Lorentz factor \( \gamma \).

• beta_e – [out] Pointer to store the sampled dimensionless velocity \( \beta = v/c \).

• localState – Pointer to the curand RNG state for the current thread.

Returns:

void

__device__ double sample_y_distr(const double Thetae, curandState *localState)

Samples the auxiliary energy variable \( y \) for the Maxwell-Jüttner distribution.

• This function implements the composition-rejection algorithm described by Canfield et al. (1987) to efficiently sample electron energies. It decomposes the distribution into a weighted sum of four components, each sampled from a Chi-Square distribution with varying degrees of freedom.

• Routine:

1. Weights: Calculates weights \( \pi_3, \pi_4, \pi_5, \pi_6 \) based on the dimensionless temperature \( \Theta_e \).

2. Composition: Randomly selects one of the four components based on these weights using a uniform random variable.

3. Sampling: Samples \( x \) from a \( \chi^2 \) distribution with 3, 4, 5, or 6 degrees of freedom.

4. Transformation: Converts \( x \) to the auxiliary variable \( y = \sqrt{x/2} \).

5. Rejection: Applies a final rejection step to correct the shape of the approximated distribution using the probability:

   \[ P = \frac{\sqrt{1 + \frac{1}{2} \Theta_e y^2}}{1 + y \sqrt{\frac{\Theta_e}{2}}} \]

Parameters:

• Thetae – Dimensionless electron temperature \( \Theta_e = k_B T_e / m_e c^2 \).

• localState – Pointer to the curand RNG state for the current thread.

Returns:

The sampled auxiliary variable \( y \), where \( \gamma_e = y^2 \Theta_e + 1 \).

__device__ double sample_mu_distr(const double beta_e, double random)

Samples the cosine of the angle between the photon and the electron. This function performs an inversion sampling of the relative orientation angle \( \mu = \cos\alpha \) between the incoming photon and the scattering electron. In the Monte Carlo scattering scheme, the probability of a collision is proportional to the relative velocity factor:

\[ P(\mu) \propto (1 - \beta_e \mu) \]

where \( \mu \in [-1, 1] \).

• Mathematical Derivation: The Cumulative Distribution Function (CDF) is obtained by integrating the probability:

  \[ R = \frac{\int_{-1}^{\mu} (1 - \beta_e x) dx}{\int_{-1}^{1} (1 - \beta_e x) dx} \]
  Solving for \( \mu \) via the quadratic formula yields the implementation:
  \[ \mu = \frac{1 - \sqrt{(1 + \beta_e)^2 - 4 \beta_e R}}{\beta_e} \]

Parameters:

• beta_e – Electron velocity in units of the speed of light ( \( v/c \)).

• random – A uniform random number in the range [0, 1].

Returns:

The sampled cosine of the relative angle \( \mu \).

__device__ void scatter_super_photon(struct of_photonSOA ph, struct of_photonSOA php, double Ne, double Thetae, double B, double Ucon[NDIM], double Bcon[NDIM], double Gcov[NDIM][NDIM], curandState *localState, unsigned long long photon_index)

Manages the Compton scattering process for a superphoton on the GPU. This function handles the full scattering pipeline for a single photon indexed in a Structure of Arrays (SOA). It transforms the photon into the local fluid frame, samples a scattering electron, executes the scattering event, and updates the photon’s coordinate-frame properties.

Routine:

• Tetrad Construction: Builds a local orthonormal tetrad using the fluid velocity \( u^\mu \) and magnetic field \( b^\mu \).

• Frame Transformation: Rotates the photon 4-momentum \( k^\mu \) from the coordinate basis to the local tetrad basis.

• Electron Sampling: Select a specific electron momentum from the Maxwell-Jüttner distribution.

• Scattering: Executes the scattering selecting the scattered photon momentum.

• Coordinate Recovery: Transforms the scattered momentum back to the global coordinate frame and updates the photon’s scattering history and energy.

Parameters:

• ph – Input photon Structure of Arrays (SOA).

• php – Output (scattered) photon Structure of Arrays (SOA).

• Ne – Local electron number density.

• Thetae – Dimensionless electron temperature \( \Theta_e \).

• B – Magnetic field strength in cgs.

• Ucon – Fluid 4-velocity \( u^\mu \).

• Bcon – Magnetic field 4-vector \( b^\mu \).

• Gcov – Covariant metric tensor \( g_{\mu\nu} \).

• localState – Pointer to the curand RNG state for this thread.

• photon_index – The index of the photon within the SOA.

Returns:

void