Aberration v2

This is a small presentation about aberration and modulation of CMB maps, and a new and improved implementation of this in pixell.

What is aberration?

Aberration is the apparent deflection of light due to the observers velocity relative to the source. It remaps directions such that \begin{align} z' = \frac{z+\beta}{1+z\beta}, \quad z = \frac{z'-\beta}{1-z'\beta} \end{align} where $z$|$z'$ is the cosine of the true|observed angle from RA = 167.919°, dec = -6.936°, our direction of motion with respect to the CMB; $\beta = 0.001235$ is our relative speed, and $\gamma = 1/\sqrt{1-\beta^2}$ (see Planck 2013 XXVII). We can therefore simulate the effect of aberration in our maps by:

  1. Get the coordinates of each pixel in our output map. This is a simple function of our projection.
  2. Use the aberration formula to calculate the corresponding true coordinates.
  3. Read off the values at those locations in a map of the unaberrated sky, e.g. the output from a CMB simulation.
  4. (For polarization, apply a Q/U rotation to take into account the effect of parallel transport)

Interpolation troubles

Despite all the trigonometry involved in steps 1, 2 and 4, the main difficulty with this procedure turns out to be step 3. The problem is that the locations we need in step 3 rarely correspond to a pixel center in our unaberrated map, so we have to interpolate from nearby pixels. In previous versions of pixell, this was implemented using bicubic spline interpolation, which is a commonly used high-quality interpolation. The result of this can be seen below. Press j/k to switch between the raw map and the aberrated map.

The deflection is quite subtle, only a few arcmin, but as you can see, CMB fetures appear to towards our direction of travel (near the middle left edge of the map) and away from its antipode (near the center of the map). The deflected map looks sensible by eye. But if we repeat the test with a map with more small-scale power, we notice an issue (here with a 10x bigger beta since I'm using smaller pixels here). Press j/k to cycle between the true sky and the aberrated one.

As is hopefully visible, the aberrated map has less power than the raw map. It appears to be slightly smoothed. It's a bit subtle by eye, but it's dramatic in the power spectrum (blue curve). Near the Nyquist limit the power loss reaches more than 30%!

Non-uniform FFT to the rescue

ducc0, which I already use for spherical harmonics transforms, recently added support for non-uniform fast fourier transforms, and used them to implement fast and accurate lensing simulations. I used these to implement a new aberration interpolator in pixell, which turned out to be 4x as fast, use much less memory, and most importantly solved the power loss, as shown in the red curve above. Press j/k to cycle between the raw map, the nuFFT interpolated aberrated map and the old aberrated map.

Modulation

Aside from displacing positions on the sky, our speed relative to the CMB restframe also results in a doppler shift. \begin{align} T' = \frac{T}{\gamma(1-z'\beta)} \end{align} where $T$ is the original temperature and $T'$ is the observed one. The observed temperature is increased in the direction of motion (z' close to 1) and decreased in the other direction, as expected from a doppler shift. To first order this multiplies the raw sky map by a dipole, creating the 3.4 mK CMB dipole from the 2.7 K CMB monopole, but all multipoles are affected.

If we ignore the contribution from the monopole, then the effect of modulation is quite subtle. Press j/k to step through gradually larger doppler boosts below.

Map temperature units

The effects of modulation are simple in terms of the blackbody temperature, but CMB telescopes don't actually observe the temperature, they observe sky's spectral radiance and just choose to express that in temperature-like units. These linearized CMB temperature units are defined by dividing the radiance excess relative to the monopole by the derivative of the Planck law with respect to temperature: \begin{align} \Delta \tilde T \equiv \frac{I_\nu-B_\nu(T_0)}{\frac{dB_\nu}{dT}(T_0)} \end{align} where $\Delta \tilde T$ is the map "temperature", $I_\nu$ is the sky's spectral radiance at frequency $\nu$, $B_\nu$ is CMB blackbody spectrum and $T_0$ is the CMB monopole temperature. If the source of the radiance is the CMB itself, this becomes \begin{align} \Delta \tilde T \equiv \frac{B_\nu(T)-B_\nu(T_0)}{\frac{dB_\nu}{dT}(T_0)} \end{align} This can be Taylor-expanded as they do in the Planck paper. I used this approach in my old implementation, but with a bug resulting in a non-linear correction that was off by a factor of about 2. In my new implementation I simply use the full formula instead of Taylor expanding it. This comes at a speed cost, but is both simpler and more accurate and is still fast enough not to be a bottleneck. Note that the full formula consumes about 5 digits of numerical precision in polarization (2.7 K vs. 10 µK being added together), so this formula needs double precision which provides 16 digits of precision (as opposed to single precision's 7 digits). The plots below show the effect working in map temperature units. Press j/k to cycle through speeds and a/b to switch between map units, real temperature and plain aberration. The result is frequency-dependent, and this is evaluated at 150 GHz.

The reason why much of the sky seems to become empty with linearized map units is that the CMB in that direction becomes too cold, with most of the blackbody spectrum falling below the detector frequency (150 GHz in this case). Similarly, the other side of the sky has a greater boost than one would expect from simple modulation due to the effect of the high-frequency Blackbody fall-off being lessened.

What about polarization?

Polarization works much like total intensity, except

  1. It has no monopole-sourced terms, since there is no polarization monopole. (The monopole still enters into equation 4, though, since polarization is just the difference between two different subsets of the total signal.)
  2. It gets a small polarization rotation due to positions being displaced on the sky. This is implemented using a function from ducc, though this one only supports double precision, so it might be worth it to implement directly at some point.

In the plot below, press j/k to change β and a/s to switch between T and Q (U looks very similar), all for linearized map units.

Deaberration

In CMB analysis we observe the boosted (= aberrated + modulated) sky, but want to recover what it would look like in the CMB rest frame, since this is where the CMB should be a statistically isotropic field Gaussian described by a simple power spectrum (ignoring lensing). Aberration and modulation are both reversible operations, and the math works out such that applying aberration and modulation with the opposite sign for β undoes their effect. We can call these operations deaberration and demodulation, and the full operation deboosting.

Things get a bit more complicated when working with linearized map units due to their non-linear relationship with the blackbody temperature, but if one is careful one can invert this too.

The plot below shows the effect of first boosting and then deboosting a map with a relatively high β = 0.1228, a case which is around 100x as challenging as our real sky. Press j/k to step through the original map, the boosted map, and the deboosted map.

Usage

These functions are implemented in pixell.aberration. The main functions are aberration.boost_map and aberration.deboost_map. Here's come example code:

from pixell import aberration
# Apply a doppler boost to map_raw, and then recover it by
# deboosting it. map_deboost should be very similar to map_raw
map_obs     = aberration.boost_map(map_raw)
map_deboost = aberration.deboost_map(map_obs)
# Use temperature instead of linearized map units
map_plain   = aberration.boost_map(map_raw, modulation="plain")
# Use aberration only, no modulation at all
map_abonly  = aberration.boost_map(map_raw, modulation="none")
# Custom aberration magnitude
map_very    = aberration.boost_map(map_raw, beta=0.5)
# Custom aberration direction
map_north   = aberration.boost_map(map_raw, dir=[np.pi/2,0])