HWP Systematics#

The description and simulation of Half-Wave Plate (HWP) systematic effects in a realistic manner including explicitly the influence of the harmonics of the rotation frequency of the HWP. However, in simple cases, a description of the HWP systematics that does not take the harmonics explicitly into account can be applied. For this reason, litebird_sim includes two ways of simulating the systematics effects of a HWP.

The simpler case is already described in here, where a given hwp mueller matrix is directly used to compute the TOD with the scan_map() routine.

The explicit harmonics case is explained below. The high-level way of choosing between one or the other description is by setting the attribute harmonic_expansion to True or False when instantiating the HWP object. See more about this in Half Wave Plate.

HWP Harmonics Formalism#

This module implements rotation frequency harmonics explicit HWP non-idealities using Mueller’s formalism (see Jones formalism below) (as described in Patanchon et al. 2023). The Mueller matrix describing the effect of a HWP in the signal received by a detector can be expanded into the harmonics of the HWP rotation frequency f, from which we pick the 0th, 2nd and 4th terms.

(1)#\[M^{\text{HWP}}_{ij} = M^{0f}_{ij}(\Theta) + M^{2f}_{ij}(\Theta, \Phi, 2\rho) + M^{4f}_{ij}(\Theta, \Phi, 4\rho)\]

where:

  • \(\Theta\) is the incidence angle between the center of the HWP and the detector position;

  • \(\Phi\) is the angle describing the azimuthal position of the detector in the focal plane;

  • \(\rho\) is the HWP rotation angle in focal plane coordinates (please note that this is not the case in Patanchon et al. 2023.) .

The \(\Theta\) dependence is included in the input matrices, which are then multiplied at each sample by the trigonometric term depending harmonically on \(\rho\) and \(\Phi\), such that:

(2)#\[\begin{split}M^{0f}_{ij} &= M_{input_{ij}}^{0f} \\ M^{2f}_{ij} &= M_{input_{ij}}^{2f} \cos(2\rho - 2\Phi - \phi_{ij}^{2f}) \\ M^{4f}_{ij} &= M_{input_{ij}}^{4f} \cos(4\rho - 4\Phi - \phi_{ij}^{4f})\end{split}\]

where \(\phi_{ij}\) are harmonic and element dependent phases obtained trough EM simulations. Because the Mueller matrix for an ideal HWP (for a detector at azimuthal angle \(\Phi = 0\)) is given by:

(3)#\[\begin{split}M_{\text{IdealHWP}} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos{4\rho} & \sin{4\rho} \\ 0 & \sin{4\rho} & -\cos{4\rho} \end{pmatrix}\end{split}\]

which, in the HWP reference frame (\(\rho = 0\)) yields the usual matrix for an ideal HWP (diagonal 1, 1 -1 terms), we can now see that an ideal HWP is obtained by setting:

(4)#\[ \begin{align}\begin{aligned}\begin{split}M^{0f}_{input_{ideal}} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\end{split}\\\begin{split}M^{2f}_{input_{ideal}} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\end{split}\\\begin{split}M^{4f}_{input_{ideal}} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 1 & 1 \end{pmatrix}\end{split}\end{aligned}\end{align} \]

The expected sin and negative cosine terms are accounted for by the phases \(\phi\) in the expressions shown above, which is (in an ideal case) \(\pi\) for the 4f UU element, and to \(\frac{\pi}{2}\) for IU and UQ (sin terms).

Changing one of the values in the matrices above represents a non-ideality in the hwp.

Jones Formalism#

In the Jones formalism, a non-ideal static HWP can be described by perturbations \(\delta_{ij}\) to each ij element of the ideal HWP Jones Matrix, i.e.:

\[\begin{split}J = \begin{pmatrix} 1+\delta_{11} & \delta_{11} \\ \delta_{21} & -1 + \delta_{22}\end{pmatrix}\end{split}\]

These perturbations can be related to the Mueller matrix elements through the relation:

\[M_{ij} = \frac{1}{2} \mathrm{Tr}(\sigma_i J \sigma_j J^{\dagger})\]

that is, for each element of the matrix:

