Disagreement at large scales between linear and nonlinear correlation_3d

[tilmantroester]

At large scales, the correlation function obtained with correlation_3d from linear and non-linear power spectra does not agree, even though the power spectra seem to be fine. The Pk2D objects also have the same k range and extrapolation settings, so it's unclear to me where this is going wrong.

import numpy as np

import pyccl as ccl
import matplotlib.pyplot as plt


cosmo = ccl.CosmologyVanillaLCDM()

assert np.allclose(
    cosmo.get_linear_power().get_spline_arrays()[1],
    cosmo.get_nonlin_power().get_spline_arrays()[1]
)

r = np.geomspace(5, 5000, 500)

u = r**2.5
plt.semilogx(r, u*ccl.correlation_3d(cosmo, r=r, a=1, p_of_k_a=cosmo.get_nonlin_power()),
             label="nonlin")
plt.semilogx(r, u*ccl.correlation_3d(cosmo, r=r, a=1, p_of_k_a=cosmo.get_linear_power()),
             label="lin")
plt.legend()
plt.xlabel("$r$ [Mpc]")
plt.ylabel(r"$\xi(r)$")
plt.savefig("xi_r.png")

plt.figure()
k = np.geomspace(1e-5, 1, 500)

plt.loglog(k, cosmo.get_nonlin_power()(k, a=1), label="nonlin")
plt.loglog(k, cosmo.get_linear_power()(k, a=1), label="lin")
plt.legend()
plt.xlabel("$k$ [Mpc$^{-1}$]")
plt.ylabel(r"$P(k)$")
plt.savefig("Pk.png")