Dipole anisotropy

The LiteBIRD Simulation Framework provides tools to simulate the signal associated with the relative velocity between the spacecraft’s rest frame with respect to the CMB. The motion of the spacecraft in the rest frame of the CMB is the composition of several components:

  1. The motion of the spacecraft around L2;

  2. The motion of the L2 point in the Ecliptic plane;

  3. The motion of the Solar System around the Galactic Centre;

  4. The motion of the Milky Way.

Components 1 and 2 are simulated by the LiteBIRD Simulation Framework using appropriate models for the motions, while components 3 and 4 are included using the Sun velocity derived by the solar dipole measured by the Planck satellite.

The motion of the spacecraft around L2 is modelled using a Lissajous orbit similar to what was used for the WMAP experiment [1], and it is encoded using the SpacecraftOrbit class.

Position and velocity of the spacecraft

The class SpacecraftOrbit describes the orbit of the LiteBIRD spacecraft with respect to the Barycentric Ecliptic Reference Frame and the motion of Barycentric Ecliptic Reference Frame with respect to the CMB; this class is necessary because the class ScanningStrategy (see the chapter Scanning strategy) only models the direction each instrument is looking at but knows nothing about the velocity of the spacecraft itself.

The class SpacecraftOrbit is a dataclass that is able to initialize its members to sensible default values, which are taken from the literature. As the LiteBIRD orbit around L2 is not fixed yet, the code assumes a WMAP-like Lissajous orbit, [1]. For the Sun velocity it assumes Planck 2018 solar dipole [6].

To compute the position/velocity of the spacecraft, you call spacecraft_pos_and_vel(); it requires either a time span or an Observation object, and it returns an instance of the class SpacecraftPositionAndVelocity:

import litebird_sim as lbs
from astropy.time import Time

orbit = lbs.SpacecraftOrbit(start_time=Time("2023-01-01"))
posvel = lbs.spacecraft_pos_and_vel(
    orbit,
    start_time=orbit.start_time,
    time_span_s=86_400.0,  # One day
    delta_time_s=3600.0    # One hour
)

print(posvel)

SpacecraftPositionAndVelocity(start_time=2023-01-01 00:00:00.000, time_span_s=86400.0, nsamples=25)

The output of the script shows that 25 «samples» have been computed; this means that the posvel variable holds information about 25 position/velocity pairs evenly spaced between 2023-01-01 and 2023-01-02: one at midnight, one at 1:00, etc., till midnight 2023-01-02. The SpacecraftPositionAndVelocity class keeps the table with the positions and the velocities in the fields positions_km and velocities_km_s, respectively, which are arrays of shape (nsamples, 3).

Here is a slightly more complex example that shows how to plot the distance between the spacecraft and the Sun as a function of time, as well as its speed. The latter quantity is of course most relevant when computing the CMB dipole.

import numpy as np
import litebird_sim as lbs
from astropy.time import Time
import matplotlib.pylab as plt

orbit = lbs.SpacecraftOrbit(start_time=Time("2023-01-01"))

posvel = lbs.spacecraft_pos_and_vel(
    orbit,
    start_time=orbit.start_time,
    time_span_s=3.15e7,  # One year
    delta_time_s=86_400.0,  # One day
)

# posvel.positions_km is a N×3 array containing the XYZ position
# of the spacecraft calculated every day for one year. We compute
# the modulus of the position, which is of course the
# Sun-LiteBIRD distance.
sun_distance_km = np.linalg.norm(posvel.positions_km, axis=1)

# We do the same with the velocities
speed_km_s = np.linalg.norm(posvel.velocities_km_s, axis=1)

# Plot distance and speed as functions of time
_, ax = plt.subplots(nrows=2, ncols=1, figsize=(7, 7))

ax[0].plot(sun_distance_km)
ax[0].set_xlabel("Day")
ax[0].set_ylabel("Distance [km]")

ax[1].plot(speed_km_s)
ax[1].set_xlabel("Day")
ax[1].set_ylabel("Speed [km/s]")

(Source code, png, hires.png, pdf)

_images/spacecraft_demo.png

Computing the dipole