\[\begin{split}M_{II} &= Re(\delta_{11} - \delta_{22} + 1) \\ M_{QI} &= Re(\delta_{11} + \delta_{12}) \\ M_{UI} &= Re(\delta_{21} - \delta_{12}) \\ M_{IQ} &= Re(\delta_{11} + \delta_{12}) \\ M_{IU} &= Re(\delta_{12} - \delta_{21}) \\ M_{QQ} &= Re(\delta_{12} - \delta_{22} + 1) \\ M_{UU} &= Re(\delta_{22} + \delta_{11} - 1) \\ M_{UQ} &= Re(\delta_{12} + \delta_{21}) \\ M_{QU} &= Re(\delta_{12} + \delta_{21})\end{split}\]

We expand each of the perturbations $delta_{pq}$ into harmonics of the rotation frequency as:

\[\delta_{pq} = \sum_n \gamma_{pq}^{n_{f}} exp(in\alpha)\]

where \(p\), \(q\) are the matrix indexes, \(n_{f} = 0,2,4\) are the harmonics and \(\alpha\) is the HWP angle in focal plane coordinates (assuming for simplicity that the azimuthal angle of the detector is set to 0, and ignoring the dependence on the incidence angle). If our input matrices have the \(\gamma\) values as complex numbers \(A e^{i \phi}\), for each harmonic, we get:

\[\delta_{pq}^{n_{f}}( n_{f} \alpha) = \sum_{n_{f}} A_{pq}^{n_{f}} e^{i\phi} e^{in\alpha}\]

Developing this and taking only the real part we get to:

\[\delta_{pq}^{n_{f}}(n_{f} \alpha) = \sum_{n_{f}} A_{pq}^{n_{f}} \ \cos(n\alpha + n\phi_{ij}^{n_{f}})\]

In terms of implementation in LBS, this is done in 3 steps:

Then, the obtain matrix is rotated from the hwp frame to the focal plane frame, in order to apply the Mueller formalism explained above. In an ideal case, the jones perturbation matrices will be composed of only zero values and the transformation will yield the ideal hwp mueller matrix.

Band Integration#

Currently, the integration in band is available only for the Jones formalism. It works by the trapezoidal rule:

\[\sum_i \sum_f \big[d_f(i) + d_{f+1}(i)\big] \frac{\nu_{f+1} - \nu_{f}}{2}\]

where \(d_{f}(i)\) is the tod computed for a sample i, and at frequency at index f, and \(\nu_{f}\) if the frequency in Hz at index f.

The steps to perform band integration are:

1 - give the path to a csv/txt file containing the jones parameters for all frequencies of the instrument to the attribute jones_per_freq_csv_path when instantiating NonIdealHWP object.

2 - set integrate_in_band to True when calling scan_map() (or hwp_harmonics.fill_tod() if working at a lower-level).

The csv/txt file should have the following columns:

freq

Jxx_0f

Phxx_0f

Jxy_0f

Phxy_0f

Jyx_0f

Phyx_0f

Jyy_0f

Phyy_0f

Jxx_2f

Phxx_2f

Jxy_2f

Phxy_2f

Jyx_2f

Phyx_2f

Jyy_2f

Phyy_2f

70

0.93

40

0.001

100

0.003

32

0.95

170

0.012

85

0.0008

95

0.0025

28

0.011

92

71

0.92

42

0.0012

102

0.0032

34

0.94

168

0.013

87

0.0009

97

0.0027

30

0.012

94

72

0.94

38

0.0009

98

0.0028

30

0.96

172

0.011

83

0.0007

93

0.0023

26

0.010

90

Examples#

Examples can be seen in the notebooks:

1 - Performing a simulation with explicit HWP harmonics expansion 2 - Performing band integration with a Jones Matrix

HWP emissivity#

The module hwp_diff_emiss implements the additive contamination coming from HWP thermal emission as described in Micheli et al. 2024. This model complements the framework of HWP systematics since this effect cannot be introduced directly with Mueller’s formalism.

At finite temperature \(T_\textsc{hwp}\), HWPs emit blackbody radiation (in \(\mathrm{W}\,\mathrm{str}^{-1}\,\mathrm{m}^{-2}\,\mathrm{Hz}^{-1}\) units) weighted by their frequency-dependent emissivity, \(\varepsilon(\nu)\)

