galsci / pysm

PySM 3: Sky emission simulations for Cosmic Microwave Background experiments
https://pysm3.readthedocs.io/
BSD 3-Clause "New" or "Revised" License
36 stars 23 forks source link

Synchrotron templates with small scales in polarization fraction tensor formalism #73

Closed b-thorne closed 3 years ago

b-thorne commented 3 years ago

Simple synchrotron small-scales

Start implementing synchrotron polarization fraction tensor approach, as described in #72.

Plan is to use Haslam 408 MHz template for intensity emission, and WMAP K-band for Q/U. Fit power law for high-ell slope to power spectra of polarization tensor quantities (logs of maps), using same range of ell as in the original PySM paper.

Work is being done in this script.

b-thorne commented 3 years ago

First a sanity check of power spectra and fits for index using Stokes params. Compare to original PySM paper figure 5.

Here we're using the same combination of masks as #72, rather than the masks described in the PySM paper, and we're using namaste rather than spice.

fitted.pdf

Screen Shot 2021-02-21 at 11 47 58 AM

Temperature

The signal-dominated part of the TT spectrum seems to have ~three regimes with differing slopes. We'll see if this is also true of the transformed quantities, but perhaps will have to choose a range to fit for gamma.

Polarization

There is a bit of a discrepancy between the polarized spectra found here and in the PySM paper. This could be due to the differences of mask / power spectrum estimation. The fitted indices come from a least squares fit over the multipole range 2 to 36, without a consideration of the errors, unlike in PySM. We get a shallower slope here, but perhaps would get a different result if we used sample variance errors like was done in the paper.

zonca commented 3 years ago

thanks @bthorne93, a notebook would be better because then I can turn into HTML and add all of it (code, explanation and figures) to the docs. But if you don't like to use notebooks, I can transform this in a notebook once it is finalized.

zonca commented 3 years ago

I am working on a new version of my dust fitting notebook that goes through the whole process, I will have it ready pretty soon.

b-thorne commented 3 years ago

If we apply the transform naively to the maps at their native resolutions, we get the following:

transformed_q_naive.pdf

The noise in the WMAP K band is large, and at higher latitudes, this means for a lot of pixels P is larger than I, and the transformation is no longer valid, as we will be taking logs of negative numbers. To remove this effect, it makes sense to smooth both maps to a common resolution, here we choose 2degrees, which is roughly the same as was done in the PySM paper.

Applying the transformation to these smoothed maps, one gets:

transformed_q_2degfwhm.pdf

This looks more reasonable, however we are still getting localized regions where the transformation fails. These anomalous regions are mostly close to the Galactic plane, and presumably represent regions of intense emission, as they pop out in both the Haslam and WMAP maps. To proceed, I have taken the remaining pixels for which P>I, and manually scaled them down to P=0.99 I. This ad hoc method only affects a small number of pixels in localized regions of sky, but perhaps we can think of something better? The result looks like:

transformed_q_2degfwhm_corrected.pdf

b-thorne commented 3 years ago

Now we want to extend to small scales by fitting power laws to the power spectra of the transformed variables, which we can refer to as lower case spectra: tt, ee, bb. We calculate power spectra on sky region defined by the joint Planck intensity and power spectrum masks (since the transform mixes intensity and polarization):

transformed_q_2degfwhm_corrected_masked.pdf

The resulting power spectra, and the power law fits on a the interval ell=10-36 are:

fitted_eebb.pdf

seclark commented 3 years ago

On the pixels where P > I: it looks like these regions mostly get masked anyway for the power spectrum fit, but another option is to mask the problematic regions before smoothing and then use the smoothing step to interpolate over them. e.g. smooth the mask, smooth the map with the masked areas set to 0, and then divide the smoothed map by the smoothed mask. Somewhat better justified than just forcing P/I = 0.99, and should work fine since these are localized regions that have angular extents smaller than the smoothing scale.

seclark commented 3 years ago

By the way at least some of those regions where the transformation fails are background radio galaxies. I think I spot Cen-A as that funny double-lobed spot that survives the 2deg smoothing.

b-thorne commented 3 years ago

Once we have fitted the power spectrum, it is not totally clear to me how to transform back to real IQU.

