Recently I have been trying to teach myself the basics of modern cosmology. To do so I have been using Barkana & Loeb.

I've been working my way through and come to virial radius. They give the expression;


Where Where $\Delta_{c}=18\pi^{2}+82(\Omega_{m}^{z}-1)-39(\Omega_{m}^{z}-1)^{2}$, and $\Omega_{m}^{z}=\frac{\Omega_{m}(1+z)^{3}}{\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}+\Omega_{k}(1+z^{2})}$.

This is fair enough, and I believe I understand both the expression and the derivation. They then go on to plot the virial radius against redshift for various masses;

enter image description here

To check my understanding, I tried to do the same, but my plot, despite having a similar scale, is nowhere near steep enough and doesn't intercept the $r_{vir}=0.01$ line in the given range at all;

enter image description here

My plot was generated via the following python code;

import numpy
import matplotlib.pyplot as plt

k=1 # for closed overdensity, by Birkhoff's theorem
omega_k = -k/(H_0**2)

def omega_mz(z):

    num = omega_m*(1+z)**3
    den = num+omega_l+omega_k*(1+z)**2

    return num/den

def Delta_c(z):

    return 18*numpy.pi**2+82*(omega_mz(z)-1)-39*(omega_mz(z)-1)**2

def r_vir(M, z):

    factor_1 = M/(1E8*(1/h))
    factor_2 = (omega_m/omega_mz(z))*(Delta_c(z)/(18*numpy.pi**2))
    factor_3 = (1+z)/10

    return 0.784*(factor_1**(1/3))*(factor_2**(-1/3))*(factor_3**(-1))/h

(fig, ax) = plt.subplots()
ax.semilogy(z_vals, r_vir(1E4, z_vals), label="M=1E4")
ax.semilogy(z_vals, r_vir(1E8, z_vals), label="M=1E8")
ax.semilogy(z_vals, r_vir(1E14, z_vals), label="M=1E14")
plt.xlim(0, 30)
ax.set(xlabel='z', ylabel='r_vir');

I've been fiddling around with this code and the scale of the data for a few hours now and just can't seem to get the two plots to agree. I feel like I must be missing something really conceptually simple (as I said, I'm new to this), and would really appreciate some help figuring it out.

p.s. I wasn't sure whether this would do better here or somewhere like the numerics SE, so feel free to migrate if you think it would be more appropriate elsewhere.

EDIT: At the suggestion of Chair, I've computed the above manually for some specific values to help identify whether the error is numerical, conceptual or with the formula itself. I computed $r_{vir}(z=0,M=10^{12}M_{\odot})\approx 2\times 10^{2}kpc$, and $r_{vir}(z=30,M=10^{12}M_{\odot})\approx 6.9kpc$. The value at $z=0$ is similar to the Barkana & Loeb plot, but the $z=30$ value is far too large. This reflects what happened with my numerically generated plot, so I'm pretty sure the problem isn't with my code. Instead I think the error must either be conceptual or with the formula itself.

  • $\begingroup$ You mentioned that you may be missing something conceptual, in which case your quesiton would be on-topic here. But questions about debugging code are off-topic. Could you test if the problem is with the equations or the code? maybe try throwing some sample values at the equations (manually, without a loop) and see if they predict the correct values. If you find out that the equations aren't right, you could ask a question here about those, and if the equations are correct, try comp science SE or stack overflow (I don't know which of them would be appropriate; look at their help centers). $\endgroup$ – user191954 Sep 10 '18 at 5:53
  • $\begingroup$ Thanks for the advice, I plugged in some numbers and believe its a conceptual error or an error with the formula rather than a numerical error (see my update for details). $\endgroup$ – Confused_Cosmologist Sep 10 '18 at 7:19
  • $\begingroup$ Looks to me that the associated caption to the B&L plot says they're plotting $n-\sigma$ fluctuations; while not reading the rest of the paper, I'd think the equation you have isn't what is plotted in the figure you reference. $\endgroup$ – Kyle Kanos Sep 10 '18 at 12:16
  • $\begingroup$ I though by $n-\sigma$ fluctuations they meant they were plotting the virial radius against redshift for masses $n-\sigma$ from the mean halo mass? In this case both plots would be of the same quantity. I used a range of masses which should be similar to the masses corresponding to the $1-\sigma$, $2-\sigma$, and $3-\sigma$ fluctuations they plotted in my plot. $\endgroup$ – Confused_Cosmologist Sep 10 '18 at 12:42
  • $\begingroup$ @Confused_Cosmologist like I said, I'm basing my remark off the caption & not the text. $\endgroup$ – Kyle Kanos Sep 10 '18 at 15:14