\[\varepsilon(\nu) I_\textsc{hwp}(\nu, T_\textsc{hwp}) = \varepsilon(\nu) \frac{2h\nu^3}{c^2} \frac{1}{e^{\frac{h\nu}{kT_\textsc{hwp}}}-1}\,,\]

If the optical properties differ between the ordinary and extraordinary axes, the HWP will also exhibit differential emissivity (Pisano et al. 2015).

In the HWP’s own reference frame, with the \(x\) and \(y\) axes aligned to the ordinary and extraordinary axes respectively, the emitted radiation can be represented as the Stokes vector

\[\begin{split}\mathbf{S}_\textsc{hwp}(\nu) = I_\textsc{hwp}(\nu, T_\textsc{hwp})\begin{pmatrix} \varepsilon(\nu) \\ \Delta\varepsilon(\nu) \\ 0 \end{pmatrix}\,,\end{split}\]

where \(\Delta\varepsilon(\nu)\) parametrizes the polarized component of the emission.

To simplify the model, we consider band-averaged quantities, which are calculated through a sensitivity calculation code. Thus, for each channel \(i\), the Stokes vector can be written as

\[\begin{split}\mathbf{S}_{\textsc{hwp},i} = I_{\textsc{hwp},i}(T_\textsc{hwp}) \begin{pmatrix} \varepsilon_i \\ \Delta\varepsilon_i \\ 0 \end{pmatrix}\,.\end{split}\]

To determine the corresponding contribution to the TOD, we rotate the Stokes vector from the HWP frame to the detector frame by an angle \((\xi-\rho)\), where \(\rho\) and \(\xi\) are the HWP and detector orientations in the focal-plane reference frame. This rotation is denoted by \(\mathcal{R}_{\xi-\rho}\), and the detected signal is then

\[d_{\textsc{hwp},i} = \begin{pmatrix} 1 & 1 & 0 \end{pmatrix} \mathcal{R}_{\xi-\rho} \mathbf{S}_{\textsc{hwp},i}\,.\]

More explicitly:

\[\begin{split}\begin{aligned} d_{\textsc{hwp},i} &= I_{\textsc{hwp},i}(T_\textsc{hwp}) \begin{pmatrix} 1 & 1 & 0 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos[2(\xi-\rho)] & \sin[2(\xi-\rho)] \\ 0 & -\sin[2(\xi-\rho)] & \cos[2(\xi-\rho)] \\ \end{pmatrix} \begin{pmatrix} \varepsilon_i \\ \Delta\varepsilon_i \\ 0 \end{pmatrix}\nonumber\\ &= I_{\textsc{hwp},i}(T_\textsc{hwp}) \begin{pmatrix} 1 & 1 & 0 \end{pmatrix} \begin{pmatrix} \varepsilon_i \\ \Delta\varepsilon_i \cos[2(\xi-\rho)] \\ -\Delta\varepsilon_i \sin[2(\xi-\rho)] \end{pmatrix}\nonumber\\ &= \varepsilon_i\,I_{\textsc{hwp},i}(T_\textsc{hwp}) + \Delta\varepsilon_i\, I_{\textsc{hwp},i}(T_\textsc{hwp}) \cos[2(\xi-\rho)]\,. \end{aligned}\end{split}\]

The first term represents the unpolarized thermal emission, while the second term is the polarized emission, modulated at twice the relative angle between the HWP and detector.

To simulate the 2f signal from a rotating, emitting HWP, one can use the method of Simulation class Simulation.add_2f(), or any of the low-level functions: add_2f_to_observations(), add_2f_for_one_detector().

Examples#

The examples below skip the simulation and observation creation for brevity. If needed, the implementation for those parts is explained in other sections of the docs.

(... importing modules, creating simulation, setting scanning strategy, instrument, etc...)

sim.set_hwp(lbs.IdealHWP(hwp_radpsec))

# creating the detectors (two mock detectors here)
dets = [
     lbs.DetectorInfo(name="det_A", sampling_rate_hz=sampling_hz),
     lbs.DetectorInfo(name="det_B", sampling_rate_hz=sampling_hz),
 ]

