General framework for baryonic effects with the halo model

[damonge]

It would be really good to have a general, flexible, and efficient implementation of halo model power spectra with baryonic effects. Most implementations split the total matter profile into several components (e.g. modified dark matter profile, bound gas, ejected gas, stars). The total matter profile is then the sum of these. Then, there are many different parametrisations for each of these components, with some implementations sharing the same functional forms, and other combining different ones. Our implementation should be able to play well with this modelling richness. I.e. we should provide a general framework in which users can easily use existing models, modify them, create new components and combine them with existing ones, etc. E.g. a user may like the gas profile from Arico et al. but would like to combine it with a new stellar profile using their own parametrisation. E.g. 2: a user may have a model in which the gas and star parameters share some physical parameters (e.g. to fix their normalisation).

Here is a bare bones proposal for what this could look like.

1. General structure: we create a HaloProfile class called BaryonsCombiner, or something like that, that takes in a list of profile objects and whose _real and _fourier methods consist of adding the result of calling these methods for each profile in the list. In addition to this combiner class, we will provide specific implementations for the ingredients of some of the most popular models, as well as examples for how to combine them.
2. Parameter consistency: the combiner class should be able to handle the free parameters of all profiles in the list gracefully. I.e. calls to its update_parameters method should update the parameters of all profiles in its list, passing only the relevant parameters, including all shared parameters.
3. Normalisation: this is potentially a big deal. The question is: what should the total mass of a halo be in the presence of baryonic effects? I think we should provide the following three options when creating a combined halo profile:
   • Simple/naive normalisation: we enforce that the total mass of the halo (i.e. the volume integral of the density profile, out to infinity) should be equal to the "mass" of the halo, with this mass given by e.g. M200c (or whatever other mass definition is being used). This is ultimately incorrect: M200c will only contain the mass out to the "virial" radius, whereas feedback blows some of it to much larger distances.
   • Strict normalisation: we enforce that the total mass should be equal to the total mass of a "dark-matter-only" profile when integrated out to infinity. I.e. baryonic effects push matter around but don't create/destroy it. This "dark-matter-only" profile would be an additional profile passed to the BaryonsCombiner object that represents the matter density in a gravity-only simulation.
   • None: in this case users are in charge of ensuring that all profiles in their list are consistently normalised. I.e. the _real method of the combined profile will just be the sum of the _reals of all its children profiles, with no additional normalisation.

It would be really good to hear people's thoughts/suggestions about the above. This is only a proposal and I'm sure others have better/alternative ideas.

A few more thoughts:

• I really like @DhayaaAnbajagane 's BaryonForge, and I think we should use several chunks of his codebase here. E.g. I'd like to implement profile operations in CCL natively! Also happy to hear your thoughts Dhayaa.
• We should keep track of how slow our implementation is. I think our first version of it will probably be quite slow, and that's fine. Once we're happy with it, we should work on ways to speed it up (not crucial for the v3 paper, mind you). A usual bottleneck is the FFTLog transforms, so coming up with alternatives would be great. Any ideas based on people's experience would be very welcome!

I'll tag a few other people who I think will be interested in this @anicola @elisachisari @nikosarcevic @rreischke (and I'm sure I'm missing others).

[tilmantroester]

Excellent proposal! One point to potentially keep in mind are some of the inherent limitations of (vanilla) halo models and if there are ways to address them at the level of this framework. For example:

• For high accuracy predictions (eg total matter power spectrum for cosmic shear), the 1-2h transition regime needs to be accounted for somehow.
• For things like cross-correlations, there needs to be some thought on how to allow profiles outside the virial radius interact with other profiles.

[DhayaaAnbajagane]

Thanks for putting all this together, @damonge ! Completely agreed that the power of this framework would be in its modularity! Having now implemented/validated the profiles from multiple "distinct" baryon halo models, it's quite obvious to me that the different profiles we use are often somewhat slightly altered versions of some common model. So being able to mix and match will be a key strength for whatever we build here.

Some loose set of thoughts on the above list:

