Non-linearity injection

Non-linearity is the effect of a non-ideal TES detector response, where the responsivity \(S\) is not constant as the optical power varies. The LiteBIRD Simulation Framework provides a non-linearity simulation module to simulate the effect of non-linearity on TODs.

The framework provides the simplest case, which is quadratic non-linearity. This case is described in Micheli+2024, where the effect of non-linearity is propagated to the estimation of the tensor-to-scalar ratio.

Considering a first order correction of the usual linear gain, a TOD \(d(t)\) is modified according to:

\[d(t) = [1+g_1 d(t)] d(t)\]

where \(g_1\) is the detector non-linearity factor in units of \(K^{-1}\).

To simulate a quadratic non-linearity, one can use the method of Simulation class Simulation.apply_quadratic_nonlin(), or any of the low-level functions: apply_quadratic_nonlin_to_observations(), apply_quadratic_nonlin_for_one_detector(). The following example shows the typical usage of the method and low-level functions:

import numpy as np
import litebird_sim as lbs
from astropy.time import Time

start_time = Time("2025-02-02T00:00:00")
mission_time_days = 1
sampling_hz = 19

dets = [
    lbs.DetectorInfo(
        name="det_A", sampling_rate_hz=sampling_hz
    ),
    lbs.DetectorInfo(
        name="det_B", sampling_rate_hz=sampling_hz
    ),
]

sim = lbs.Simulation(
    base_path="nonlin_example",
    start_time=start_time,
    duration_s=mission_time_days * 24 * 3600.0,
    random_seed=12345,
)

sim.create_observations(
    detectors=dets,
    split_list_over_processes=False,
)

# Creating fiducial TODs
sim.observations[0].nl_2_self = np.ones_like(sim.observations[0].tod)
sim.observations[0].nl_2_obs = np.ones_like(sim.observations[0].tod)
sim.observations[0].nl_2_det = np.ones_like(sim.observations[0].tod)

If the non-linearity parameter is not read from the IMo, one has to specify \(g_1\) using the g_one_over_k argument as in the following example:

# Define non-linear parameters for the detectors. We choose the same
# value for both detectors in this example, but it is not necessary.
sim.observations[0].g_one_over_k = np.ones(len(dets)) * 1e-3

# Applying non-linearity using the `Simulation` class method
sim.apply_quadratic_nonlin(component = "nl_2_self",)

# Applying non-linearity on the given TOD component of an `Observation`
# object
lbs.non_linearity.apply_quadratic_nonlin_to_observations(
    observations=sim.observations,
    component="nl_2_obs",
)

# Applying non-linearity on the TOD arrays of the individual detectors.
for idx, tod in enumerate(sim.observations[0].nl_2_det):
    lbs.non_linearity.apply_quadratic_nonlin_for_one_detector(
        tod_det=tod,
        g_one_over_k=sim.observations[0].g_one_over_k[idx],
    )

In particular, the effect of detector non-linearity must be included to assess its impact when coupled with other systematic effects. As described in Micheli+2024, a typical case is the coupling with HWP synchronous signal (HWPSS) appearing at twice the rotation frequency of the HWP. This kind of signal can be produced by non-idealities of the HWP, such as its differential transmission and emission.

In that case, the usual TOD \(d(t)\) will contain an additional term, and can be written as:

\[d(t) = d(t) = I + \mathrm{Re}[\epsilon_{\mathrm{pol}}e^{4i\chi}(Q+iU)]+A_2 \cos(2 \omega_{HWP} t) + N\]

where \(A_2\) is the amplitude of the HWPSS and \(\omega_{HWP}\) is the rotation speed of the HWP. In presence of detector non-linearity, the 2f signal is up-modulated to 4f, affecting the science band.

The framework provides an independent module to introduce this signal in the simulation, adding it to the TODs. To simulate the 2f signal from a rotating, non-ideal 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().

If the 2f amplitude is not read from the IMo, one has to specify \(A_2\) using the amplitude_2f_k argument. The argument optical_power_k enables the inclusion of the integrated nominal optical power expected for each channel. See the following example:

import numpy as np
import litebird_sim as lbs
from astropy.time import Time
from litebird_sim.pointings import get_hwp_angle


telescope = "LFT"
channel = "L4-140"
detlist = ["000_001_017_QB_140_T", "000_001_017_QB_140_B"]
imo_version = "vPTEP"
start_time = Time("2025-02-02T00:00:00")
mission_time_days = 1

imo = lbs.Imo(flatfile_location=lbs.PTEP_IMO_LOCATION)


sim = lbs.Simulation(
    base_path="nonlin_example",
    start_time=start_time,
    imo=imo,
    duration_s=mission_time_days * 24 * 3600.0,
    random_seed=12345,
)

# Load the definition of the instrument
sim.set_instrument(
    lbs.InstrumentInfo.from_imo(
        imo,
        f"/releases/{imo_version}/satellite/{telescope}/instrument_info",
    )
)

dets = []
for n_det in detlist:
    det = lbs.DetectorInfo.from_imo(
        url=f"/releases/{imo_version}/satellite/{telescope}/{channel}/{n_det}/detector_info",
        imo=imo,)
    det.sampling_rate_hz = 1
    dets.append(det)

sim.create_observations(
    detectors=dets,
    split_list_over_processes=False,
)

sim.set_scanning_strategy(imo_url=f"/releases/{imo_version}/satellite/scanning_parameters/")

sim.set_hwp(
    lbs.IdealHWP(sim.instrument.hwp_rpm * 2 * np.pi / 60,),
)

sim.prepare_pointings()
sim.precompute_pointings()

# Creating fiducial TODs
sim.observations[0].tod_2f = np.zeros_like(sim.observations[0].tod)

# Define differential emission parameters for the detectors.
sim.observations[0].amplitude_2f_k = np.array([0.1, 0.1])
sim.observations[0].optical_power_k = np.array([1.0, 1.0])

# Adding 2f signal from HWP differential emission using the `Simulation` class method
sim.add_2f(component="tod_2f")

Refer to the API reference for the full list of non-linearity simulation parameters.

API reference

litebird_sim.non_linearity.apply_quadratic_nonlin(tod, g_one_over_k: float)

Apply a quadratic nonlinearity to some time-ordered data

This functions modifies the values in tod by adding a parabolic function of the TOD itself. The g_one_over_k argument must be a N_dets array containing the amplitude of the nonlinear gain to be applied.

litebird_sim.non_linearity.apply_quadratic_nonlin_for_one_detector(tod_det, g_one_over_k)

litebird_sim.non_linearity.apply_quadratic_nonlin_for_one_sample(data, g_one_over_k)

litebird_sim.non_linearity.apply_quadratic_nonlin_to_observations(observations: Observation | List[Observation], component: str = 'tod', g_one_over_k: float | None = None)

Apply a quadratic nonlinearity to some time-ordered data

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

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

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

# Ask `apply_quadratic_nonlin_to_observations` to store the nonlinear TOD
# in `observations.nl_tod`
apply_quadratic_nonlin_to_observations(sim.observations, component="nl_tod")