Band averaged polarised intensity

[AlecThomson]

Hey @joshoewahp,

Congrats on the v2 release!

I'm just poking around in the plotting and RM-synthy bits. Is it right that the linear pols are still extracted from the mean values across the band? If so, I'm happy to help get them from the FDF instead.

Are you able to direct me to the critical point where the band-averaged values are computed? Am I right in guessing its here?

dstools/dstools/dynamic_spectrum.py

Line 599 in 1e227ea

self.flux[stokes] = np.nanmean(ydata, axis=avg_axis)

[AlecThomson]

Oh also noting that it might be neat to add a de-biasing routine to counter the Ricean bias in polarised intensity

[joshoewahp]

Thanks!

At the moment they are extracted from a band-average after de-rotating at the most significant RM component in the FDF spectrum.

The current workflow prior to that step is:

• use RM-lite to compute the 2D FDF,
• select the maximum RM from the timesample at the Stokes I band-averaged lightcurve peak,
• de-rotate the complex Q + iU dynamic spectrum using that RM,
• form Q and U lightcurves by averaging over the frequency axis of the de-rotated Q and U DS (this is the step you've linked),
• form polarised intensity and PA lightcurves from those.

Very happy to get your input on a better approach!

One thing I'm keen to explore is how effective it is to allow each time-sample to select the ideal RM to allow for a time-variable RM as your RM-lite docs show. Though after a quick test with a ~200 mJy pulse, this seems to introduce too much RM noise into the recovered PA. Will follow up with some examples when I get a moment but would be good to chat further.

Oh also noting that it might be neat to add a de-biasing routine to counter the Ricean bias in polarised intensity

Good call, should be easy enough to add this in once we settle on a good approach for the above!

[AlecThomson]

So there's probably a few avenues to take, depending on what information you'd like dstools to ultimately extract.

Keeping things at the most basic, here's my line of thinking. RM-Synthesis effectively does the de-rotation for you (in an optimal sense) across all samples RMs. So all the information you need is already there after doing it. By taking the peak in the FDF (and assuming the Faraday simple case) we get the band-averaged value in-place without any worries of bandwidth depolarisation.

In general, I'm wary of doing a derotation. Or, at the very least, only presenting values from derotating. I think it'd generally be better to do a rotation after the fact, so that the user has control of how the derotation is done. This would allow for, e.g. a Galactic foreground correction that is constant with time. EDIT: Also, deriving an RM to derotate and subtract will add to the error budget of the residual spectrum. I'm guessing that's why your test had PAs going all over the place

Some things to consider along the way. It's generally more robust to the fit the peak of the FDF - I have some helper tools in rm_lite for this. The min/max/delta RM should also be selected based on the frequency coverage of a given observation as well to produce optimal results.

[krose098]

This is the way I am currently using rm_lite to calculate the corrected, de-biased polarised intensity from an DynamicSpectrum object :

import numpy as np
from rm_lite.utils.synthesis import make_phi_arr, freq_to_lambda2, rmsynth_nufft

data =  ds.data['Q'].real + 1j * ds.data['U'].real
data = data.T
freq_hz = ds.freq*1e6
phis = make_phi_arr(500, 0.1)
lam_sq_0_m2 = float(np.mean(freq_to_lambda2(freq_hz)))
psi0_time = ds.time % 180

fdf_spectrum = rmsynth_nufft(
    complex_pol_arr=data, 
    lambda_sq_arr_m2=freq_to_lambda2(freq_hz), 
    phi_arr_radm2=phis, 
    weight_arr=np.ones_like(freq_hz), 
    lam_sq_0_m2=lam_sq_0_m2,
)

peak_pi_spectrum = np.max(np.abs(fdf_spectrum), axis=0)
dy_Q = np.std(fdf_spectrum.real) 
dy_U = np.std(fdf_spectrum.imag)
noise_QU = (dy_Q + dy_U)/2

peak_pi_spectrum_corr = peak_pi_spectrum.copy()
peak_pi_spectrum_snr = peak_pi_spectrum / noise_QU
bias_correction_level = 5 # sigma
correction_idx = peak_pi_spectrum_snr > bias_correction_level
peak_pi_spectrum_corr[correction_idx] = np.sqrt(peak_pi_spectrum[correction_idx]**2 - 2.3 * noise_QU**2)

It should be fairly straightforward to implement this or something similar into the creation of the DynamicSpectrum.

[joshoewahp]

So there's probably a few avenues to take, depending on what information you'd like dstools to ultimately extract.

Keeping things at the most basic, here's my line of thinking. RM-Synthesis effectively does the de-rotation for you (in an optimal sense) across all samples RMs. So all the information you need is already there after doing it. By taking the peak in the FDF (and assuming the Faraday simple case) we get the band-averaged value in-place without any worries of bandwidth depolarisation.

In general, I'm wary of doing a derotation. Or, at the very least, only presenting values from derotating. I think it'd generally be better to do a rotation after the fact, so that the user has control of how the derotation is done. This would allow for, e.g. a Galactic foreground correction that is constant with time. EDIT: Also, deriving an RM to derotate and subtract will add to the error budget of the residual spectrum. I'm guessing that's why your test had PAs going all over the place

OK that sounds sensible for computing polarised intensity. I'm still a bit confused about getting the appropriate intrinsic PA (i.e. de-rotated to $\lambda^2=0$), as we still need to choose a single RM for this right?

Below is an example comparing a few approaches to computing the polarised intensity / PA.

First following your 2D RM synthesis example—with bottom panel showing L in light blue as the magnitude of the complex FDF cube evaluated at $\phi_{max}$ for each timestep, and middle panel showing PA computed directly from Q/U (black), de-rotating to $\lambda^2=0$ with the RM at Stokes I peak of $961.8,rad/m^2$ (red), and de-rotating each timestep at the max RM(t) shown in the top panel (blue).

Then for comparison, following the approach I described above by first de-rotating with an RM of $961.8,rad/m^2$, then forming Q/U lightcurves by averaging across the band, and forming L (dark blue, bottom panel) and PA (orange, middle panel) from those lightcurves.

We get similar results which is great—but given L determined via RM synthesis incorporates the whole FDF, is there an analogous way to compute the intrinsic PA using more information than the FDF peak?

Some things to consider along the way. It's generally more robust to the fit the peak of the FDF - I have some helper tools in rm_lite for this. The min/max/delta RM should also be selected based on the frequency coverage of a given observation as well to produce optimal results.

OK perhaps fitting the FDF peak in each timestep will provide a more stable derotated PA. Would also be cool to look at fitting the 2D FDF directly to track smoothly varying RM?

[AlecThomson]

@joshoewahp would you be happy top share the data above? I'm keen to have a play with a few things.

Another idea to throw out for consideration is the use of the 1D tools. The top-level 1D functions would do all the nice things for us, including error propagation. The downside is that we'd lose the efficiency of doing the 2D rmsyth in one go. That is, we'd have to call the 1D tool on each time step. This could be parallelised in some form, but would still be less efficient.

The other option would be to build out some 2D top-level tools to do this for us.