Vtyjkyuk

The simple explanation that textbooks and the internet say is that "gravitional potential is a scalar quantity hence can be added algebraically".

However, I'm not sure if it is that simple. Take for example, at point O where it is on the horizontal axis connecting the centre of the Moon and the centre of the Earth such that the gravitational field strength at that point is zero. Any point above or below point O will have a net gravitational field strength that is the vector sum of both Earth's (g_earth) and Moon's gravitational field strength (g_moon). This vector sum will be smaller than (g_earth)+(g_moon) but larger than (g_earth) or (g_moon) alone.

Now back to the definition of gravitational potential at a point in a gravitational field, which is the work done (by an external agent) per unit mass in moving a mass from infinity to that point. This work done will be the integral of the gravitational field strength with respect to the distance away from the source of the field.

Therefore, since the gravitational field strength at point O due to the Moon and the Earth is zero, and the fact that gravitational field strength above and below point O isn't (g_earth)+(g_moon), hence the work done per unit mass in moving a mass from infinity to point O will not be as simple as adding the gravitational potential due to the Moon and Earth algebriaically.

Am I correct?

It looks like you are confusing vector magnitudes with potentials.

The potentials are indeed additive because the forces are additive, which is confirmed experientially.

If $vec F=vec F_1+vec F_2$, and we know that for a conservative force $vec F=-fracdUdx hat x$ then we have

$$-fracdUdx hat x=-fracdU_1dx hat x-fracdU_2dx hat x$$

Then we can integrate both sides with respect to $x$ to get

$$U=U_1+U_2+U_0$$

Where $U_0$ is some constant. So you see the potentials do add given that the forces add, which it seems like you agree is true. This argument works in more than one dimension as well. Notice also that none of this depends on how the magnitudes of $F_1$ and $F_2$ compare to each other or their sum.

We can also use what you reference about doing work from infinity:

$$U=-int_infty^O vec F cdot d vec x=-int_infty^O vec F_1 cdot d vec x-int_infty^O vec F_2 cdot d vec x=U_1+U_2$$

If you think about it these two methods aren't that different. The latter one just determines what $U_0$ is by setting the potential to $0$ at infinity.

Point O is a point of unstable equilibrium. This means that the potential energy is actually at a local maximum at point O. Forces just determine the slope/gradient of the potential energy, not values of the potential energy. This is evident from the above work, and knowing that you can always add a constant to the potential energy without changing the physics.

The physical reason for why we can successfully add gravitational potentials in Newtonian gravity is because it possesses no graviton self-interactions. Since (Newtonian) gravity does not interact with itself, the resultant theory is linear: so adding $2$ solutions gives us another solution to the differential equation of motion.

But Mercury's observed perihelion precession cannot be accounted for using Newtonian gravity precisely because of this absence of non-linearities. On the other hand, general relativity (GR) has gravitational self-interactions and is therefore non-linear; it correctly describes the precession. The idea that gravity interacts with itself in GR could be traced back to the idea of upholding the weak equivalence principle.

Deriving Newtonian gravity from GR involves studying the approximate version of GR called linearized general relativity where the metric is broken into a small perturbation around flat space: $g_mu nu = eta_mu nu + h_mu nu$. The equations of motion for $h$ are linear, and $h_00 = -2phi$, where $phi$ is the Newtonian gravitational potential.

In a similar vein, electromagnetism is linear (assuming all sources and currents are held fixed, so there are no back reaction effects).