The CMB dipole is caused by a Doppler shift of the frequencies observed while looking at the CMB blackbody spectrum, according to the formula

(1)\[T(\vec\beta, \hat n) = \frac{T_0}{\gamma \bigl(1 - \vec\beta \cdot \hat n\bigr)},\]

where \(T_0\) is the temperature in the rest frame of the CMB, \(\vec \beta = \vec v / c\) is the dimensionless velocity vector, \(\hat n\) is the direction of the line of sight, and \(\gamma = \bigl(1 - \vec\beta \cdot \vec\beta\bigr)^2\).

However, CMB experiments usually employ the linear thermodynamic temperature definition, where temperature differences \(\Delta_1 T\) are related to the actual temperature difference \(\Delta T\) by the relation

(2)\[\Delta_1 T = \frac{T_0}{f(x)} \left(\frac{\mathrm{BB}(T_0 + \Delta T)}{\mathrm{BB}(T_0)} - 1\right) = \frac{T_0}{f(x)} \left(\frac{\exp x - 1}{\exp\left(x\frac{T_0}{T_0 + \Delta T}\right) - 1} - 1\right),\]

where \(x = h \nu / k_B T\),

\[f(x) = \frac{x e^x}{e^x - 1},\]

and \(\mathrm{BB}(\nu, T)\) is the spectral radiance of a black-body according to Planck’s law:

\[\mathrm{BB}(\nu, T) = \frac{2h\nu^3}{c^2} \frac1{e^{h\nu/k_B T} - 1} = \frac{2h\nu^3}{c^2} \frac1{e^x - 1}.\]

There is no numerical issue in computing the full formula, but often models use some simplifications, to make the math easier to work on the blackboard. The LiteBIRD Simulation Framework implements several simplifications of the formula, which are based on a series expansion of (1); the caller must pass an object of type DipoleType (an enum class), whose value signals which kind of approximation to use:

  1. The most simple formula uses a series expansion of (1) at the first order:

    \[\Delta T(\vec\beta, \hat n) = T_0 \vec\beta\cdot\hat n,\]

    which is associated to the constant DipoleType.LINEAR.

  2. The same series expansion for (1), but stopped at the second order (DipoleType.QUADRATIC_EXACT):

    \[\Delta T(\vec\beta, \hat n) = T_0\left(\vec\beta\cdot\hat n + \bigl(\vec\beta\cdot\hat n\bigr)^2\right),\]

    which discards a \(-T_0\,\beta^2/2\) term (monopole).

  3. The exact formula as in (1) (DipoleType.TOTAL_EXACT).

  4. Using a series expansion to the second order of (2) instead of (1) and neglecting monopoles (DipoleTotal.QUADRATIC_FROM_LIN_T):

    \[\Delta_2 T(\nu) = T_0 \left(\vec\beta\cdot\hat n + q(x) \bigl(\vec\beta\cdot\hat n\bigr)^2\right),\]

    where the dependence on the frequency ν is due to the presence of the term \(x = h\nu / k_B T\) in the equation. This is the formula to use if you want the leading frequency-dependent term (second order) without the boosting induced monopoles.

  5. Finally, linearizing (1) through (2) (DipoleTotal.TOTAL_FROM_LIN_T):

    \[\Delta T = \frac{T_0}{f(x)} \left(\frac{\mathrm{BB}\left(T_0 / \gamma\bigl(1 - \vec\beta\cdot\hat n\bigr)\right)}{\mathrm{BB}(T_0)} - 1\right) = \frac{T_0}{f(x)} \left(\frac{\mathrm{BB}\bigl(\nu\gamma(1-\vec\beta\cdot\hat n), T_0\bigr)}{\bigl(\gamma(1-\vec\beta\cdot\hat n)\bigr)^3\mathrm{BB}(T_0)}\right).\]

    In this case too, the temperature variation depends on the frequency because of (2). This is the formula that is typically used by CMB experiments.

You can add the dipole signal to an existing TOD through the function add_dipole_to_observations(), as the following example shows:

from astropy.time import Time
import numpy as np
import litebird_sim as lbs
import matplotlib.pylab as plt