(obs,) = sim.create_observations(
   detectors=dets,
)

sim.prepare_pointings()

# Define differential emission amplitude for the detectors, in the same temperature units you used for the TOD.
# If just one value is set, all the detectors will have the same 2f amplitude.
sim.observations[0].amplitude_2f_k = np.array([0.1, 0.2])

# Adding 2f signal from HWP differential emission using the `Simulation` class method
# By default, the TOD is added to ``Observation.tod``. If you want to add it to some
# other field of the :class:`.Observation` class, use `component`

sim.add_2f()

API reference#

litebird_sim.hwp_harmonics.hwp_harmonics.set_band_params_for_one_detector(hwp: NonIdealHWP, band_params: HWPJonesParams, bandcenter_ghz: float, bandwidth_ghz: float, bandpass: dict[str, object] | None, include_beam_throughput: bool) → tuple[HWPJonesParams, ndarray[tuple[Any, ...], dtype[float64]]]#

Load frequency-dependent parameters and normalized bandpass profiles for a single detector based on HWP Jones coefficients.

This function retrieves the necessary HWP coefficients, filters them by the detector’s frequency range, and calculates the integrated bandpass profile (bpi) weighted by the derivative of the Planck function.

Parameters:

  • hwp – The HWP object

  • bandcenter_ghz – Central frequency of the detector band.

  • bandwidth_ghz – Width of the detector band.

  • bandpass – Dictionary containing bandpass profile metadata.

  • include_beam_throughput – Whether to include beam effects in the profile.

Returns:

  • An instance of HWPJonesParams.

  • A numpy array representing the normalized bandpass profile (bpi).

Return type:

A tuple containing

Raises:

NotImplementedError – If the HWP calculus is set to Mueller instead of Jones.

litebird_sim.hwp_harmonics.hwp_harmonics.fill_tod_with_hwp_harmonics(hwp: NonIdealHWP, observation: Observation, tod: ndarray, maps: HealpixMap | dict[str, ~litebird_sim.maps_and_harmonics.HealpixMap | ~litebird_sim.input_sky.SkyGenerationParams] | ~litebird_sim.maps_and_harmonics.SphericalHarmonics | dict[str, ~litebird_sim.maps_and_harmonics.SphericalHarmonics | ~litebird_sim.input_sky.SkyGenerationParams] | None=None, pointings: ndarray | None = None, hwp_angle: ndarray | None = None, save_tod: bool = True, input_names: list[str] | None = None, pointings_dtype=<class 'numpy.float64'>, apply_non_linearity: bool = False, add_2f_hwpss: bool = False, mueller_phases: dict[str, ~numpy.ndarray] | None=None, integrate_in_band: bool = False, nthreads: int | None = None, include_beam_throughput: bool = False)#

Fill a TOD for one observation, using HWP rotation speed harmonics calculus.

Parameters:

  • hwp – HWP, optional Half-wave plate (HWP) model. If None, HWP effects are ignored unless the Observation object contains HWP data.

  • observation (Observation) – container for the TOD,detectors and pointings. If the TOD is not required, you can avoid allocating observations.tod by setting allocate_tod=False in Observation.

  • tod – np.ndarray Time-ordered data (TOD) array of shape (n_detectors, n_samples) that will be filled with the simulated sky signal.

  • maps – HealpixMap | SphericalHarmonics | dict[str, HealpixMap] | dict[str, SphericalHarmonics] Sky model to be scanned. In dictionary form, the keys must match the entries of input_names.

  • pointings (optional) – if not present, it is either computed on the fly (generated by lbs.get_pointings() per detector), or read from observations.pointing_matrix (if present). pointings must be a np.array of shape (N_det, N_samples, 3).

  • hwp_angle (optional) – 2ωt, hwp rotation angles (radians). If pointings is passed, hwp_angle must be passed as well, otherwise both must be None. If not passed, it is computed on the fly (generated by lbs.get_pointings() per detector). hwp_angle must be a np.array of dimensions (N_samples).

  • save_tod (bool) – if True, obs.tod is saved in observations.tod and locally as a .npy file; if False, obs.tod gets deleted.

  • input_names – array-like of str or None, default=None Per-detector keys to select entries in maps when maps is a dictionary.

  • pointings_dtype – if pointings is None and is computed within fill_tod, this is the dtype for pointings and tod (default: np.float32).

  • apply_non_linearity (bool) – applies the coupling of the non-linearity systematics with hwp_sys

  • add_2f_hwpss (bool) – adds the 2f hwpss signal to the TOD

  • mueller_phases (dict) – the non ideal phases for the mueller matrix elements. When None is given, the temporary values from Patanchon et al. [2021] are used.

  • integrate_in_band – bool, default=False Whether to integrate the signal over the detector’s frequency band. Only implemented for the Jones formalism.

  • nthreads (int) – number of threads to use in ducc’s Healpix methods. If None, the function reads from the OMP_NUM_THREADS environment variable.

  • include_beam_throughput – bool, default=False Whether to include beam throughput in the bandpass in the bandpass profile.

