SPT-scaled version of cib for websky-0.4

[marcelo-alvarez]

websky-0.4 cib power is higher at ell=3000 than predicted by the model constrained by SPT (which is lower than that constrained by Planck at 217 GHz), so it will be useful to add an option for scaling websky-0.4 CIB to match the SPT ell=3000 power.

[zonca]

the easiest implementation, if it is sufficient, is a table of correction factors to be applied to each Websky CIB map. The frequencies are listed at https://github.com/galsci/pysm/pull/129/files#diff-5df48318d41e2d6bb56076c2ac45a0322e1a3841397efb008b61da46ca0a5473R47-R93 If you provide the table of correction factors, I can take care of implementing them in WebSkyCIB.

[sriniraghunathan]

[sriniraghunathan]

scaling factors are: 1.83, 1.75, 1.59, 1.46 for 93, 145, 225, 278 GHz. But I need to redo this for the above bands that Andrea posted. Like I mentioned in the email, using SPT-based scaling for highest frequencies, like ~857 GHz or above, is going to be tricky.

[marcelo-alvarez]

Thanks for providing the scalings. I think we just need to agree on a sensible extrapolation and interpolation to other frequencies.

How about 1+(2.12-1)/(1+exp[(nu-235)/135])?

[zonca]

I think it would be useful to have at least 1 point at low frequency to make sure the extrapolation is reasonable, maybe 30 GHz? [edit] Sorry nevermind, of course it depends on SPT bands.

Anyway @sriniraghunathan please let us know what you think about @marcelo-alvarez 's curve, if it seems fine, I can proceed to implement it.

[sriniraghunathan]

[sriniraghunathan]

So I picked all CIB maps from cori:/global/cscratch1/sd/xzackli/for_nersc_www/cib/ and convert them from MJy/Sr to uK_CMB using the below conversion. The maps do not exactly match what I got yesterday with that interpolation that Andrea mentioned. As a result the scaling factors are now slightly lower. Please see the attached plot with Marcelo's closed form and a slightly modified fit. The fits are clearly off at high frequencies and also at low frequencies which clearly indicates we cannot trust SPT scaling at those bands.

def websky_cib_unit_conversion(nu, nu0 = 100e9):

    x = nu/56.8e9
    t_cmb = 2.7255 #kelvin
    t1 = 1.05e3
    t2 = (np.exp(x) - 1.)**2.
    t3 = np.exp(-x)
    t4 = (nu / nu0)**-4

    return t1 * t2 * t3 * t4

[marcelo-alvarez]

What are the factors?

We should get something into a branch so we can finally do the actual test where we run PySM for whatever frequencies you used, have it produce the CIB with the scaling, mask it appropriately, and get the correct power spectrum. If it doesn't work then we can iterate to see what didn't work.

[sriniraghunathan]

bands = [ 30.0, 47.4, 70.0, 143, 189, 217, 275, 306, 525, 545, 584, 643, 729, 817, 857, 906, 994, 1080] scaling = [1.64652581, 1.71635281, 1.74148928, 1.64030465, 1.5274546, 1.45415951, 1.30351212, 1.22688181, 0.81361881, 0.78696348, 0.73953811, 0.67810142, 0.60727171, 0.55287663, 0.5329301, 0.51182782, 0.48148208, 0.4591423 ]

these are the factors. Not sure why they are slightly lower than before. But maybe it does not matter too much. Edit: Of course, I should also emphasise again that do not trust very low / high frequeency bands here.

[sriniraghunathan]

and @zonca please note that these scaling are for power. For maps, you must take square roots.

[marcelo-alvarez]

scaling factors are: 1.83, 1.75, 1.59, 1.46 for 93, 145, 225, 278 GHz. But I need to redo this for the above bands that Andrea posted.

You provided the above factors before and now have revised values. I can update the fit I did above based on the new values at the SPT bands you are using. The benefit of the form I developed is that it goes through your points and approaches unity at high frequency where websky is already in broad agreement with Planck and Herschel measurements.

Please post the corrected scalings (at the frequencies corresponding to SPT band centers) so I can redo the fit and we can then discuss whether it is acceptable or not to open a PR to merge in the scalings and/or include them in the current open PR. Thanks.

