Higher Order Multipoles Suppressed for ConvolvedFFTPower

[isands17]

Describe the bug I am using ConvolvedFFTPower to estimate the redshift-space power spectrum for a lognormal galaxy catalog with uniform randoms, and the BOSS survey window and redshift cuts applied. I expect the resulting power spectrum multipoles to roughly match the Kaiser predictions for even multipoles; the monopole is pretty close, but the quadrupole and hexadecapole are much smaller than expected. Below, the solid lines are the results from ConvolvedFFTPower, and the dashed lines are the Kaiser predictions:

Steps To Reproduce

1. Create data and randoms catalogs: data = LogNormalCatalog(Plin=Plin, nbar=3e-4, BoxSize=3e3, Nmesh=256, bias=b1, seed=49) randoms = UniformCatalog(nbar=3e-3, BoxSize=3e3, seed=98)

2. Add RSD: data['Position'] = data['Position'] + data['VelocityOffset'] * data['LOS']. The LOS is a unit vector between each object in the data catalog and an "observer" at a specified location (in this case, the origin of the box).

3. Transform Cartesian coords to sky coords (also requires an observer location; here, I use the center of the box to make sure the entire BOSS window is contained in the simulation).

4. Apply BOSS survey window (using healpy mask) and cut redshift to range 0.2 < z < 0.5 for both data and randoms.

5. Create redshift histogram and calculate n(z)

6. Create fkp catalog and convert to mesh, then call ConvolvedFFTPower: mesh = fkp.to_mesh(Nmesh=256, compensated=True, BoxSize=3000., resampler='tsc',nbar='NZ', fkp_weight='WEIGHT_FKP', comp_weight='Weight') r = ConvolvedFFTPower(mesh, poles=[0,1,2,3,4], kmin = 0.)

Actual behaviour The resulting power spectrum multipoles are suppressed for l = 2, 4.

Expected behaviour Power spectrum multipoles l = 2, 4 should be closer to Kaiser predictions.

Screenshots If applicable, add screenshots to illustrate your problem. For example, this can be a screenshot of the relevant sections of your notebook.

Version

0.3.14

Additional context

This issue is independent of mesh size; the plot above is from a run with nMesh = 256, but is nearly identical to a version I ran on a cluster with nMesh = 1024. It also seems to be independent of redshift cuts and the survey window applied, as long as n(z) is close to 3e-4 (the nbar used to generate lognormal catalog). Also, if I shift all of the objects so that they're far from the origin of the box (e.g. data['Position'] += [0,0,10000]), I recover the results from FFTPower for a fixed LOS [0,0,1]. Below, the solid lines are poles from ConvolvedFFTPower, dotted lines are from FFTPower, and dashed are Kaiser predictions:

I think this issue could be related to how the line-of-sight is defined, specifically where the observer is placed. The observer is used in three places:

1. Adding RSD (where LOS is between each object and the origin of the box, (0,0,0)).
2. Transforming Cartesian coordinates to sky coordinates (here, the observer is at (1500, 1500, 1500)).
3. ConvolvedFFTPower itself.

However, it's not clear to me where ConvolvedFFTPower places the observer, or how the LOS is defined in the algorithm! Is it possible to specify the observer location for ConvolvedFFTPower? Or is the higher-order multipole suppression being caused by something else entirely?

Thank you!

[rainwoodman]

Thanks for the analysis.

I believe that ConvolvedFFTPower assumes the observer is at 0, 0, 0 -- actually, if the input is in RA, DEC, them the observer must be sitting at 0, 0, 0. Internally, I think the location of observer is implied by how the Ylm factors are defined (but I may be wrong). If the observer is not at 0, 0, 0, the factors will likely get very complicated.

I noticed that your step 2 and 1 are using inconsistent observer locations -- should the sky coordinate center on the intersection of all of the LOS vectors? Could you try moving the observer at 0, 0, 0 in step 2?

[isands17]

If I move the observer to (0,0,0) in step 2, all of the poles are overestimated:

The issue with doing this is that the Cartesian positions get transfomed to Ra, Dec that don't cover the full range of the survey window function (after selecting for the window, we're left with only ~150,000 galaxies in the data catalog!)

The easiest fix seems to be shifting the Cartesian coordinates by -[1500, 1500, 1500] so that the full survey window fits inside the simulation box (so Cartesian coordinates now range from (-1500, 1500) in each dimension). This gives us poles that seem somewhat closer:

However, I'm not convinced this is related to the LOS used to transform Cartesian coords to sky coords at all. I did a test where I didn't apply the window function and didn't apply redshift cuts, and didn't transform the coordinates. It's just a log-normal data catalog and a uniform randoms catalog going into ConvolvedFFTPower with n(z) = 3e-4, and we see very similar behavior in the higher-order multipoles:

If the box coordinates are from (-1500, 1500) instead of (0,3000), does ConvolvedFFTPower still put the observer at (0,0,0)?

[rainwoodman]

I think the observer is always at 0, 0, 0. Thinking about this a bit more, ConvolvedFFTPower doesn't care about the box coordinate -- it only takes the RA DEC and Z -- the observer sits at Z=0 which is usually 0, 0, 0.

The direction of the RSD LOS matters because it introduces an anisotropy to the density field. If the observer is not co-centered with the RSD LOS, then that anisotropy is wiped out.

