I've being working with one of the 'mini' snapshots for TNG100, and found that NeutralHydrogenAbundance was not a data field for the mini snaps. I just wanted to ask if there was any method or work-around to compute neutral hydrogen abundance with the data supplied in the mini snaps?
Dylan Nelson
30 Jul '21
Yes, this field was removed from mini snapshots since it can be, in theory, reconstructed from the other available fields, following Rahmati+13. I post some code which may be a useful starting point:
def neutral_fraction(nH, sP, temp=1e4, redshift=None):
""" The neutral fraction from Rahmati+ (2012) Eqn. A8. """
# recombination rate from Rahmati+ (2012) Eqn. A3, also Hui & Gnedin (1997). [cm^3 / s] """
lamb = 315614.0/temp
alpha_A = 1.269e-13*lamb**1.503 / (1+(lamb/0.522)**0.47)**1.923
# photoionization rate
if redshift is None:
redshift = sP.redshift
_, gamma_UVB_z = uvbPhotoionAtten(np.log10(nH), np.log10(temp), redshift)
# A6 from Theuns 98
LambdaT = 1.17e-10*temp**0.5*np.exp(-157809.0/temp)/(1+np.sqrt(temp/1e5))
A = alpha_A + LambdaT
B = 2*alpha_A + gamma_UVB_z/nH + LambdaT
return (B - np.sqrt(B**2-4*A*alpha_A))/(2*A)
Ryan Flynn
2 Aug '21
Thank you! One question about the code you sent - what is the function uvbPhotoionAtten? Is that also from the Rahmati paper?
Dylan Nelson
2 Aug '21
Yes gamma_UVB_z is the value of the hydrogen photoionization rate, at the redshift of interest. Table 2 of Rahmati+ (should use FG09 to be most consistent with TNG).
Apologies on re-opening this thread, but I've also been working on reconstructing this field using the particle data.
After replicating the above code and using the FG09 model for gamma_UVB, I've created a plot of temperature vs neutral hydrogen abundance as a check for the most massive subhalo in TNG100-1 at snapshot 99, the left panel being the actual field in the TNG catalog and the right being my reconstruction. Although the shape of the curve looks correct, the two big issues are 1) there is a systematic offset in the y-direction of about 0.3 dex in my approximation, and 2) we aren't resolving the same level of scatter as seen in the original quantity. I'm wondering if I have missed a step or have done a wrong calculation, perhaps in calculating nH?
Dylan Nelson
20 May '22
Perhaps plot vs density instead of vs temperature? Perhaps the discrepancy is arising only for gas with density above the star-formation threshold?
Austen Gabrielpillai
6 Jun '22
Using a smaller subhalo for faster plotting and calculations, I've attached the following plot of the the two scatters overlaid on top of each other for a better, direct comparison. The left being neutral hydrogen abundance vs. temperature, and the right is neutral hydrogen abundance vs. density for both the TNG quantity and our approximation.
The offset between the two does not seem to be correctable by a factor, as taking the ratio between the truth and approximation shows that the ratio varies between 0.6 and 244. It also does not seem like an entirely simple fix as we can see different behavior in the right plot between Density [10^-5, 10^-3] and NeutralHydrogenAbundance [10^-2, 10^0].
Note that my temperature calculations come from the FAQ,
XH = 0.76 # the hydrogen mass fraction
gamma = 5.0/3.0 # the adiabatic index
KB = 1.3807e-16 # the Boltzmann constant in CGS units [cm^2 g s^-2 K^-1]
mp = 1.6726e-24 # the proton mass [g]
mu = (4*mp) / (1 + 3*XH + 4*XH*subhalo['ElectronAbundance'])
temp = (5/3 - 1) * (subhalo['InternalEnergy'] / KB)* mu* 1e10
my nH calculations are based on the following,
# density convert to Msun / kpc^-3 Msun to g kpc^-3 to cm^-3 divide by proton mass (g) to get number density
nH = 0.76 * np.double(subhalo['Density']) * (1.0e10 / 0.6774) * 0.6774**3 * 1.989e33 * 3.086e21**-3 / (mp)
and the neutral_fraction from the code snippet above in the thread.
Reviewing Rahmati+13. are we supposed to be using the total photoionization rate (given by A1) and not just the UVB ionization rate, or does the UVB give a good enough approximation?
Dylan Nelson
6 Jun '22
Perhaps do a temperature vs density plot to check? If by "Density" on the x-axis you aim to mean "Hydrogen Number Density [atoms / cm^3]", then something doesn't seem right. The effective ISM model should kick in at ~0.1 cm^-3, but this looks to be the feature at 1e-3 in your plot.
Austen Gabrielpillai
8 Jun '22
Apologies, by 'Density' I am referring to the direct quantity from the particle data. Attached is a new plot, where the first panel is temperature vs. the catalog's density parameter (in units of (10^10 M⊙/ h) / (ckpc / h)^3 ) (left) and temperature vs. hydrogen number density (in units of atoms / cm^3). The effective ISM model feature that you mentioned does appear to occur in the n_H vs. temperature plot at 0.1 cm^-3.
I've being working with one of the 'mini' snapshots for TNG100, and found that NeutralHydrogenAbundance was not a data field for the mini snaps. I just wanted to ask if there was any method or work-around to compute neutral hydrogen abundance with the data supplied in the mini snaps?
Yes, this field was removed from mini snapshots since it can be, in theory, reconstructed from the other available fields, following Rahmati+13. I post some code which may be a useful starting point:
Thank you! One question about the code you sent - what is the function uvbPhotoionAtten? Is that also from the Rahmati paper?
Yes
gamma_UVB_zis the value of the hydrogen photoionization rate, at the redshift of interest. Table 2 of Rahmati+ (should use FG09 to be most consistent with TNG).Apologies on re-opening this thread, but I've also been working on reconstructing this field using the particle data.
After replicating the above code and using the FG09 model for gamma_UVB, I've created a plot of temperature vs neutral hydrogen abundance as a check for the most massive subhalo in TNG100-1 at snapshot 99, the left panel being the actual field in the TNG catalog and the right being my reconstruction. Although the shape of the curve looks correct, the two big issues are 1) there is a systematic offset in the y-direction of about 0.3 dex in my approximation, and 2) we aren't resolving the same level of scatter as seen in the original quantity. I'm wondering if I have missed a step or have done a wrong calculation, perhaps in calculating nH?
Perhaps plot vs density instead of vs temperature? Perhaps the discrepancy is arising only for gas with density above the star-formation threshold?
Using a smaller subhalo for faster plotting and calculations, I've attached the following plot of the the two scatters overlaid on top of each other for a better, direct comparison. The left being neutral hydrogen abundance vs. temperature, and the right is neutral hydrogen abundance vs. density for both the TNG quantity and our approximation.
The offset between the two does not seem to be correctable by a factor, as taking the ratio between the truth and approximation shows that the ratio varies between 0.6 and 244. It also does not seem like an entirely simple fix as we can see different behavior in the right plot between Density [10^-5, 10^-3] and NeutralHydrogenAbundance [10^-2, 10^0].
Note that my temperature calculations come from the FAQ,
my nH calculations are based on the following,
and the neutral_fraction from the code snippet above in the thread.
Reviewing Rahmati+13. are we supposed to be using the total photoionization rate (given by A1) and not just the UVB ionization rate, or does the UVB give a good enough approximation?
Perhaps do a temperature vs density plot to check? If by "Density" on the x-axis you aim to mean "Hydrogen Number Density [atoms / cm^3]", then something doesn't seem right. The effective ISM model should kick in at ~0.1 cm^-3, but this looks to be the feature at 1e-3 in your plot.
Apologies, by 'Density' I am referring to the direct quantity from the particle data. Attached is a new plot, where the first panel is temperature vs. the catalog's density parameter (in units of (10^10 M⊙/ h) / (ckpc / h)^3 ) (left) and temperature vs. hydrogen number density (in units of atoms / cm^3). The effective ISM model feature that you mentioned does appear to occur in the n_H vs. temperature plot at 0.1 cm^-3.