The Sunyaev–Zeldovich Effect

Let's make some maps of the thermal Sunyaev-Zeldovich (SZ) effect. First, let's load up the example halo catalog included in this package. These will be automatically downloaded the first time you load them.

using XGPaint, Plots

# example , ra and dec in radians, halo mass in M200c (Msun)
ra, dec, redshift, halo_mass = XGPaint.load_example_halos()
print("Number of halos: ", length(halo_mass))

Number of halos: 2455621

This small catalog is limited to a relatively small patch of the sky. Before we generate some SZ maps, let's take a look at the halo mass distribution.

b = 10.0 .^ (11:0.25:16)
histogram(halo_mass, bins=b, xaxis=(:log10, (1e11, 1e16)),
    yscale=:log10, label="", xlabel="Halo mass (solar masses)", ylabel="counts")

Also notice the steep dropoff in counts at $10^{12} M_{\odot}$, corresponding to the halo resolution of the Websky simulation from which these halos were taken. Fortunately, the SZ effect is dominated by the most massive halos.

To allow for safe threaded painting on a single map, we'll sort the halo catalog by declinations.

ra, dec, redshift, halo_mass = sort_halo_catalog(ra, dec, redshift, halo_mass);

(Float32[0.03305517, -0.052936707, 0.07140066, -0.06599737, -0.05942857, -0.05192846, 0.022953046, 0.01955623, 0.05080256, 0.00839793  …  0.060685053, -0.061119348, -0.040161338, 0.047465798, -0.041020952, -0.048511826, -0.018889412, 0.058788735, -0.050192643, -0.0027171609], Float32[-0.11871151, -0.118711114, -0.118711084, -0.11871107, -0.1187108, -0.118710734, -0.11871068, -0.118710674, -0.11871066, -0.11871049  …  0.10907136, 0.1090714, 0.10907147, 0.109071486, 0.1090715, 0.109071545, 0.10907169, 0.10907171, 0.10907183, 0.10907187], Float32[2.583177, 0.6798465, 1.4767964, 2.4593022, 1.3147808, 1.7026237, 1.313409, 1.2228533, 2.061896, 1.7019778  …  0.4401552, 3.107578, 1.2217613, 2.3513637, 2.2199023, 1.7342392, 0.6818484, 0.8421148, 1.2661853, 2.0941606], Float32[1.8216356f12, 3.4135917f12, 4.958732f12, 3.3517478f12, 2.6270643f12, 2.1461618f12, 1.5084386f12, 4.1571304f12, 1.3568221f12, 2.2222684f12  …  4.438609f12, 1.3447792f12, 1.5955485f12, 1.7920556f12, 8.5669514f12, 2.6231164f12, 2.2672384f12, 2.409821f13, 5.5713213f12, 2.7740707f12])

This package relies on a precomputed interpolation table for performance when making large sky maps of the SZ effect. For this tutorial, we'll load a profile and the associated interpolator from disk.

model, interp = XGPaint.load_precomputed_battaglia()
print(model)

Battaglia16ThermalSZProfile{Float64, Cosmology.FlatLCDM{Float64}}(0.15804878048780488, Cosmology.FlatLCDM{Float64}(0.6774, 0.6924089058605902, 0.3075, 9.109413940981078e-5))

Note

If you want to generate your own model (such as varying cosmology), you would instead configure the appropriate model and then build an interpolator. The interpolator generation is multithreaded and takes about five minutes on 16 cores.

model = Battaglia16ThermalSZProfile(Omega_c=0.2589, Omega_b=0.0486, h=0.6774)
interp = build_interpolator(model, cache_file="cached_b16.jld2", overwrite=true)

This will save the results to the cache_file as well. If you want to load a result from disk, you can specify overwrite=false^[1].

Now, we'll set up the map to generate. We will construct a small CAR (Clenshaw-Curtis variant) patch on the sky. You have to construct a workspace for each new sky patch.

using Pixell
box = [4.5   -4.5;           # RA
       -3     3] * Pixell.degree  # DEC
shape, wcs = geometry(CarClenshawCurtis{Float64}, box, 0.5 * Pixell.arcminute)
m = Enmap(zeros(shape), wcs)

# construct on workspace on this patch
workspace = profileworkspace(shape, wcs)

XGPaint.CarClenshawCurtisProfileWorkspace{Pixell.Enmap{Float64, 2, Matrix{Float64}, Pixell.CarClenshawCurtis{Float64}}}(Enmap(shape=(1080, 720),wcs=CarClenshawCurtis{Float64}(naxis=2,cdelt=(-0.008333333333333333, 0.008333333333333333),crval=(0.0, 0.0),crpix=(541.0, 361.00000000000006))), Enmap(shape=(1080, 720),wcs=CarClenshawCurtis{Float64}(naxis=2,cdelt=(-0.008333333333333333, 0.008333333333333333),crval=(0.0, 0.0),crpix=(541.0, 361.00000000000006))), Enmap(shape=(1080, 720),wcs=CarClenshawCurtis{Float64}(naxis=2,cdelt=(-0.008333333333333333, 0.008333333333333333),crval=(0.0, 0.0),crpix=(541.0, 361.00000000000006))), Enmap(shape=(1080, 720),wcs=CarClenshawCurtis{Float64}(naxis=2,cdelt=(-0.008333333333333333, 0.008333333333333333),crval=(0.0, 0.0),crpix=(541.0, 361.00000000000006))))

Now it's time to apply the model to the map.

@time paint!(m, model, workspace, interp, halo_mass, redshift, ra, dec)
plot(log10.(m), c = :thermal)

Healpix

To generate Healpix maps, you'll need the Healpix.jl package.

using Healpix

nside = 4096
m_hp = HealpixMap{Float64,RingOrder}(nside)
w0 = XGPaint.HealpixProfileWorkspace(nside,
    exp(minimum(interp.ranges[1])), π/180, interp)
ws = [XGPaint.HealpixProfileWorkspace(nside, w0.θmin, w0.θmax,
    interp, w0.posmap) for _ in 1:Threads.nthreads()];


@time XGPaint.paint!(m_hp, model, ws, interp, halo_mass, redshift, ra, dec)
# Healpix.saveToFITS(m_hp, "!y.fits", typechar="D")

plot(m_hp)

1For maps with many pixels and halos, building the interpolator is only about 10% of the map generation cost. If you have a use case where you instead need to vary cosmology on very small patches, contact the developers!