Using just the LogNormal without RSD, (I assume that's your figure 3), then there is no anisotropy to begin with -- hence the power at L=2 is 0.

When you do the shifting (Figure 2), did you make sure the randoms and the data cover the same footprint? The fact that you are getting some L=4 seems to suggest that they are misaligned (coming from the rectangular shape?) A quick way to check this is to do 2D histogram in RA, DEC, and then a 1D histogram in Z.

[isands17]

Sorry, I should have clarified– figure 3 does have RSD, it just doesn't have a window function, but there should be anisotropy. The window function just selects for RA, DEC within a certain range, and cutting for redshift restricts the range to 0.2 < z < 0.5 but I don't think that would remove the anisotropy?

Here are my histograms for LOS defined with the observer at (0,0,0) (and with RSD). It seems like the Cartesian to Sky coordinate transformation ends up measuring the redshift as much larger than the z = 0.38 used to generate the linear power spectrum (the distribution is centered around ~ z = 2, so when we cut for 0.2 < z < 0.5, not many galaxies are left). However, the randoms and data redshift distributions match (scaled by a factor of 10, which I expect since the randoms catalog is 10 times larger):

Also, the 2d RA/DEC histogram has the same overall shape, but different densities within the window, which could potentially be a problem (although there are 10 times more randoms than data):

However, if I set the observer to (1500, 1500, 1500) (for both RSD generation and transforming the coordinates), the RA/DEC histograms cover the entire survey window, although the density is similar to before. I'm still seeing the same issue with the redshift histogram as well:

[isands17]

We've tried to simplify the problem to see if we can track it down; we removed all window and redshift selection functions, leaving just a lognormal catalog, uniform randoms, and RSD. Then, we add [0,0,1e5] to the position coordinates to recover the fixed LOS limit, with LOS = [0,0,1]. The resulting multipoles seem reasonable: solid lines are measured multipoles, dashed lines are Kaiser predictions, and solid red is the linear power spectrum.

But if we flip the setup, so that we subtract [0,0,1e5] from all of the coordinates, and change the LOS to [0,0,-1], the poles are suppressed or even negative:

However, if we subtract [0,0,1e5] but use the LOS pointing in the opposite direction [0,0,1] we get the same power spectrum from adding [0,0,1e5]: This seems odd to me– I don't see why coordinates shifted in the negative direction would take a positive LOS.

A screenshot of the relevant code in my Jupyter notebook is below. This is for the fixed-LOS test in the negative direction.

Thanks so much!

[hsgg]

Thank you @rainwoodman for your previous responses! The first two plots in @isands17 's last post clearly show the problem we are having. I can reproduce them locally on my computer as well. If the simulation box is in the +z direction, then we get good agreement between the simulation and flat-sky Kaiser prediction, whereas if the box is in the -z direction, we do not. It is my understanding that ConvolvedFFTPower uses the Yamamoto estimator that selects the correct LOS in either case, and we should get agreement with the flat-sky Kaiser prediction in both cases. What are we missing? Any help appreciated.

[rainwoodman]

@eelregit knows the math of the estimator better and probably can offer some insight into the strange asymmetry between +z and -z.

Thanks for the investigation. The flat sky limit behaviour respect to LOS is puzzling. How do we picture a negative l=2 pole? Does it correspond to flipping the sign of the RSD offset?

I am also curious what happens with LOS=[0, 0, -1], but your global shift remains [0, 0, 1e5]. It could be as simple as an extra |abs| around a dot product that it shouldn't be. It is also encouraging that at least there is one limit it is working correctly.

Also, the code snippet is illustrative, but could you also post the full source code (RA DEC conversion + ConvFFTPower) -- such that we reproduce the problem exactly as your set up -- and we'll eventually set up a unit test based on your case, after we figure this all out.

[isands17]

When LOS = [0,0,-1] and the global shift is [0,0,1e5], I get a similar multipoles to the ones where LOS = [0,0,-1] and the shift was [0,0,-1e5]. Apparently the direction of the LOS is important, not the direction of the shift.

Here's the full code snippet. I actually never convert to RA and DEC, since I do not apply a window function in this example.

[eelregit]

For wide angle case you use the unit vector pointing to the survey, so negative shift and positive LOS won't happen together. Likewise in flat sky limit, it shouldn't happen either.

Mathematically the equation should be [(v_vec dot r_hat) r_hat], so @rainwoodman we should update the documentation about. The code should be fine.

[eelregit]

I have edited to correct myself in the above comments many times. Sorry if you are confused by the multiple notifications.

[hsgg]

Mathematically the equation should be [(v_vec dot r_hat) r_hat], so @rainwoodman we should update the documentation about. The code should be fine.

Thank you @eelregit , I believe that indeed is the problem! Can you confirm @isands17 ?

Then it is indeed a documentation error, e.g., the second code box at https://nbodykit.readthedocs.io/en/latest/cookbook/lognormal-mocks.html#Redshift-space-Power-Spectrum would need to be updated.

[isands17]

Thank you @eelregit! That did fix the problem in the fixed-LOS limit. However, we are still seeing suppressed higher-order multipoles in the general/ wide angle cases– I'll keep looking into it.

[eelregit]

Have you checked the bias? The bias parameter is not the linear bias, I think, but a coefficient c in exp(c \times delta_lin). So maybe you can adjust your Kaiser b1 to match P_0, and then see if P_2 and P_4 work as well.

[isands17]

I think I tracked it down:

data['VelocityOffset']*data['LOS']

produces a 3d vector– it isn't a proper dot product. So

(data['VelocityOffset']*data['LOS'])*data['LOS']

seemed to fix the issue with the fixed-LOS, but not the general case.

With the proper math implemented, the multipole suppression seems to be fixed:

velocityOffset= np.asarray(data['VelocityOffset'])
LOS = np.asarray(data['LOS'])
product = (velocityOffset*LOS).sum(1)
data['Position'] = data['Position'] + LOS*product[:,np.newaxis]

Thank you all for your help!