TESSA — physics and algorithms behind stellar magnetic imaging

TESSA (Tomographic Exploration of Stellar Surface Analysis Interface) is my implementation of Zeeman–Doppler imaging (ZDI): the reconstruction of a rotating star’s large-scale surface magnetic field from a time series of polarized line profiles.

This note describes the physical assumptions, mathematical representation, and numerical algorithm behind the code. The central problem is an inverse one: many different surface maps can reproduce a finite set of disk-integrated spectra. TESSA therefore combines an explicit forward model with maximum-entropy regularization to select the simplest map supported by the observations.

The information carried by Stokes V

In a magnetic field, a spectral line is split by the Zeeman effect. For an effective Landé factor $g_{\rm eff}$ and central wavelength $\lambda_0$, the characteristic wavelength displacement is

\[\Delta\lambda_Z = 4.6686\times10^{-12}\,g_{\rm eff}\,\lambda_0^2 B,\]

where $\lambda_0$ is in nm and $B$ is in gauss. In the weak-field regime the splitting is smaller than the intrinsic line width, but the oppositely polarized Zeeman components produce an antisymmetric circular-polarization signature. The local Stokes V profile is approximately

\[V_{\rm loc}(\lambda) =-\Delta\lambda_Z B_{\rm los} \frac{\partial I_{\rm loc}}{\partial\lambda},\]

so its amplitude and sign trace the line-of-sight component $B_{\rm los}$. TESSA normally works with Least Squares Deconvolution (LSD) profiles, in which hundreds or thousands of photospheric lines are combined into high signal-to-noise mean Stokes I and V profiles.

The spatial information comes from rotation. A surface element approaching the observer contributes on the blue side of the line, while a receding element contributes on the red side. As the star rotates, a magnetic region leaves a phase-dependent track across the velocity-resolved Stokes V profile. ZDI uses these tracks to infer longitude and, more weakly, latitude.

Small-scale mixed-polarity fields cancel in disk-integrated Stokes V. The recovered map should therefore be interpreted as the large-scale magnetic topology visible at the spectral resolution and rotational velocity of the observations.

Surface geometry and rotational encoding

TESSA divides the stellar surface into approximately equal-area latitude rings. For $N_{\rm rings}$ rings, the number of longitude cells at colatitude $\theta_k$ is chosen as

\[N_\phi^{(k)} \simeq 2N_{\rm rings}\sin\theta_k, \qquad \theta_k=\frac{\pi(k+1/2)}{N_{\rm rings}}.\]

This produces roughly $2N_{\rm rings}^2$ surface elements without strongly oversampling the polar regions. At every observation time, the code computes the visibility $\mu=\cos\theta_{\rm view}$, projected area, line-of-sight magnetic field, and projected rotational velocity of every cell.

The local intensity is weighted by a linear limb-darkening law,

\[L(\mu)=1-\varepsilon(1-\mu),\]

and rapid rotators can additionally include a Roche surface and gravity darkening. Solar-like differential rotation is represented by

\[\Omega(\theta)=\Omega_{\rm eq}-d\Omega\cos^2\theta,\]

so each latitude accumulates a slightly different phase over a multi-epoch data set. Geometry for all phases and cells is evaluated as a set of vectorized matrices and reused throughout the inversion.

Local line formation and disk integration

The default weak-field model represents local Stokes I with a Voigt profile. Its Gaussian width describes thermal and turbulent broadening, while its Lorentzian component represents the extended wings. Each surface element is Doppler-shifted by its projected rotation velocity and weighted by brightness, limb darkening, and projected area.

For phase $j$, the schematic disk-integrated profiles are

\[I_j(v)=\frac{\sum_{k\in\mathrm{visible}} w_{jk}\,I_{\mathrm{loc},k}(v-v_{jk})} {\sum_{k\in\mathrm{visible}}w_{jk}},\] \[V_j(v)=\frac{\sum_{k\in\mathrm{visible}} w_{jk}\,V_{\mathrm{loc},k}(v-v_{jk},B_{\mathrm{los},jk})} {\sum_{k\in\mathrm{visible}}w_{jk}},\]

where $w_{jk}$ contains the projected area, limb-darkening factor, gravity-darkening factor, and local brightness. Instrumental broadening is applied to the resulting profile.

For fields or lines outside the weak-field approximation, TESSA also provides an Unno–Rachkovsky model. This uses the analytic polarized radiative-transfer solution for a Milne–Eddington atmosphere, allowing Stokes I and V to respond nonlinearly to the magnetic field. Additional Hα models support emission profiles, self-absorption, and numerical templates for chromospheric applications.

Representing a vector magnetic field

A free three-component magnetic vector in every surface pixel would introduce too many poorly constrained parameters. TESSA instead expands the surface field in spherical harmonics and separates it into poloidal and toroidal components:

\[B_r(\theta,\phi)= \sum_{\ell=1}^{\ell_{\max}}\sum_{m=0}^{\ell} \operatorname{Re}\!\left[\alpha_{\ell m}Y_\ell^m(\theta,\phi)\right],\] \[B_\theta=-\sum_{\ell,m}\operatorname{Re}\!\left[ \beta_{\ell m}Z_{\ell m}+\gamma_{\ell m}X_{\ell m}\right], \qquad B_\phi=-\sum_{\ell,m}\operatorname{Re}\!\left[ \beta_{\ell m}X_{\ell m}-\gamma_{\ell m}Z_{\ell m}\right].\]

Here $Y_\ell^m$ is the normalized spherical harmonic; $X_{\ell m}$ and $Z_{\ell m}$ contain its azimuthal and meridional derivatives. The coefficients $\alpha_{\ell m}$ describe the radial poloidal field, $\beta_{\ell m}$ the horizontal poloidal field, and $\gamma_{\ell m}$ the toroidal field.

The truncation degree $\ell_{\max}$ sets the smallest permitted angular scale. Typical ZDI data constrain only low-order modes, so values around 10–15 are normally sufficient. TESSA can fit a full vector field or impose poloidal, poloidal-plus-toroidal, or potential-field constraints.

Associated Legendre functions and their derivatives are precomputed on the surface grid. Field evaluation then becomes matrix multiplication, and the same basis matrices provide analytic derivatives of the field with respect to every harmonic coefficient.

The inverse problem and maximum entropy

Let $\mathbf d$ be the observed profile samples, $\mathbf F(\mathbf x)$ the synthetic profiles, and $\mathbf x$ the brightness pixels and magnetic harmonic coefficients. TESSA measures the data misfit with

\[\chi_r^2(\mathbf x)=\frac{1}{N_d} \sum_{i=1}^{N_d} \frac{[d_i-F_i(\mathbf x)]^2}{\sigma_i^2}.\]

Because the inversion is underdetermined, minimizing $\chi^2$ alone would fit noise and create unsupported surface structure. TESSA uses the Skilling–Bryan maximum-entropy method: among the maps that reach a target $\chi_r^2$—normally close to one—it selects the map with the largest entropy,

\[\max_{\mathbf x} S(\mathbf x) \quad\text{subject to}\quad \chi_r^2(\mathbf x)=\chi_{\rm aim}^2.\]

Brightness pixels are positive and use a standard relative entropy around a default image. Magnetic harmonic coefficients can be positive or negative, so they use a signed entropy. For a coefficient $a_k$ and default scale $q$, this takes the form

\[S_{\rm mag}=\sum_k w_k\left[ \psi_k-2q-a_k\ln\left(\frac{\psi_k+a_k}{2q}\right) \right], \qquad \psi_k=\sqrt{a_k^2+4q^2}.\]

The default scale and harmonic weights penalize unnecessary magnetic complexity, especially poorly constrained high-order modes.

At each iteration TESSA computes the model, residuals, entropy gradient, and response matrix

\[R_{ik}=\frac{\partial F_i}{\partial x_k}.\]

The Skilling–Bryan algorithm constructs a small search subspace from the entropy and $\chi^2$ gradients and their curvature-weighted directions. It then solves for a bounded step that approaches the target $\chi^2$ while increasing entropy. Iteration stops when the requested fit is reached and the entropy and constraint gradients become locally parallel.

In practical terms, one inversion cycle is:

evaluate visibility and Doppler shifts at all observed phases;
synthesize disk-integrated Stokes I and V profiles;
calculate residuals and the profile Jacobian;
update brightness pixels and/or magnetic harmonic coefficients in the MEM search subspace;
repeat until the target fit and convergence criteria are satisfied.

Products and interpretation

TESSA returns fitted line profiles and maps of $B_r$, $B_\theta$, and $B_\phi$. From the harmonic coefficients it also derives magnetic-energy fractions in poloidal, toroidal, axisymmetric, dipolar, quadrupolar, and higher-order components. Optional potential-field source-surface extrapolation provides a visualization of the large-scale field above the reconstructed surface.

The result remains conditional on the assumed inclination, $v\sin i$, rotation period, differential rotation, line model, and regularization. Incomplete phase coverage creates longitude-dependent sensitivity; low inclination introduces north–south ambiguity; and maximum entropy suppresses weak small-scale features. These are properties of the inverse problem rather than merely numerical limitations.

Implementation and status

The physical model, inversion engine, and plotting layer are written in Python. A FastAPI backend exposes configuration, observations, run control, and results to a browser interface. A structured config.json records the stellar parameters and model choices, while command-line execution supports reproducible batch analyses.

TESSA is under active development. The repository is currently private while the interface, chromospheric models, and validation suite are being refined.

Method references

Donati, J.-F. & Brown, S. F. (1997), Zeeman–Doppler imaging of active stars. V.
Skilling, J. & Bryan, R. K. (1984), Maximum entropy image reconstruction: general algorithm.
Unruh, Y. C. & Collier Cameron, A. (1995), The sensitivity of Zeeman Doppler imaging to field orientation.
Hobson, M. P. & Lasenby, A. N. (1998), The entropic prior for distributions with positive and negative values.
Folsom, C. P. et al. (2018), The large-scale magnetic field and differential rotation of ξ Bootis A.