Choices for the concentration model

[aesar77]

Hi,

In cluster_lens.ipynb notebook tutorial I get this warning:

ERROR The Eisenstein & Hu 98 transfer function cannot be computed for Ob0 = 0. /lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model diemer19, mass 3.37e+14, mdef vir. warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \

I figured it was caused by the function concent_m_w_scatter from slsim.Halos.halo_population.py when I call slsim.Deflectors.ClusterDeflectors.draw_deflector():

`def concent_m_w_scatter(m, z, sig):
"""Concentration parameter of halos in B. Diemer and A. V. Kravtsov, 2015
with updated parameters of Diemer & Joyce 2019 arXiv:1407.4730 [astro-
ph.CO]. arXiv:1809.07326 [astro-ph.CO].

:param m: halo mass
:type  m: nd.array
:param z: redshift
:type  z: float
:param sig: intrinsic scatter of logarithmic concentration parameter
:type  sig: float
:return: con_halo: ndarray, concentration parameter of halos
"""
con_mean = concentration.concentration(m, "vir", z, model="diemer19")
sca = np.random.lognormal(0.0, sig, len(m))
con = con_mean * sca
con[con < 1.0] = 1.0  # TODO check
return con`

I checked the concentration.py from the Colossus package. I figured that, if the mass halo definition mdef ("vir") is not a native definition of the given model ("diemer19") then it tries to compute c and M numerically, instead of just returning the known model concentration. In conentration.py,

`models = OrderedDict()
"""
Dictionary containing a list of models.

An ordered dictionary containing one :class:ConcentrationModel entry for each model.
"""

models['bullock01'] = ConcentrationModel()
models['bullock01'].mdefs = ['200c']

models['duffy08'] = ConcentrationModel()
models['duffy08'].mdefs = ['200c', 'vir', '200m']

models['klypin11'] = ConcentrationModel()
models['klypin11'].mdefs = ['vir']

models['prada12'] = ConcentrationModel()
models['prada12'].mdefs = ['200c']
models['prada12'].universal = True

models['bhattacharya13'] = ConcentrationModel()
models['bhattacharya13'].mdefs = ['200c', 'vir', '200m']

models['dutton14'] = ConcentrationModel()
models['dutton14'].mdefs = ['200c', 'vir']

models['diemer15_orig'] = ConcentrationModel()
models['diemer15_orig'].mdefs = ['200c']
models['diemer15_orig'].universal = True
models['diemer15_orig'].depends_on_statistic = True

models['diemer15'] = ConcentrationModel()
models['diemer15'].mdefs = ['200c']
models['diemer15'].universal = True
models['diemer15'].depends_on_statistic = True

models['klypin16_m'] = ConcentrationModel()
models['klypin16_m'].mdefs = ['200c', 'vir']

models['klypin16_nu'] = ConcentrationModel()
models['klypin16_nu'].mdefs = ['200c', 'vir']

models['ludlow16'] = ConcentrationModel()
models['ludlow16'].mdefs = ['200c']
models['ludlow16'].universal = False

models['child18'] = ConcentrationModel()
models['child18'].mdefs = ['200c']

models['diemer19'] = ConcentrationModel()
models['diemer19'].mdefs = ['200c']
models['diemer19'].universal = False
models['diemer19'].depends_on_statistic = True

models['ishiyama21'] = ConcentrationModel()
models['ishiyama21'].mdefs = ['500c', '200c', 'vir']
models['ishiyama21'].universal = False`

we see that diemer19 has no virial mass definition; however, there are other models including this definition. So, my question is, is there a particular reason why this combination was chosen in
concent_m_w_scatter function? Can I just change it to "200c"?

[sibirrer]

Thank you for raising the issue. I assign @furcelay who coded up the cluster module. I think what you would need to do is to define an AstroPy cosmology with a Ob >0 which might not be the default one (a bit surprised thought), or there might not be a good translation between AstroPy and Colossus

[aesar77]

FlatLambdaCDM model from the astropy.cosmology model is imported into the notebook, and set as FlatLambdaCDM(H0=70, Om0=0.3). When I checked the Astropy documentation, I see that the default value for Ob0 is 0, unless it is explicitly set in FlatLambdaCDM. I think it might be causing other issues in cluster_deflectors.py

def set_cosmo(self): "Set the cosmology in Colossus to match the astropy.cosmology instance.""" params = dict( flat=(self.cosmo.Ok0 == 0.0), H0=self.cosmo.H0.value, Om0=self.cosmo.Om0, Ode0=self.cosmo.Ode0, Ob0=self.cosmo.Ob0 if self.cosmo.Ob0 is not None else 0.04897, Tcmb0=self.cosmo.Tcmb0.value if self.cosmo.Tcmb0.value > 0 else 2.7255, Neff=self.cosmo.Neff, sigma8=0.8102, ns=0.9660499, ) colossus_cosmo.setCosmology(cosmo_name="halo_cosmo", **params)

Because I think here it is thought that the default value is None and in that case set to 0.04897, but it is not None hence set to 0. In Colossus cosmology models Ob0=0.04897 corresponds to WMAP1-ML model with H0=68.00, Om0=0.3136 while sigma8=0.8102 correspons to Planck-18 model and ns=0.9960 corresponds to WMAP7-only. Tcmb0 and Neff values match between astropy cosmology and colossus cosmology. Now if the reason a custom Colossus cosmology is defined here with these parameters to match astropy cosmology, I am confused about why we set the astropy cosmology with these particular values FlatLambdaCDM(H0=70, Om0=0.3) that does not match any models when combined with the rest of the parameter choices.

[sibirrer]

we mostly use FlatLambdaCDM(H0=70, Om0=0.3) as these are the only relevant parts to get the expansion history to calculate angular diameter distances (that is the only thing really needed in SLSim/lensing) and hence structure formation parameters, such as Ob have not been plugged in. Might be good if you could updated the astropy instance in this particular notebook

[furcelay]

Hi @aesar77 and @sibirrer, thank you for noticing this problem. I'm also surprised that Ob is zero by default. I'll make an issue from this discussion and modify the cosmology to include the right Ob. Thanks!

[furcelay]

Regarding the choice of the concentration model, I used the one from the Hammoks pipeline (at slsim/Halos/halo_population.py) to avoid having multiple definitions of this function. As you pointed out, the diemer19 model has no virial mass definition, so we should change it to '200c' (likely the best option) or change the concentration model. However, we would need to assess how this affects the Hammocks pipeline.I remember that Masamune Oguri, or one of his students, contributed to that part of the code. @sibirrer, is that correct? Should we include them in the discussion?

[sibirrer]

@furcelay, @kta-cosmo coded this part up, so indeed having him chime in is great. Certainly can matter (and is tricky) to have the right mass definitions for these models