--enable-tides Angular Momentum is not conserved

[veome22]

Describe the bug
When using the --enable-tides TRUE option, the total angular momentum (m_TotalAngularMomentum) of binaries is not conserved during tidal synchronization. For some binaries, the increase in angular momentum is >2% per time-step.

Label the issue
urgency_moderate - This is a moderately urgent issue
severity_moderate - This is a moderately severe bug

To Reproduce
Using COMPAS v02.41.02, run
./COMPAS --number-of-systems 1 --random-seed 2 --eccentricity-distribution FLAT --enable-tides TRUE --initial-mass-function UNIFORM

In BaseBinaryStar.cpp, print the value of m_TotalAngularMomentum before and after OmegaAfterSynchronisation() at each time-step.

Expected behavior
The tidal synchronization prescription should conserve total angular momentum when determining the final spins of the binary. The value of m_TotalAngularMomentum should, therefore, be the same during each call of OmegaAfterSynchronisation().

Screenshots
Detailed Output plot for the example binary:

Relative difference in Total and Orbital angular momentum for the binary. Note that the difference occurs when (either) one of the binary components evolves off the Main Sequence. [Edit: J_orbital plot has been updated to correctly reflect the difference after synchronization.]

Versioning (please complete the following information):

• OS: Linux 8.8 on HPC
• COMPAS v02.41.02

[jeffriley]

This statement in the tides code in BaseBinaryStar::EvaluateBinary(() (just after the root finder runs)

m_SemiMajorAxis = std::cbrt(G_AU_Msol_yr * (m_Star1->Mass() + m_Star2->Mass()) / m_Omega / m_Omega); // re-calculate semi-major axis

@ilyamandel @veome22 Can you see anything wrong with that equation? We assume circular here.

We use the root finder to calculate omega given

(a) the mass of each constituent star
(b) the moment of inertia of each constituent star
(c) the system total angular momentum

The root finder finds omega, and the omega found satisfies equation (4) in the attached document (and equation (3) if we let Lf = Li).

We then use the equation above to calculate a new semi-major axis using the omega found. But, if we then calculate the total angular momentum using the new semi-major axis, and eccentricity = 0, we get a slight difference from the total angular momentum we came in with - sometimes the difference is very small (< 0.03%), other times it is a few percent.

So either the equation above is wrong for calculating semi-major axis, or the numbers (total angular momentum? In that how it was calculated is not the same as how we calculate it to check that it was conserved here) that we come in with are slightly wrong sometimes.

[veome22]

I have been checking to see if the binary angular momentum computed using the Semi-Major Axis and Omega are consistent with each other. Right after the binary is synchronized, the relationship holds pretty well. This suggests to me that the formula itself is probably fine.

However, I noticed that BEFORE the binary is synchronized, Semi-Major Axis and Omega DO NOT give equivalent values of angular momentum. I'm guessing this is because at the beginning of the synchronization call, the semi-major axis has been updated to reflect stellar evolution, but Omega has not.

But since CalculateAngularMomentum() uses the Omegas to determine the total angular momentum, the value at the beginning is probably inconsistent with the actual angular momentum of the binary (as determined by semi-major axis). This is then fed into the OmegaAfterSynchronisation() function.

Do you think that could be an issue?

[jeffriley]

We don't evolve omega for individual stars, or stars that are constituents of a binary, unless CHE is enabled (and only then for stars that are CH stars), or, now, when tides are enabled. When CHE is not enabled, and tides are not enabled, the total angular momentum is just the orbital angular momentum of the binary. So yes, before the binary is synchronised due to tides, if it is not evolving chemically homogeneously, omega for the constituent stars has not been evolved, and the total angular momentum for the binary is calculated with non-rotating constituent stars. Yes, that could be an issue.

I'm wondering if the inconsistency of total angular momentum being for non-rotating stars sometimes, and for rotating stars at other times, is an issue we need to address - even in the documentation. @ilyamandel?

[jeffriley]

I believe the problem here is that BaseBinaryStar::CalculateAngularMomentum() uses just the mass of a star to calculate the moment-of-inertia of the star, whereas the tides code in BaseBinaryStar::EvaluateBinary() calls CalculateMomentOfInertiaAU(), and CalculateMomentOfInertiaAU() on the giant branch uses 0.1 * (mass - coremass) (when remnant mass is zero), so we get a different value for the moment-of-inertia, and so the angular momentum is also a different value.

Should we change BaseBinaryStar::CalculateAngularMomentum() to call CalculateMomentOfInertiaAU()?

@ilyamandel thoughts?

[ilyamandel]

The moments of inertia are calculated in BaseBinaryStar::CalculateAngularMomentum() from the gyration radii as follows:

double Is1  = ks1 * m1 * R1 * R1;
double Is2  = ks2 * m2 * R2 * R2;

Assuming consistent inputs are passed, this should be identical to, say,

MainSequence::CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return CalculateGyrationRadius() * m_Mass * m_Radius * m_Radius; }

[Of course, it's a good idea to check if they are consistent... And yes, I would certainly prefer using actual stellar values to passing them in as BaseBinaryStar::CalculateAngularMomentum() does -- the former is more consistent with the OOP paradigm -- but I presume there was a reason for the latter choice. ;-) ]

[ilyamandel]

Ahh, sorry, I see your point -- CalculateGyrationRadius() for HG, FGB (and, by inheritance, CHeB) stars is not defined consistently with
GiantBranch::CalculateMomentOfInertia(const double p_RemnantRadius = 0.0) const { return (0.1 * (m_Mass - m_CoreMass) * m_Radius * m_Radius) + (0.21 * m_CoreMass * p_RemnantRadius * p_RemnantRadius); } // k2 = 0.1 and k3 = 0.21 as defined in Hurley et al. 2000, after eq 109

I agree, this is definitely a problem (and likely the problem). Well spotted, @jeffriley! (I didn't check that code when reviewing the tides PR... :( )

I personally think the simplest solution is to do away with gyration radii altogether (except in internal calculations, where relevant) and just work with moments of inertia -- including passing them in directly where gyration radii are currently passed in, as in BaseBinaryStar::CalculateAngularMomentum() . Otherwise, we end up having to recompute an effective gyration radius for giant stars to be consistent with GiantBranch::CalculateMomentOfInertia().