Lunisolar calendar: what phase was the moon Y years ago on the spring equinox? [closed]

Question

How would I go around answering these questions (e.g., for $Y=10,000$ -- IOW, historical observations are not available and the time is too long ago to rely on a fixed calendar):

1. How many days ago ($D$) was the spring equinox $Y$ years ago?
2. What phase (how many days from conjunction) was the moon on that day $D$?

Here is my "preliminary research":

Methodology

I suppose I could solve the 3 body problem (Sun, Earth, Moon) numerically backwards in time, but I suspect that the numerical errors and gravitation perturbations from other planets will accumulate rather quickly (how quickly?)

I could do arithmetic with the durations of lunation and year, but it appears that they are not as stable as I wish they were (are they?)

Software

Online calculators (NOAA, Stellafane) do not go that far back. Presumably they rely on a version of Almagest but, again, it is not clear if these formulas are reliable for 3-5-10k year calculations.

Question:

How would I answer 1 & 2 above?

I am interested in both methodology (how it is done in general) and implementation (e.g., a program that would compute those numbers for me, given the value of $X$). It's okay to answer only a part of this, I will ask another, follow-up question.

PS. I asked this question on Math.SE and was told to "ask an astronomer".

Answer 1

This is a challenging problem and I can't give you a definitive answer, but here are some thoughts, mostly on the methodology rather than implementation.

First, you are right that fundamentally the way to do this is to integrate the equations of motion backward in time, taking into account the relevant forces. You are also right that only considering the Sun, Moon, and Earth is inadequate. The farther back in time you want to go, the more precise your model needs to be, so you need to include the influence of the other planets.

However, just considering the gravitational forces of the various bodies, but treating them as spherical masses, is also inadequate for your question. For different reasons for the two parts, you have to take into account the fact that the Earth isn't a perfect sphere.

For the equinox, this is because of precession. For determining the date of the equinox, it's not enough to know where the Earth is in its orbit; you also have to know where the pole is pointing. There are detailed models of the Earth's precession available, but I don't know how far back they go before uncertainties become large.

For the moon phase, you have to take into account the fact that interactions between the Moon and the Earth's equatorial bulge are gradually transferring angular momentum from the Earth's rotation into the Moon's orbital motion, causing the Moon to slowly move away from the Earth (and thus causing the length of the lunar month to increase slowly). Again, this has been measured and could be incorporated into a model. The rate of orbital period change is such that 10,000 years changes the Moon's period by several seconds. That doesn't sound like much, but the cumulative effect is several full orbits' difference over that time period. In other words, you'd be completely wrong about the Moon's position in its orbit on a given date if you didn't take it into account.

I suspect that there are other effects that would need to be taken into account as well to get an accurate prediction in the distant past, but those are the two that come to mind.

Regarding implementations, currently JPL Horizons is what I would consider the state of the art (though there may be others of comparable accuracy). Their time horizon is 10,000 years (i.e. you can't request positions farther back than that). I don't know if that limit is set by accuracy (which will certainly get worse as you go farther back) or by computational resources.