[marcelo-alvarez]

bands = [ 30.0, 47.4, 70.0, 143, 189, 217, 275, 306, 525, 545, 584, 643, 729, 817, 857, 906, 994, 1080]

scaling = [1.64652581, 1.71635281, 1.74148928, 1.64030465, 1.5274546, 1.45415951, 1.30351212, 1.22688181, 0.81361881, 0.78696348, 0.73953811, 0.67810142, 0.60727171, 0.55287663, 0.5329301, 0.51182782, 0.48148208, 0.4591423 ]

these are the factors. Not sure why they are slightly lower than before. But maybe it does not matter too much. Edit: Of course, I should also emphasise again that do not trust very low / high frequeency bands here.

Sorry, I think our comments are crossing paths. I think it would be ideal if you provide the scalings at the frequencies corresponding to the SPT band centers used for their latest model predictions. This is where we believe the scalings the most. Then I can fit a curve through that with correct low and high frequency asymptotics and we can consider if that works.

[marcelo-alvarez]

@zonca if this issue needs to be resolved soon, then adopting the scaling factors @sriniraghunathan posted above makes sense to me, provided it is clear in the documentation that the scalings are not valid at nu > 300 GHz (I am not worried about very low frequencies as the constraints are very weak there).

If it is necessary at some point in the future to generate maps at nu > 300 GHz that are consistent with the scaling at nu < 300 GHz, this can be done by smoothly extrapolating the scaling factor to nu > 300 GHz appropriately.

[sriniraghunathan]

@marcelo-alvarez if we are interested in 93/145/225/278 GHz, then those maps not available in that cori folder. So, there is nothing to update and what I provided yesterday should be good.

[marcelo-alvarez]

@srini I was able to get what we need from the information you provided, thanks.

@zonca we can go with the following functional form (which matches Srini's scalings at the SPT band center frequencies corresponding to the constraints he used):

def power_scaling(nu):
    return 1 + (1.84-1)/(1+np.exp((nu-227)/75))

The map scaling is square root of the power scaling. For your convenience I also provide the map scaling factors:

print('frequencies:', websky_freqs)
map_scalings = np.sqrt(power_scaling(websky_freqs))
print('map scalings:', map_scalings)
frequencies: [18.7000 21.6000 24.5000 27.3000 30.0000 35.9000 41.7000 44.0000 47.4000
 63.9000 67.8000 70.0000 73.7000 79.6000 90.2000 100.0000 111.0000
 129.0000 143.0000 153.0000 164.0000 189.0000 210.0000 217.0000 232.0000
 256.0000 275.0000 294.0000 306.0000 314.0000 340.0000 353.0000 375.0000
 409.0000 467.0000 525.0000 545.0000 584.0000 643.0000 729.0000 817.0000
 857.0000 906.0000 994.0000 1080.0000]
scalings: [1.3382 1.3375 1.3368 1.3361 1.3354 1.3338 1.3321 1.3314 1.3303 1.3245
 1.3229 1.3220 1.3205 1.3179 1.3127 1.3075 1.3010 1.2888 1.2780 1.2696
 1.2596 1.2346 1.2114 1.2033 1.1858 1.1575 1.1358 1.1153 1.1033 1.0957
 1.0735 1.0639 1.0500 1.0335 1.0163 1.0077 1.0059 1.0036 1.0016 1.0005
 1.0002 1.0001 1.0000 1.0000 1.0000]

Here is a plot that shows the agreement at the percent level to what Srini provided at the SPT bands: Please let me know if you need anything else, otherwise I will close this issue. Thanks.

[zonca]

Very good thanks, I'll close the issue when I'll have completed the implementation

[zonca]

@sriniraghunathan @marcelo-alvarez @xzackli should I apply this correction also to the radio galaxies map?

[sriniraghunathan]

@zonca - No this is just for CIB. I have not checked radio galaxies but could check those tomorrow.

[marcelo-alvarez]

@zonca I think it is ok to include the radio galaxy intensity maps as they are. They are abundance matched to available data as described in https://arxiv.org/abs/2110.15357.