1. General structure: Indeed I've found a common class keeps things reasonably clean while providing the user with some guidance on the types of profiles they should add together. Like you referenced above, one thing that was imperative to doing this is also being able to do arithmetic operations at just the class level, i.e. "Gas = BoundGas + EjectedGas" as an example. I have routines for this that work on the _real method and can be easily extended to work with fourier methods as well (I just did former since none of my profiles have any analytic fourier methods, so _fourier is always the FFTlog of _real). This helps clean up things significantly. There are some nice __repr__ hacking I did too so the class __str__ tells you the list of operations being composed together.
2. Parameter consistency Logistically this one is a bit of a pain since people have used the same model and params with different parameter names. The flip-side is it is indeed possible to normalize the naming so some params are shared across models. Regarding consistent updates, I have a simple routine that can recursively update params across models. The piece I haven't done is to check the updates are done correctly when using all versions of "composed profiles", where by composed I mean the "Gas = A + B" model in Point 1 above. It has all worked for me thus far but I haven't checked all edge-cases extensively.
3. Normalisation: In my experience, doing the strict normalization (i.e. the numerical integral) is not a significant time sink as far as the total runtime goes, so that could be the default option. If a model is properly specified, then M_hydro == M_dmo in the model (when eval'ed to infinity) and the (re)normalization would just be very close to 1. Therefore, our normalization step would really be most helpful for cases where the model isn't consistently specified. I don't have a sense for what use-cases exist for such mis-specified models. But it may be necessary to provide it if we want to do a simple mix/match approach.
4. Re: Normalization --- For the Mead and Arico models we have M_infinity = M200c. This is not the case for the Schneider model, and so when using the latter in a halo model P(k) you actually need to modify the HMCalculator class a bit. This is because the linear correction from Alex Mead is being implemented assuming M_total = M200c (which is reasonable) but as a result we don't obey the relation b(k=0) = 1. I have a class that makes a simple, necessary modification and show that we get the right behavior even for the Schneider19 model. I bring this up as this is the other spot where normalization of the profiles plays a role.

Happy to contribute any code from BaryonForge to be replicated/used towards these CCL efforts! That entire codebase was developed with the goal of synergizing with CCL, so this plan is along the lines of what I'd thought about in the past! My current implementation of the halo model Pk (including existing timings) is all summarized in this notebook.

[damonge]

Thanks both! (@tilmantroester , no idea why I didn't tag you above. d'uh... sorry about that).

Regarding this point:

For things like cross-correlations, there needs to be some thought on how to allow profiles outside the virial radius interact with other profiles.

What do you mean exactly? In principle we do not truncate above the virial radius (unless a given profile has this by construction). Or do you mean allowing for non-zero variance outside the virial radius in the one-halo term?

[rteyssier]

Happy to help in my capacity of old useless scientist ;-) I am a big fan of Baryon forge too

…

On Fri, Feb 21, 2025 at 12:30 PM David Alonso ***@***.***> wrote: Thanks both! ***@***.*** <https://github.com/tilmantroester> , no idea why I didn't tag you above. d'uh... sorry about that). Regarding this point: For things like cross-correlations, there needs to be some thought on how to allow profiles outside the virial radius interact with other profiles. What do you mean exactly? In principle we do not truncate above the virial radius (unless a given profile has this by construction). Or do you mean allowing for non-zero variance outside the virial radius in the one-halo term? — Reply to this email directly, view it on GitHub <#1227 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA72CGK5SC5UN57O3U2S6E32Q5PEHAVCNFSM6AAAAABXSZBHJ2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNZVGE2DKNZYGM> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***> [image: damonge]*damonge* left a comment (LSSTDESC/CCL#1227) <#1227 (comment)> Thanks both! ***@***.*** <https://github.com/tilmantroester> , no idea why I didn't tag you above. d'uh... sorry about that). Regarding this point: For things like cross-correlations, there needs to be some thought on how to allow profiles outside the virial radius interact with other profiles. What do you mean exactly? In principle we do not truncate above the virial radius (unless a given profile has this by construction). Or do you mean allowing for non-zero variance outside the virial radius in the one-halo term? — Reply to this email directly, view it on GitHub <#1227 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AA72CGK5SC5UN57O3U2S6E32Q5PEHAVCNFSM6AAAAABXSZBHJ2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNZVGE2DKNZYGM> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[nikosarcevic]

Hi, @damonge, @anicola @rreischke @DhayaaAnbajagane @rteyssier @tilmantroester et al.

fantastic suggestions!

I will give it more thought and get back to you.

In the meantime, for the rest of you (Elisa nda David know about this): there is a fully functioning Baryon Halo Model implementation that lives buried on this branch, done by myself and @elisachisari .
Example notebook is this one.

I have designed it so you can initialize the class with some params (can be extended) and you can get either the individual profiles, cross Pks or a full one.
It is slow and perhaps could have been done better, but at the same time, I am sure there are some architecture solutions that we can use in our next baryon stages. Hope you find it helpful!
I would appreciate if you can have a look and see if we can use it s a whole or at leat parts of it.

Best
Niko

[tilmantroester]

What do you mean exactly? In principle we do not truncate above the virial radius (unless a given profile has this by construction). Or do you mean allowing for non-zero variance outside the virial radius in the one-halo term?

It's been a while but IIRC, the issue we had in HMx is that for common profiles (like NFW) there is a truncation, which causes weird effects when combined with profiles outside the truncation radius. In HMx the solution was to put the ejected component into the 2h term, effectively having different profiles for the 1h and 2h terms. But there might be better solutions to this by now that I don't know about!

[DhayaaAnbajagane]

Just linking Sec 6.2.6 of the Hmx paper that shows past Tilman's (and Alex's) discussion of this point:)

I'm not aware of solutions to this issue that have come out since the HMx paper was written (though I also don't have a sense for how prevalent this DM-gas anticorrelation issue is across the parameter space). Personally, I feel the easiest thing to do is not have a sharp truncation in these profiles. We'd therefore be adopting the Schneider19 normalization formalism for the profiles in the other models as well. But my perspective is also influenced by the fact that adding the ejected mass to the 2-halo term causes problems when using the resulting halo model for baryonifying simulations (modifying the 2-halo term causes issues on large-scales when baryonifying sims)

[DhayaaAnbajagane]

Summarizing here (for record sake) a message from CCL slack to say I finally fixed up the last part of my HMx implementation in BaryonForge, and at this point the codebase contains every piece of the fiducial HMx model. We therefore have everything we need to do a like-for-like transfer of pyhmcode functionality (at least for the fiducial model) into CCL.

(This last step I implemented was about enabling the Mead2020 calculation where the ejected gas is added to just the 2-halo term and doesn't impact the 1-halo term at all.)