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.


closed as off-topic by user191954, Frobenius, Qmechanic Sep 10 '18 at 11:44

  • This question does not appear to be about physics within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\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