start_time = Time("2022-01-01")
time_span_s = 120.0  # Two minutes
sampling_hz = 20

sim = lbs.Simulation(start_time=start_time, duration_s=time_span_s, random_seed=12345)

# We pick a simple scanning strategy where the spin axis is aligned
# with the Sun-Earth axis, and the spacecraft spins once every minute
sim.set_scanning_strategy(
    lbs.SpinningScanningStrategy(
        spin_sun_angle_rad=np.deg2rad(0),
        precession_rate_hz=0,
        spin_rate_hz=1 / 60,
        start_time=start_time,
    ),
    delta_time_s=5.0,
)

# We simulate an instrument whose boresight is perpendicular to
# the spin axis.
sim.set_instrument(
    lbs.InstrumentInfo(
        boresight_rotangle_rad=0.0,
        spin_boresight_angle_rad=np.deg2rad(90),
        spin_rotangle_rad=np.deg2rad(75),
    )
)

# A simple detector looking along the boresight direction
det = lbs.DetectorInfo(
    name="Boresight_detector",
    sampling_rate_hz=sampling_hz,
    bandcenter_ghz=100.0,
)

(obs,) = sim.create_observations(detectors=[det])

sim.compute_pointings()

# Simulate the orbit of the spacecraft and compute positions and
# velocities
orbit = lbs.SpacecraftOrbit(obs.start_time)
pos_vel = lbs.spacecraft_pos_and_vel(orbit, obs, delta_time_s=60.0)

t = obs.get_times()
t -= t[0]  # Make `t` start from zero

# Create two plots: the first shows the full two-minute time span, and
# the second one shows a zoom over the very first points of the TOD.
_, axes = plt.subplots(nrows=2, ncols=1, figsize=(10, 12))

# Make a plot of the TOD using all the dipole types
for type_idx, dipole_type in enumerate(lbs.DipoleType):
    obs.tod[0][:] = 0  # Reset the TOD
    lbs.add_dipole_to_observations(obs, pos_vel, dipole_type=dipole_type)

    axes[0].plot(t, obs.tod[0], label=str(dipole_type))
    axes[1].plot(t[0:20], obs.tod[0][0:20], label=str(dipole_type))

axes[0].set_xlabel("Time [s]")
axes[0].set_ylabel("Signal [K]")
axes[1].set_xlabel("Time [s]")
axes[1].set_ylabel("Signal [K]")
axes[1].legend()

(Source code)

The example plots two minutes of a simulated timeline for a very simple instrument, and it zooms over the very first points to show that there is indeed some difference in the estimate provided by each method.

Methods of class simulation

The class Simulation provides two simple functions that compute poisition and velocity of the spacescraft Simulation.compute_pos_and_vel(), and add the solar and orbital dipole to all the observations of a given simulation Simulation.add_dipole().

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

start_time = Time("2025-01-01")
time_span_s = 1000.0
sampling_hz = 10.0

sim = lbs.Simulation(
    start_time=start_time,
    duration_s=time_span_s,
    random_seed=12345,
)

# We pick a simple scanning strategy where the spin axis is aligned
# with the Sun-Earth axis, and the spacecraft spins once every minute
sim.set_scanning_strategy(
    lbs.SpinningScanningStrategy(
        spin_sun_angle_rad=np.deg2rad(0),
        precession_rate_hz=0,
        spin_rate_hz=1 / 60,
        start_time=start_time,
    ),
    delta_time_s=5.0,
 )

# We simulate an instrument whose boresight is perpendicular to
# the spin axis.
sim.set_instrument(
    lbs.InstrumentInfo(
        boresight_rotangle_rad=0.0,
        spin_boresight_angle_rad=np.deg2rad(90),
        spin_rotangle_rad=np.deg2rad(75),
    )
)

# A simple detector looking along the boresight direction
det = lbs.DetectorInfo(
    name="Boresight_detector",
    sampling_rate_hz=sampling_hz,
    bandcenter_ghz=100.0,
)

sim.create_observations(detectors=det)

sim.compute_pointings()

sim.compute_pos_and_vel()

sim.add_dipole()