Raises:

NotImplementedError – If integrate_in_band is True and the HWP calculus is set to Mueller.

litebird_sim.hwp_harmonics.mueller_methods.compute_signal_for_one_detector(tod_det, m0f, m2f, m4f, rho, psi, mapT, mapQ, mapU, cos2Xi2Phi, sin2Xi2Phi, phi, xi, apply_non_linearity, g_one_over_k, add_2f_hwpss, amplitude_2f_k, phases_2f, phases_4f)#

Single-frequency case: compute the signal for a single detector, looping over (time) samples.

litebird_sim.hwp_harmonics.jones_methods.compute_signal_for_one_detector(tod_det, deltas_j0f, deltas_j2f, rho, psi, mapT, mapQ, mapU, cos2Xi2Phi, sin2Xi2Phi, phi, xi, apply_non_linearity, g_one_over_k, add_2f_hwpss, amplitude_2f_k)#

Single-frequency case: compute the signal for a single detector, looping over (time) samples.

litebird_sim.hwp_harmonics.jones_methods.integrate_inband_signal_for_one_detector(tod_det, freqs, band, deltas_j0f, deltas_j2f, mapT, mapQ, mapU, rho, psi, phi, xi, cos2Xi2Phi, sin2Xi2Phi, apply_non_linearity, g_one_over_k, add_2f_hwpss, amplitude_2f_k)#

Multi-frequency case: band integration of the signal for a single detector, looping over (time) samples.

litebird_sim.hwp_harmonics.common.integrate_inband_signal_for_one_sample(T, Q, U, freqs, band, mII, mQI, mUI, mIQ, mIU, mQQ, mUU, mUQ, mQU, cos2Xi2Phi, sin2Xi2Phi, cos2Psi2Phi, sin2Psi2Phi)#

Multi-frequency case: band integration with trapezoidal rule, \(\sum (f(i) + f(i+1)) \cdot (\nu_(i+1) - \nu_i)/2\) for a single (time) sample.

litebird_sim.hwp_diff_emiss.add_2f(tod, hwp_angle, pol_angle_rad: ndarray, amplitude_2f_k: float | ndarray)#

Add the HWP differential emission to some time-ordered data

This functions modifies the values in tod by adding the contribution of the HWP synchronous signal coming from differential emission. The amplitude_2f_k argument must be a N_dets array containing the amplitude of the HWPSS.

litebird_sim.hwp_diff_emiss.add_2f_to_observations(observations: ~litebird_sim.observations.Observation | list[~litebird_sim.observations.Observation], hwp: ~litebird_sim.hwp.HWP | None = None, component: str = 'tod', amplitude_2f_k: float | None = None, pointings_dtype=<class 'numpy.float64'>)#

Add the HWP differential emission to some time-ordered data

This is a wrapper around the add_2f() function that applies to the TOD stored in observations, which can either be one Observation instance or a list of observations.

By default, the TOD is added to Observation.tod. If you want to add it to some other field of the Observation class, use component:

for cur_obs in sim.observations:
    # Allocate a new TOD for the 2f alone
    cur_obs.hwp_2f_tod = np.zeros_like(cur_obs.tod)

    # Ask `add_2f_to_observations` to store the 2f
    # in `observations.hwp_2f_tod`
    add_2f_to_observations(sim.observations, component="hwp_2f_tod")