If we follow the same sort of prescription as PySM then we would add a high-pass filtered small scale realization to the low-pass filtered logged quantities. However, since we had to apply the low-pass filter in real space in order to deal with the noise, how should we define the high-pass filter?

I think we could go about this in a slightly different way which I've outlined here, however, that relies on an ongoing piece of development to the CMBLensing.jl code, which is not quite ready yet.

https://www.notion.so/Polarization-tensor-transforms-ad6975b15e5449e29d775d5f52ff818c

zonca commented 3 years ago

in the dust model in #72 I apply another low-pass filter quite sharp using a sigmoid directly to the pol tensor fraction maps. then use the complement to shape the simulated small scales power spectrum.

NicolettaK commented 3 years ago

Ciao, as discussed a the call this week, I have uploaded a notebook with my proposal on how to modify the PySM template for synchrotron spectral index taking into account S-PASS information. It is here: https://gist.github.com/NicolettaK/f801f299a7de116651758ea286eaed66 Let me know if you have any comment/questions!

giuspugl commented 3 years ago

Thanks @NicolettaK ! few questions :

NicolettaK commented 3 years ago

@giuspugl:

giuspugl commented 3 years ago

@zonca how would you foresee to proceed for synchrotron ? Shall we open a new PR with the implementation of spatially varying beta_s proposed by @NicolettaK ?

NicolettaK commented 3 years ago

There are a couple of things we need to decide:

brandonshensley commented 3 years ago

There are a couple of things we need to decide:

  • Which average value of beta_s do we want? Current PySM template is at -3, but Planck and low frequency observations point towards steeper values

Agreed that this needs to be re-evaluated. Do you have a recommendation @NicolettaK?

  • which angular resolution do we want in the new template? Original PySM is at 5° but we are adding gaussian small scales and we can do it up to whatever ell

Is this in reference to the amplitude map or the beta map? Agree that we should make the cutoff at the smallest scale we think we trust and then add a random realization for the smaller scales for both maps.

  • Shall we go with only one realization or we want multiple ones?

The idea should be that all new models are set up for generating random realizations of the small scales on the fly, including spectral parameters. We can work in analogy to d10 and d11, where the first is just one specific realization of the second.

zonca commented 3 years ago

@zonca how would you foresee to proceed for synchrotron ? Shall we open a new PR with the implementation of spatially varying beta_s proposed by @NicolettaK ?

yes, let's make a dedicated notebook in a new Pull Request that downloads all inputs on the fly like #82 (and stores them in data/ subfolder) and produces outputs alm and maps (in data/). We will discuss later how to organize the Synchrotron models.

review-notebook-app[bot] commented 3 years ago

I've committed what I have so far on the synchrotron templates. I showed the polarization part on the call on Wednesday, and have tried to extend this to intensity. However, I was finding it harder to get the intensity part working.

Check out this pull request on  ReviewNB

See visual diffs & provide feedback on Jupyter Notebooks.


Powered by ReviewNB

NicolettaK commented 3 years ago

@brandonshensley for the average value of beta_s I would go for -3.1, I think this is the best fit value from Planck+WMAP.

For the angular resolution, I was referring to the beta map, up to which scales do we want to have it?

giuspugl commented 3 years ago

We decided yesterday to inject smaller angular scales up to the ell =1000 , according to the plot below image @NicolettaK can we inject smaller scales to the Beta_s map up to this ell ? What do you think ?

NicolettaK commented 3 years ago

Hi, we can decide to inject gaussian small scales on the beta_s map up to any ell we decide. Although I'm not fully convinced we should go below the degree scale, we really don't have any observational data at those scales for beta_s that confirms that there are indeed variations.

giuspugl commented 3 years ago

I went through the notebook posted by @bthorne93. I think i found the reason of the inconsistency observed in the power spectra: image

Finally, I consider two different evaluation of Namaster spectra.

Results w/ Joint Mask

a 2deg apodized synchrotron mask (referred as JointMask ) masking point sources and Gal. Spur, and the HFI mask GAL080 as the ones used in #82. image image

The resulting power spectra are very much affected by the point sources at high Gal. latitudes, making less clear the how smaller angular scales are injected for the intensity. Notice also the increase at higher ell due to the noise in the input maps.