for i in range(5):
    print(f"{sim.observations[0].tod[0][i]:.5e}")

3.44963e-03
3.45207e-03
3.45413e-03
3.45582e-03
3.45712e-03

Note that even if Simulation.compute_pos_and_vel() is not explicitly invoked, Simulation.add_dipole() takes care of that internally initializing SpacecraftOrbit and computing positions and velocities.

API reference

class litebird_sim.spacecraft.SpacecraftOrbit(start_time: Time, earth_l2_distance_km: float = 1496509.30522, radius1_km: float = 244450.0, radius2_km: float = 137388.0, ang_speed1_rad_s: float = 4.0266897619299316e-07, ang_speed2_rad_s: float = 3.955023673808657e-07, phase_rad: float = -0.8367806565761614, solar_velocity_km_s: float = 369.816, solar_velocity_gal_lat_rad: float = 0.842173724, solar_velocity_gal_lon_rad: float = 4.6080357444)

Bases: object

A dataclass describing the orbit of the LiteBIRD spacecraft

This structure has the following fields:

  • start_time: Date and time when the spacecraft starts its nominal orbit

  • earth_l2_distance_km: distance between the Earth’s barycenter and the L2 point, in km

  • radius1_km: first radius describing the Lissajous orbit followed by the spacecraft, in km

  • radius2_km: second radius describing the Lissajous orbit followed by the spacecraft, in km

  • ang_speed1_rad_s: first angular speed of the Lissajous orbit, in rad/s

  • ang_speed2_rad_s: second angular speed of the Lissajous orbit, in rad/s

  • phase_rad: phase difference between the two periodic motions in the Lissajous orbit, in radians

  • solar_velocity_km_s: velocity of the Sun as estimated from Planck 2018 Solar dipole (see arxiv: 1807.06207)

  • solar_velocity_gal_lat_rad: galactic latitude direction of the Planck 2018 Solar dipole

  • solar_velocity_gal_lon_rad: galactic longitude direction of the Planck 2018 Solar dipole

The default values are the nominal numbers of the orbit followed by WMAP, described in Cavaluzzi, Fink & Coyle (2008).

ang_speed1_rad_s: float = 4.0266897619299316e-07
ang_speed2_rad_s: float = 3.955023673808657e-07
earth_l2_distance_km: float = 1496509.30522
phase_rad: float = -0.8367806565761614
radius1_km: float = 244450.0
radius2_km: float = 137388.0
solar_velocity_ecl_xyz_km_s = array([-359.00346797,   52.57540642,  -71.52769031])
solar_velocity_gal_lat_rad: float = 0.842173724
solar_velocity_gal_lon_rad: float = 4.6080357444
solar_velocity_km_s: float = 369.816
start_time: Time
class litebird_sim.spacecraft.SpacecraftPositionAndVelocity(orbit: SpacecraftOrbit, start_time: Time, time_span_s: float, positions_km=None, velocities_km_s=None)

Bases: object

Encode the position/velocity of the spacecraft with respect to the Solar System

This class contains information that characterize the motion of the spacecraft. It is mainly useful to simulate the so-called CMB «orbital dipole» and to properly check the visibility of the Sun, the Moon and inner planets. The coordinate system used by this class is the standard Barycentric Ecliptic reference frame.

The fields of this class are the following:

  • orbit: a SpacecraftOrbit object used to compute the positions and

    velocities in this object;

  • start_time: the time when the nominal orbit started (a astropy.time.Time

    object);

  • time_span_s: the time span covered by this object, in seconds;

  • positions_km: a N×3 matrix, representing a list of N XYZ vectors

    encoding the position of the spacecraft in the Barycentric Ecliptic reference frame (in kilometers);

  • velocities_km_s: a N×3 matrix, representing a list of N XYZ vectors

    encoding the linear velocity of the spacecraft in the Barycentric Ecliptic reference frame (in km/s).

compute_velocities(time0: Time, delta_time_s: float, num_of_samples: int)

Perform a linear interpolation to sample the satellite velocity in time

Return a N×3 array containing a set of num_of_samples 3D vectors with the velocity of the spacecraft (in km/s) computed every delta_time_s seconds starting from time time0.

