HWP_sys

This module implements HWP non-idealities both using Jones formalism (as described in Giardiello et al. 2021) and the Mueller one. In the Jones formalism, a non-ideal HWP is described by

(1)\[\begin{split}J_{HWP} = \begin{pmatrix} 1+h_{1} & \zeta_{1} e^{i \chi_1}\\ \zeta_{2} e^{i \chi_2}& -(1+h_{2}) e^{i \beta} \\ \end{pmatrix}\end{split}\]

where:

  • \(h_1\) and \(h_2\) are the efficiencies, describing the deviation from the unitary transmission of light components \(E_x\), \(E_y\). In the ideal case, \(h_1 = h_2 = 0\);

  • \(\beta=\phi-\pi\), where \(\phi\) is the phase shift between the two directions. It accounts for variations of the phase difference between \(E_x\) and \(E_y\) with respect to the nominal value of \(\pi\) for an ideal HWP. In the ideal case, \(\beta=0\);

  • \(\zeta_{1,2}\) and \(\chi_{1,2}\) are amplitudes and phases of the off-diagonal terms, coupling \(E_x\) and \(E_y\). In practice, if the incoming wave is fully polarized along x(y), a spurious y(x) component would show up in the outgoing wave. In the ideal case, \(\zeta_{1,2}=\chi_{1,2}=0\).

In the Mueller formalism, we have a general matrix

(2)\[\begin{split}J_{HWP} = \begin{pmatrix} M^{TT} & M^{TQ} & M^{TU} \\ M^{QT} & M^{QQ} & M^{QU} \\ M^{UT} & M^{UQ} & M^{UU} \\ \end{pmatrix}\end{split}\]

which, in the ideal case, would be

(3)\[\begin{split}J_{HWP} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -1 \\ \end{pmatrix}\end{split}\]

In the most general case, the Jones non-ideal parameters and the Mueller matrix elements can vary inside a bandpass.