Results w/ GAL080

image

Vice versa we can easily infer with Gal080 how the small scales are injected, being a more regularized mask involving only the Galactic midplane. the polarization fraction as estimated from the input maps is highly contaminated by noise i therefore lowpass filtered the map so that the smaller angular scales (mainly due to noise are excluded) and we got similar distributions with the output maps with small angular scales: image

Finally, i checked the E-to-B ratio and the case with small scales injected shows a ratio < 1.

image

/!\ Warnings

brandonshensley commented 3 years ago

The power spectra of the map we recover after masking point sources looks great. I am not following the argument as to why Haslam is preferable to WMAP for the intensity map? This paper on the NCP is one piece of evidence that scaling Haslam would be an ok thing to do if we need to. Agree about the need to carefully document whatever procedure we settle on.

giuspugl commented 3 years ago

I am not following the argument as to why Haslam is preferable to WMAP for the intensity map? @bthorne93 would you mind to motivate this choice ?

Screen Shot 2021-10-01 at 9 09 42 AM

the two maps even morphologically look different ...

giuspugl commented 3 years ago

Although I'm not fully convinced we should go below the degree scale, we really don't have any observational data at those scales for beta_s that confirms that there are indeed variations.

I see, so perhaps we can produce a beta_s map up to the degree scale as @NicolettaK prescribes. For smaller scales in the synchrotron maps, we will just have the scales injected in the LogPol tens maps. How does that sound ?

brandonshensley commented 3 years ago

Although I'm not fully convinced we should go below the degree scale, we really don't have any observational data at those scales for beta_s that confirms that there are indeed variations.

I see, so perhaps we can produce a beta_s map up to the degree scale as @NicolettaK prescribes. For smaller scales in the synchrotron maps, we will just have the scales injected in the LogPol tens maps. How does that sound ?

I think we should have small scales in the beta maps too--it is more likely that they are there than that they aren't, and we have a reasonable estimate of the power spectrum of the beta_s map. It sounds right to place the transition between data and Gaussian realization at the degree scale.

brandonshensley commented 3 years ago

I am not following the argument as to why Haslam is preferable to WMAP for the intensity map? @bthorne93 would you mind to motivate this choice ?

Screen Shot 2021-10-01 at 9 09 42 AM

the two maps even morphologically look different ...

I see, WMAP 30 GHz probably has significant AME in it, too. The scaled Haslam map definitely looks better by eye, so agreed that it makes the most sense to just use that.

NicolettaK commented 3 years ago

The transition to gaussian for beta_s will be at 5 deg, as this is the angular resolution of the only beta_s full sky map that we have (the one currently in PySM)

b-thorne commented 3 years ago

Apologies, just taking a look at what @giuspugl has posted.

The main reason for using Haslam, as opposed to the WMAP Kband or Planck component-separated maps, is due to the high amount of AME and free-free present in the latter two. Since the log polarization tensor transformation mixes intensity and polarization, I also had to make an assumption about how to scale the Haslam template to 23 GHz, for which I chose a spatially-constant index of -3, but which could be a different choice made to be consistent with our modeling of beta_s.

In the notebook, I compared using two masks: GAL20 (which I should really have called GAL80 to be consistent with @giuspugl) and the JointMask, which are quite different:

Screen Shot 2021-10-04 at 2 37 34 PM

I found that it was difficult to get well-modulated intensity small-scales using the Haslam template with either of these masks. @giuspugl, were you able to make this work with the recaled Haslam template and either mask? It looks like your power spectra above are using the WMAP Kband intensity map.

giuspugl commented 3 years ago

I found that it was difficult to get well-modulated intensity small-scales using the Haslam template with either of these masks. @giuspugl, were you able to make this work with the recaled Haslam template and either mask? It looks like your power spectra above are using the WMAP Kband intensity map.

No, I used the Kband intensity map for sake of comparison, but all the analysis shown above is done by assuming the HASLAM template and by adopting the same preprocessing and rescaling as @bthorne93 did. see the notebook in https://gist.github.com/giuspugl/c44d061bde1bf1391c3d10ae29acb1e0

zonca commented 3 years ago

continues in #98