litebird_sim.spacecraft.compute_l2_pos_and_vel(time0: Time, earth_l2_distance_km: float = 1496509.30522)

Compute the position and velocity of the L2 Sun-Earth point at a given time.

The L2 point is not calculated using Lagrange’s equations; instead, its distance from the Earth must be provided as the parameter earth_l2_distance_km. The default value is a reasonable estimate. The L2 point is assumed to lie along the line that connects the Solar System Barycenter with the Earth’s gravitational center.

The return value is a 2-tuple containing two NumPy arrays:

  1. A 3D array containing the XYZ components of the vector specifying the position of the L2 point, in km

  2. A 3D array containing the XYZ components of the velocity vector of the L2 point, in km/s

The two vectors are always roughly perpendicular, but they are not exactly due to the irregular motion of the Earth (caused by gravitational interactions with other Solar System bodies, like the Moon and Jupiter).

litebird_sim.spacecraft.compute_lissajous_pos_and_vel(time0, earth_angle_rad, earth_ang_speed_rad_s, radius1_km, radius2_km, ang_speed1_rad_s, ang_speed2_rad_s, phase_rad)

Compute the position and velocity of the spacecraft assuming a Lissajous orbit

The position and velocity are calculated in a reference frame centered on the L2 point, whose axes are aligned with the Solar System Barycenter. This means that the position and velocity of the spacecraft with respect to the Solar System Barycenter itself can be calculated by summing the result of this function with the result of a call to compute_l2_pos_and_vel().

litebird_sim.spacecraft.compute_start_and_span_for_obs(obs: Observation | List[Observation]) → Tuple[Time, float]

Compute the start time and the overall duration in seconds of a set of observations.

The code returns the earliest start time of the observations in obs as well as their overall time span. Gaps between observations are neglected.

litebird_sim.spacecraft.cycles_per_year_to_rad_per_s(x: float) → float

Convert an angular speed from cycles/yr to rad/s

litebird_sim.spacecraft.get_ecliptic_vec(vec)

Convert a coordinate in a XYZ vector expressed in the Ecliptic rest frame

litebird_sim.spacecraft.spacecraft_pos_and_vel(orbit: SpacecraftOrbit, obs: Observation | List[Observation] | None = None, start_time: Time | None = None, time_span_s: float | None = None, delta_time_s: float = 86400.0) → SpacecraftPositionAndVelocity

Compute the position and velocity of the L2 point within some time span

This function computes the XYZ position and velocity of the second Sun-Earth Lagrangean point (L2) over a time span specified either by a Observation object/list of objects, or by an explicit pair of values start_time (an astropy.time.Time object) and time_span_s (length in seconds). The position is specified in the standard Barycentric Ecliptic reference frame.

The position of the L2 point is computed starting from the position of the Earth and moving away along the anti-Sun direction by a number of kilometers equal to earth_l2_distance_km.

The result is an object of type SpacecraftPositionAndVelocity.

If SpacecraftOrbit.solar_velocity_km_s > 0 also the Sun velocity in the rest frame of the CMB is added to the total velocity of the spacecraft.

litebird_sim.spacecraft.sum_lissajous_pos_and_vel(pos, vel, start_time_s, end_time_s, radius1_km, radius2_km, ang_speed1_rad_s, ang_speed2_rad_s, phase_rad)

Add the position and velocity of a Lissajous orbit to some positions/velocities

The pos and vel arrays must have the same shape (N×3, with N being the number of position-velocity pairs); the 3D position and velocity of the Lissajous orbit will be added to these arrays.

The times array must be an array of float values, typically starting from zero, and they must measure the time in seconds. All the other parameters are usually taken from a SpacecraftOrbit object.

This function has no return value, as the result is stored in pos and vel.

class litebird_sim.dipole.DipoleType(value)

Bases: IntEnum

Approximation for the Doppler shift caused by the motion of the spacecraft

LINEAR = 0

Linear approximation in β using thermodynamic units:

\[\Delta T(\vec\beta, \hat n) = T_0 \vec\beta\cdot\hat n\]
QUADRATIC_EXACT = 1

