Generating rotating star clusters

[dipto4]

Hello, I was wondering if there was any way to generate rotating (uniform mass) king profile models in Agama? I read through the Agama handbook but was still a bit confused. I assume you would have to add an angular momentum component in the z-direction to the distribution function and sample from that but I am not sure how to do that in Agama.

Thank you for your help!

[eugvas]

a good question, which doesn't have a single answer. If your goal is to create an N-body snapshot of a King model, the easiest way would be to create a spherical, non-rotating model first, and then impart some amount of rotation by flipping the sign of angular momentum for some particles (a trick proposed in Lynden-Bell 1962). Details may differ - one may simply switch signs of velocity for a fixed fraction p of stars with "wrong" sign of angular momentum, up to p=1, but this introduces a discontinuity in the distribution function at Lz=0. One may do something more sophisticated, e.g. let the flipping probability change smoothly from zero to pmax as Lz increases from 0 to Lcirc. All these exercises will keep the system spherically symmetric, which is somewhat implausible - rotating systems are usually flattened to some extent. More general models for rotating non-spherical King-like models exist, e.g. http://adsabs.harvard.edu/abs/2008ApJ...689.1005B, https://adsabs.harvard.edu/abs/2012A&A...540A..94V, but they seem rather complicated in practice. Another family of models using action-space DFs is proposed in http://adsabs.harvard.edu/abs/2017MNRAS.469.4740J and is closer in spirit to the iterative construction of self-consistent models specified by action-based DFs (e.g. DoublePowerLaw), but this particular DF family is not [yet?] implemented in Agama. Perhaps @azocchi may give some additional insight on this question.

[dipto4]

Thanks a lot for the detailed explanation! The method for N-body models is quite intriguing. In order to implement it, do I just take a random fraction 'p' of stars and simply change the signs of the vx, vy, and vz components or is there something else to that?

[eugvas]

yes, just that, but only for stars with a "wrong" sign of angular momentum (otherwise it will remain non-rotating!)

[dipto4]

Thanks! The method worked really well! Just one more question: for the sophisticated method of flipping that you mentioned earlier, is there any reason why the flipping probability increases as Lz increases?

[eugvas]

If one flips velocities of all particles with even slightly negative Lz, the velocity distribution of the resulting model will be discontinuous at Vphi=0, which is probably not very realistic. OTOH if the flipping probability is gradually increasing as Lz moves away from zero, then the velocity distribution will be skewed but continuous. The width of the "transition" from no flipping to full flipping might need to be made a function of energy, e.g. full flipping occurring at Lz greater than a given fraction of Lcirc(E) - that's up to experimentation.

[eugvas]

so did you succeed in getting what you wanted? I've received a similar request from another colleague, and might implement the following feature some time in the future (but not too soon): first create a QuasiSpherical DF f(J) for the given combination of density and potential, then add a reweighting multiplicative factor g(J) = [1 + k Jphi/Jcirc(E) ] or g(J) = [1 + k tanh(Jphi/Jphi0)] or something similar, with the coefficient |k|<1 prescribing how much more weight is assigned to positive Jphi (=Lz) at the expense of subtracting the same amount from negative Jphi.

[dipto4]

I ended up using the spin-flip model that you described. Thank you so much for that! For the simulations that I am interested in, that seemed to be sufficient. Would be very interested in the implementation you just mentioned! I know that rotating king models are shown to be initially axisymmetric. Would your implementation also flatten the system?

[eugvas]

no, this will preserve the even part of the DF and hence leave the density profile unchanged. you're right that a more realistic model would need to be flattened, as discussed in these papers, but that would involve quite a bit more work: first choosing a suitable DF, then determining the density/potential profile iteratively, as in the example_self_consistent_model_flattened.py, which uses a DoublePowerLaw DF for a similar purpose.

[dipto4]

I saw that example_self_consistent_model_flattened.py uses the best fit parameters from example_doublepowerlaw.exe Would you recommend using example_doublepowerlaw with a Dehnen model?

[eugvas]

yes, why not. My strategy for constructing flattened (and possibly rotating) models following a particular density profile is the following: 1) first consider a spherical density and potential pair and construct a QuasiSpherical DF (isotropic or not); then determine the most suitable parameters of a DoublePowerLaw DF that correspond to this DF (this step is performed by example_doublepowerlaw.exe) 2) flattening is controlled by the value of "mixing coefficients" coefJzIn/coefJzOut in the DoublePowerLaw DF; a spherical model has coefJz=coefJphi, and increasing coefJz makes oblate models (flattened along the z direction). Since coefJr+coefJz+coefJphi=3 by convention, and the velocity anisotropy (radial vs.tangential) is controlled by the ratio coefJr/coefJphi, this ratio needs to be preserved when changing coefJz. Unfortunately there is no general way to determine the axis ratio z/R as a function of coefJz (inner and outer controlled separately), so it's a trial and error - create a model, export to an N-body snapshot, compute its shape (measureshape.py), adjust coefJz as needed.