Polarized GRRT Code: Relativistic Radiative Modeling

• Polarized GRRT codes are computational frameworks that simulate electromagnetic propagation and polarization evolution in curved spacetimes.
• They employ high-accuracy numerical methods, such as adaptive Runge–Kutta schemes, to handle stiff Faraday effects and complex plasma interactions.
• Coupling with GRMHD simulations enables realistic imaging of polarized emissions from jets, accretion flows, and black hole environments.

Polarized General Relativistic Radiative Transfer (GRRT) codes implement the covariant propagation and polarization evolution of electromagnetic radiation through magnetized, relativistic plasmas in curved spacetimes. Such codes are essential for modeling the observable signatures of accretion flows, jets, and black hole environments where both strong gravity and polarized emission are present. The computational frameworks described in recent literature, including RaptorP, ipole, Blacklight, and others, provide robust, physically rigorous platforms for simulating polarized radiative transfer in stationary or time-dependent GRMHD backgrounds, recovering linear and circular polarization effects critical for interpreting modern VLBI and X-ray polarization data.

1. Fundamental Equations and Covariant Formulation

Covariant polarized GRRT is governed by a system of coupled ordinary differential equations along null geodesics $x^\mu(\lambda)$, parameterized by affine parameter $\lambda$. The formalism typically splits into two parts:

• Geodesic propagation:

$$\frac{dk^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha k^\beta = 0$$

$$\frac{df^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha f^\beta = 0$$

Here, $k^\mu$ is the photon four-momentum, $f^\mu$ is the polarization vector (with $k_\mu f^\mu = 0$), and $\Gamma^\mu_{\alpha\beta}$ is the Christoffel symbol for the metric $g_{\mu\nu}(x)$.

• Polarized transport of the Stokes vector $S = (I, Q, U, V)^T$ in the plasma frame:

$\lambda$0

with $\lambda$1 the vector of synchrotron (or cyclo-synchrotron) emissivities $\lambda$2, and $\lambda$3 the $\lambda$4 propagation matrix:

$\lambda$5

The coefficients $\lambda$6 represent absorption/extinction, while $\lambda$7 govern Faraday rotation and conversion.

The full covariant formulation often employs the Hermitian polarization tensor $\lambda$8 (the coherency tensor), which is evolved along the geodesic and projected into the comoving fluid tetrad at each step to recover standard Stokes parameters.

2. Numerical Algorithms and Stiffness Handling

GRRT codes implement direct integration schemes, most commonly Runge–Kutta (RK4 or RKF45), for the geodesic and polarization ODEs. Operator splitting strategies ("vacuum" and "plasma" steps) are standard; for regions where Faraday terms are stiff ($\lambda$9), codes diagonalize $$\frac{dk^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha k^\beta = 0$$0 and apply semi-analytic or matrix exponential solutions based on the Landi Degl'Innocenti formalism (see eqs. in ipole and RaptorP). This prevents numerical instability and preserves polarization evolution accuracy regardless of local plasma depth.

Adaptive step-size control is implemented to enforce error tolerances such that $$\frac{dk^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha k^\beta = 0$$1 for geodesic propagation and $$\frac{dk^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha k^\beta = 0$$2 per pixel for Stokes evolution. Many codes (RaptorP, Blacklight, BHOSS) provide fallback schemes: explicit integration where $$\frac{dk^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha k^\beta = 0$$3 are small, implicit or semi-analytic sub-cycling where either absorption or Faraday depth becomes large.

3. Coupling to MHD Backgrounds and Emission Physics

At each spatial location along a geodesic, the code interpolates primitive GRMHD variables (density $$\frac{dk^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha k^\beta = 0$$4, internal energy $$\frac{dk^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha k^\beta = 0$$5, flow four-velocity $$\frac{dk^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha k^\beta = 0$$6, magnetic field $$\frac{dk^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha k^\beta = 0$$7) from simulation snapshots (often from codes like PLUTO, HARM, or Athena++). The local emission and transfer coefficients $$\frac{dk^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha k^\beta = 0$$8, $$\frac{dk^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha k^\beta = 0$$9, $$\frac{df^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha f^\beta = 0$$0 are computed in the fluid frame using analytic formulas for the chosen electron distribution—thermal (Maxwell–Jüttner), power-law, or kappa distribution models (e.g. Melrose 1971, Mahadevan 1996, Dexter 2016).

Projection into orthonormal tetrads comoving with the plasma is performed at each step, as required for the correct calculation of local pitch angles $$\frac{df^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha f^\beta = 0$$1 and Doppler aberration. Emissivity is Doppler-boosted and transformed to the observer frequency via $$\frac{df^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha f^\beta = 0$$2.

4. Data Structures, Parallelization, and Performance

Typical GRRT codes maintain:

• Snapshot objects: GRMHD data and optional tracers (jet regions, temperature cuts).
• Ray objects: Arrays of $$\frac{df^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha f^\beta = 0$$3 for each sub-step; accumulators for $$\frac{df^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha f^\beta = 0$$4.
• Coefficient calculators: Return transfer coefficients given local plasma state and frequency.
• Main drivers: Distributed over MPI ranks by image-plane domain decomposition, with OpenMP threading over spatial rows/pixels.

Optimization frequently includes

• CUDA-porting of integration and interpolation hotspots (RaptorP): $$\frac{df^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha f^\beta = 0$$5 million geodesic steps/s and $$\frac{df^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha f^\beta = 0$$6 million transfer steps/s on a 4-GPU node, with $$\frac{df^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha f^\beta = 0$$7 accuracy.
• Overlapping I/O and computation: async snapshot loading while rays are processed.
• Dynamic load balancing to mitigate variable ray lengths.

Typical throughput enables $$\frac{df^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha f^\beta = 0$$8 polarimetric images in seconds to minutes on multi-core nodes; $$\frac{df^\mu}{d\lambda} + \Gamma^\mu_{\alpha\beta} k^\alpha f^\beta = 0$$9 images in $k^\mu$0 minutes across $k^\mu$1 cores.

5. Validation and Code-to-Code Benchmarking

Validation strategies uniformly encompass:

• Analytic solutions: Uniform/slab tests against Burn (1966), Landi Degl’Innocenti (1985), Chandrasekhar scattering matrices. Machine-precision agreement in constant-coefficient scenarios.
• Parallel transport and frame transformation: Conservation of the Walker–Penrose constant in Kerr spacetime; EVPA recovery in standard disk models.
• Inter-code benchmarks: Pixel-wise Stokes images compared among RaptorP, IPOLE, BHOSS, grtrans, and Odyssey. Normalized mean-squared errors (NMSE) for realistic GRMHD snapshots are generally $k^\mu$2 for $k^\mu$3 and $k^\mu$4, $k^\mu$5 for $k^\mu$6 (Prather et al., 2023).
• Error thresholds and best practices: $k^\mu$7 total-flux error, sub-degree systematic EVPA errors, and full code-interchangeability for EHT-class analyses.

6. Example Applications and Observational Signatures

Polarized GRRT codes have modeled:

• Multi-frequency polarized imaging of relativistic jets: RaptorP recovers frequency-dependent morphologies (ambient-dominated at low $k^\mu$8, spine/kink-dominated at high $k^\mu$9), and limb-brightening/EVPA transitions from poloidal to toroidal field signatures.
• QPO analysis: Light curves in $f^\mu$0, degree of polarization, and EVPA show quasi-periodic oscillations at rotation periods of kink structures (20–46 mHz), with spectral peaks matching simulated jet physics.
• Sgr A*, M87, and blazar polarization: Faraday rotation and conversion are key to reproducing observed LP, CP fractions, and rotation measures as a function of black hole spin and plasma temperature (Tsunetoe et al., 2020, Shcherbakov et al., 2010).
• Polarization loops and EVPA swings: clock-wise rotation of $f^\mu$1–$f^\mu$2 in light curves uniquely identifies underlying magnetic topology.

7. Code Availability and Implementation Considerations

Codes such as RaptorP (https://github.com/Lace-t/RaptorP), ipole (https://github.com/moscibrodzka/ipole), Blacklight (public domain), grtrans (https://github.com/jadexter/grtrans), Coport (https://github.com/JieweiHuang/Coport), Arcmancer, Pandurata, and GYOTO 2.0 are publicly available, doctrined for modular adaptation. Users can build executables with MPI+OpenMP+CUDA (where applicable), supply GRMHD snapshots, select electron distributions and frequencies, launch multi-frequency imaging, and generate FITS/PNG outputs for scientific postprocessing.

Installation and workflow typically involve:

• Compilation of core and GPU modules
• Configuration of observer/camera properties
• Selection of physical models (electron distribution, emission physics)
• Ray-launch through image plane
• Adaptive integration and storage of final emergent Stokes parameters

Genuine limitations include: fast-light approximations; handling of nonthermal distributions and scattering; need for special treatment near event horizons and in highly Faraday-thick media (where implicit solvers are advised). Planned extensions include support for multi-species plasmas, arbitrary metrics, advanced mesh refinement, and full Compton scattering kernels.

---

Polarized GRRT codes represent the state-of-the-art computational bridge between relativistic plasma dynamics and direct observational diagnostics, delivering polarimetric spectra and images with error quantified beneath current instrument uncertainties (Prather et al., 2023). Their modular architectures and validated algorithms form the foundation of modern high-energy, VLBI, and X-ray polarimetry modeling pipelines.