Second-order approximation in β using thermodynamic units:

\[\Delta T(\vec\beta, \hat n) = T_0\left(\vec\beta\cdot\hat n + \bigl(\vec\beta\cdot\hat n\bigr)^2\right)\]
QUADRATIC_FROM_LIN_T = 3

Second-order approximation in β using linearized units:

\[\Delta_2 T(\nu) = T_0 \left(\vec\beta\cdot\hat n + q(x) \bigl(\vec\beta\cdot\hat n\bigr)^2\right)\]
TOTAL_EXACT = 2

Exact formula in true thermodynamic units:

\[\frac{T_0}{\gamma \bigl(1 - \vec\beta \cdot \hat n\bigr)}\]
TOTAL_FROM_LIN_T = 4

Full formula in linearized units (the most widely used):

\[\Delta T = \frac{T_0}{f(x)} \left(\frac{\mathrm{BB}\left(T_0 / \gamma\bigl(1 - \vec\beta\cdot\hat n\bigr)\right)}{\mathrm{BB}(T_0)} - 1\right) = \frac{T_0}{f(x)} \left(\frac{\mathrm{BB}\bigl(\nu \gamma(1-\vec\beta\cdot\hat n), T_0\bigr)}{\bigl(\gamma(1- \vec\beta\cdot\hat n)\bigr)^3\mathrm{BB}(t_0)}\right).\]
litebird_sim.dipole.add_dipole(tod, pointings, velocity, t_cmb_k: float, frequency_ghz: ndarray, dipole_type: DipoleType)

Add the CMB dipole to some time-ordered data

This functions modifies the values in tod by adding the contribution of the CMB dipole. Use dipole_type to specify which kind of approximation to use for the dipole component. The pointings argument must be a N×3 matrix containing the pointing information, where N is the size of the tod array. The velocity argument is usually computed through spacecraft_pos_and_vel(). Finally, t_cmb_k is the temperature of the monopole and frequency_ghz is an array containing the frequencies of each detector in the TOD.

litebird_sim.dipole.add_dipole_for_one_detector(tod_det, theta_det, phi_det, velocity, t_cmb_k, nu_hz, f_x, q_x, dipole_type: DipoleType)
litebird_sim.dipole.add_dipole_to_observations(obs: Observation | List[Observation], pos_and_vel: SpacecraftPositionAndVelocity, pointings: ndarray | List[ndarray] | None = None, t_cmb_k: float = 2.72548, dipole_type: DipoleType = DipoleType.TOTAL_FROM_LIN_T, frequency_ghz: None | ndarray = None, component: str = 'tod')

Add the CMB dipole to some time-ordered data

This is a wrapper around the add_dipole() function that applies to the TOD stored in obs, 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 dipole alone
    cur_obs.dipole_tod = np.zeros_like(cur_obs.tod)

# Ask `add_dipole_to_observations` to store the dipole
# in `obs.dipole_tod`
add_dipole_to_observations(sim.observations, component="dipole_tod")
litebird_sim.dipole.calculate_beta(theta, phi, v_km_s)

Return a 2-tuple containing β·n and β

litebird_sim.dipole.compute_dipole_for_one_sample_linear(theta, phi, v_km_s, t_cmb_k)
litebird_sim.dipole.compute_dipole_for_one_sample_quadratic_exact(theta, phi, v_km_s, t_cmb_k)
litebird_sim.dipole.compute_dipole_for_one_sample_quadratic_from_lin_t(theta, phi, v_km_s, t_cmb_k, q_x)
litebird_sim.dipole.compute_dipole_for_one_sample_total_exact(theta, phi, v_km_s, t_cmb_k)
litebird_sim.dipole.compute_dipole_for_one_sample_total_from_lin_t(theta, phi, v_km_s, t_cmb_k, nu_hz, f_x, planck_t0)
litebird_sim.dipole.compute_scalar_product(theta, phi, v)

Return the scalar (dot) product between a given direction and a velocity

litebird_sim.dipole.planck(nu_hz, t_k)

Return occupation number at frequency nu_hz and temperature t_k