The class hwp_sys.HwpSys is a contanier for the parameters of the HWP systematics. It defines three methods:

  • hwp_sys.HwpSys.set_parameters() which sets the defaults and handles the interface with the parameter file of the simulation. The relevant section is tagged by [hwp_sys]. The HWP parameters can be passed both in the Jones and Mueller formalism. This choice is regulated by the flag mueller_or_jones, which has to be set to either "mueller" or "jones". In case the Jones formalism is chosen, it is converted automatically into the Mueller one through the function .hwp_sys.JonesToMueller. There is also the possibility to pass precomputed input maps (as a numpy array) through the maps argument. Otherwise, input maps are internally computed through the litebird_sim map based simulator (mbs). The argument integrate_in_band sets whether band integration is performed in the TOD computation; if built_map_on_the_fly = True, the map-making can be performed internally (instead of using the litebird_sim binner); correct_in_solver sets whether non-ideal parameters can be used in the map-making (map-making assuming a non-ideal HWP, generally using different HWP non-ideal parameters than the one used in the TOD, representing our estimate of their true value); integrate_in_band_solver regulates whether band integration is performed in the map-making (to compute the \(B^T B\) and \(B^T d\) terms, see below).

  • hwp_sys.HwpSys.fill_tod() which fills the tod in a given Observation. The pointings angles passed have to include no rotating HWP, since the effect of the rotating HWP to the polarization angle is included in the TOD computation. The TOD is computed performing this operation:

    \(d_{obs}\left(t_{i}\right)\,=\,\frac{\int d\nu\,\frac{\partial BB(\nu,T)}{\partial T_{CMB}}\,\tau\left(\nu\right)\,M_{i}^{TX}(\nu)\left(m_{CMB}+m_{FG}\left(\nu\right)\right)}{\int d\nu \frac{\partial BB(\nu,T)}{\partial T_{CMB}}\,\tau \left(\nu\right)}\),

    where \(\tau(\nu)\) is the bandpass, \(\frac{\partial BB(\nu,T)}{\partial T_{CMB}}\) converts from CMB thermodynamic temperature to differential source intensity (see eq.8 of https://arxiv.org/abs/1303.5070) and \(M_{i}^{TX}(\nu)\) is the Mueller matrix element including the non-ideal HWP.

    If built_map_on_the_fly = True, the code computes also

    \(m_{out} = {\,\left(\sum_{i} B_{i}^{T} B_{i} \right)^{-1} \left( \sum_{i} B_{i}^{T} d_{obs}(t_{i}) \right)}\),

    where the map-making matrix \(B^X = \left(\frac{\int d\nu \,\frac{\partial BB(\nu,T)}{\partial T_{CMB}}\,\tau_{s}\left(\nu\right)\,M_{i,s}^{TX}(\nu)}{\int d\nu \frac{\partial BB(\nu,T)}{\partial T_{CMB}}\,\tau_{s}\left(\nu\right)}\,\right)\).

    \(\tau_s(\nu)\) and \(M_{i,s}^{TX}(\nu)\) are the estimate of the bandpass and Mueller matrix elements used in the map-making.

  • hwp_sys.HwpSys.make_map() which can bin the observations in a map. This is available only if built_map_on_the_fly variable is set to True. With this method, it is possible to include non-ideal HWP knowledge in the map-making procedure, so use that instead of the general litebird_sim binner if you want to do so.

Defining a bandpass profile in hwp_sys

It is possible to define more complex bandpass profiles than a top-hat when using hwp_sys. This can be done both for the TOD computation (\(\tau\)) and the map-making procedure (\(\tau_s\)). All you have to do is create a dictionary with key “hwp_sys” in the parameter file (a toml file) assigned to the simulation:

sim = lbs.Simulation(
    parameter_file=toml_filename,
    random_seed=0,
 )

The dictionary under the key “hwp_sys” will also contain the paths to the files from which the HWP parameters are read in the multifrequency case (under the keys “band_filename/band_filename_solver”), or their values in the single frequency one. See the notebook hwp_sys/examples/simple_scan for more details. To define the bandpasses to use, you need to have a dictionary with key “bandpass” (for \(\tau\)) or “bandpass_solver” (for \(\tau_s\)) under “hwp_sys”:

paramdict = {...
              "hwp_sys": {...
                 "band_filename": path_to_HWP_param_file,
                 "band_filename_solver": path_to_HWP_solver_param_file,
                 "bandpass": {"band_type": "cheby",
                              "band_low_edge": band_low_edge,
                              "band_high_edge": band_high_edge,
                              "bandcenter_ghz": bandcenter_ghz,
                              "band_ripple_dB": ripple_dB_tod,
                              "band_order": args.order_tod},
                 "bandpass_solver": {...},
               ...}}

The above example is for a bandpass with Chebyshev filter, but there are other parameters to define different bandpass profile. It is important to define the “band_type”, which can be “top-hat”, “top-hat-exp”, “top-hat-cosine” and “cheby” (see the bandpass module for more details) and the band edges, which define the frequency range over which the bandpass transmission is close or equal to 1. If not assigned, the “band_type” is automatically set to “top-hat” and the band edges will correspond to the limits of the frequency array used (which, in the hwp_sys module, is read from the HWP parameter files). There are default values also for the parameters defining the specific bandpass profiles (see the hwp_sys/hwp_sys/bandpass_template_module code).

There is also the possibility to read the bandpass profile from an external file, which has to be a .txt file with two columns, the frequency and the bandpass transmission. It is important that the frequency array used for “bandpass/bandpass_solver” coincides with the ones passed in the “band_filename/band_filename_solver” file. Here is how to pass the bandpass file:

paramdict = {...
              "hwp_sys": {...
                 "band_filename": path_to_HWP_param_file,
                 "band_filename_solver": path_to_HWP_solver_param_file,
                 "bandpass": {"bandpass_file": path_to_bandpass_file },
                 "bandpass_solver": {"bandpass_file": path_to_bandpass_solver_file},
               ...}}

You can find more examples for the bandpass construction in the hwp_sys/examples/simple_scan notebook.

API reference

HWP_sys

litebird_sim.hwp_sys.hwp_sys.JonesToMueller(jones)

Converts a Jones matrix to a Mueller matrix. The Jones matrix is assumed to be a 2x2 complex matrix (np.array). The Mueller matrix is a 4x4 real matrix. Credits to Yuki Sakurai for the function.

litebird_sim.hwp_sys.hwp_sys.get_mueller_from_jones(h1, h2, z1, z2, beta)

Converts the (frequency-dependent) input Jones matrix to a Mueller matrix. Returns Mueller matrix (3x3xNfreq), V-mode related terms are discarded, given the assumption of vanishing circular polarization.

Inputs: \(h_1\), \(h_2\), \(\zeta_1\), \(\zeta_2\), \(\beta\) (i.e. systematics of the HWP, not the full Jones matrix) Returns: \(M^{II}\), \(M^{QI}\), \(M^{UI}\), \(M^{IQ}\), \(M^{IU}\), \(M^{QQ}\), \(M^{UU}\), \(M^{UQ}\), \(M^{QU}\) (single/multi-frequency Mueller matrix terms)

litebird_sim.hwp_sys.hwp_sys.compute_signal_for_one_sample(T, Q, U, mII, mQI, mUI, mIQ, mIU, mQQ, mUU, mUQ, mQU, c2ThPs, s2ThPs, c2PsXi, s2PsXi, c2ThXi, s2ThXi, c4Th, s4Th)

Bolometric equation, tod filling for a single (time) sample

litebird_sim.hwp_sys.hwp_sys.compute_signal_for_one_detector(tod_det, mII, mQI, mUI, mIQ, mIU, mQQ, mUU, mUQ, mQU, pixel_ind, theta, psi, xi, maps)

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

litebird_sim.hwp_sys.hwp_sys.integrate_inband_signal_for_one_sample(T, Q, U, freqs, band, mII, mQI, mUI, mIQ, mIU, mQQ, mUU, mUQ, mQU, c2ThPs, s2ThPs, c2PsXi, s2PsXi, c2ThXi, s2ThXi, c4Th, s4Th)

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_sys.hwp_sys.integrate_inband_signal_for_one_detector(tod_det, freqs, band, mII, mQI, mUI, mIQ, mIU, mQQ, mUU, mUQ, mQU, pixel_ind, theta, psi, xi, maps)

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

litebird_sim.hwp_sys.hwp_sys.compute_TQUsolver_for_one_sample(mIIs, mQIs, mUIs, mIQs, mIUs, mQQs, mUUs, mUQs, mQUs, c2ThPs, s2ThPs, c2PsXi, s2PsXi, c2ThXi, s2ThXi, c4Th, s4Th)

Single-frequency case: computes \(A^T A\) and \(A^T d\) for a single detector, for one (time) sample.

litebird_sim.hwp_sys.hwp_sys.compute_atd_ata_for_one_detector(atd, ata, tod, mIIs, mQIs, mUIs, mIQs, mIUs, mQQs, mUUs, mUQs, mQUs, pixel_ind, theta, psi, xi)

Single-frequency case: compute \(A^T A\) and \(A^T d\) for a single detector, looping over (time) samples.

litebird_sim.hwp_sys.hwp_sys.integrate_inband_TQUsolver_for_one_sample(freqs, band, mIIs, mQIs, mUIs, mIQs, mIUs, mQQs, mUUs, mUQs, mQUs, c2ThPs, s2ThPs, c2PsXi, s2PsXi, c2ThXi, s2ThXi, c4Th, s4Th)

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_sys.hwp_sys.integrate_inband_atd_ata_for_one_detector(atd, ata, tod, freqs, band, mIIs, mQIs, mUIs, mIQs, mIUs, mQQs, mUUs, mUQs, mQUs, pixel_ind, theta, psi, xi)

Multi-frequency case: band integration of \(A^T A\) and \(A^T d\) for a single detector, looping over (time) samples.

class litebird_sim.hwp_sys.hwp_sys.HwpSys(simulation)

Bases: object

A container object for handling tod filling in presence of hwp non-idealities following the approach of Giardiello et al. 2021 https://arxiv.org/abs/2106.08031 :param simulation: an instance of the class Simulation :type simulation: Simulation

set_parameters(nside: int | None = None, mueller_or_jones: str = 'mueller', Mbsparams: MbsParameters | None = None, integrate_in_band: bool | None = None, built_map_on_the_fly: bool | None = None, correct_in_solver: bool | None = None, integrate_in_band_solver: bool | None = None, Channel: FreqChannelInfo | None = None, maps: None | ndarray = None, parallel: bool | None = None)

It sets the input paramters reading a dictionary sim.parameters with key “hwp_sys” and the following input arguments

Parameters:
  • nside (integer) – nside used in the analysis

  • mueller_or_jones (str) – “mueller” or “jones” (case insensitive) it is the kind of HWP matrix to be injected as a starting point if ‘jones’ is chosen, the parameters \(h_1\), \(h_2\), \(\beta\), \(\zeta_1\), \(\zeta_2\) are used to build the Jones matrix \(\begin{pmatrix} 1 + h_1 & \zeta_1 \\ \zeta_2 & - (1 + h_2) e^{i \beta} \\ \end{pmatrix}\) and then converted to Mueller. \(\zeta_1\), \(\zeta_2\) are assumed to be complex \(h_1\), \(h_2\), \(\beta\) are assumed to be real \(\beta\) is assumed to be in degrees (later converted to radians. To reproduce the ideal HWP case, set all Jones parameters to 0. If Mueller parameters are used, set \(M^{II/QQ} = 1\), \(M^{UU} = -1\) and all the others to 0.

  • Mbsparams (Mbs) – an instance of the Mbs class Input maps needs to be in galactic (mbs default)

  • integrate_in_band (bool) – performs the band integration for tod generation

  • built_map_on_the_fly (bool) – fills \(A^T A\) and \(A^T d\)

  • correct_in_solver (bool) – if the map is computed on the fly, fills \(A^T A\) using map-making (solver) HWP parameters

  • integrate_in_band_solver (bool) – performs the band integration for the map-making solver

  • Channel (FreqChannelInfo) – an instance of the FreqChannelInfo class

  • maps (float) – input maps (3, npix) coherent with nside provided, Input maps needs to be in galactic (mbs default) if maps is not None, Mbsparams is ignored (i.e. input maps are not generated)

  • parallel (bool) – uses parallelization if set to True

fill_tod(obs: Observation | List[Observation], hwp_radpsec: float, pointings: ndarray | List[ndarray] | None = None, save_tod: bool = False)

It fills tod and/or \(A^T A\) and \(A^T d\) for the “on the fly” map production

Args:

obs (Observation): container for tod. If the tod is

not required, you can avoid allocating obs.tod i.e. in lbs.Observation use allocate_tod=False.

pointings (float): pointing for each sample and detector

generated by lbs.get_pointings() Optional if already allocated in obs. When generating pointing information, set the variable hwp to None since the hwp rotation angle is added to the orientation angle within the fill_tod function.

hwp_radpsec (float): hwp rotation speed in radiants per second

save_tod (bool): if True, tod is saved in obs.tod; if False,

tod gets deleted

make_map(obss)

It generates “on the fly” map. This option is only availabe if built_map_on_the_fly is set to True.

Parameters:
  • class (obss list of) – Observations: only necessary for the communicator

  • pointings (float) – pointing for each sample and detector generated by lbs.get_pointings

  • hwp_radpsec (float) – hwp rotation speed in radiants per second

Returns:

rebinned T,Q,U maps

Return type:

map (float)

Bandpass template

litebird_sim.hwp_sys.bandpass_template_module.lowpass_chebyshev(freqs, f0, order=3, ripple_dB=0.2)

Define a lowpass with chebyshev prototype

Parameters:
  • freqs (np.array) – frequency in GHz

  • f0 (float) – low-frequency edge of the band in GHz

  • order (float) – chebyshev filter order (default: 3)

  • ripple_dB (float) – maximum ripple amplitude in decibels (default: 0.2)

Returns:

The bandpass transmission

litebird_sim.hwp_sys.bandpass_template_module.find_central_frequency(freqs, bandpass)
Find the effective central frequency of

a bandpass profile as defined in https://arxiv.org/abs/1303.5070

Parameters:
  • freqs (np.array) – frequency in GHz

  • bandpass (np.array) – transmission profile

Returns:

central frequency

litebird_sim.hwp_sys.bandpass_template_module.add_high_frequency_transmission(freqs, bandpass, location=3, transmission=0.5)

Add high frequency leakage

Parameters:
  • freqs (np.array) – frequency in GHz

  • bandpass (np.array) – transmission profile

  • location (float) – multiple of the central frequency of the bandpass profile where add the leakage (default: 3)

  • transmission (float) – relative amplitude of the high frequency leakage with respect to the nominal band (default: 0.5)

Returns:

Frequency range and bandpass transmission

litebird_sim.hwp_sys.bandpass_template_module.beam_throughtput(freqs)

Beam throughtput factor

Parameters:

freqs (np.array) – frequency in GHz

Returns:

\(1/\nu^2\)

litebird_sim.hwp_sys.bandpass_template_module.bandpass_profile(freqs: array | None = None, bandpass: dict | None = None, include_beam_throughput: bool | None = None)

Returns the bandpass profile either reading it from a file or computing it through the class BandPassInfo

Parameters:
  • freqs (np.array) – the array of frequencies considered

  • bandpass (dict) – dictionary with either the key “bandpass_file” pointing to a txt file with a column of freqs and a column of bandpass profile (attention: the frequencies read from the file have to coincide with the input array freqs) or with a key “band_type” (same naming convention as bandtype in BandPassInfo). In the second case, you should also define the edges of the bandpass, band_high_edge and band_low_edge, that would otherwise be set to the lowest and highest value of freqs.

  • include_beam_throughput (bool) – if True, divides by \(\nu^2\)

Returns:

Array of frequency range (same as freqs) and bandpass profile