I have read that in his Mecanique Analytique [1788], Lagrange sought a “coordinate invariant expression for mass times acceleration”.
The discussion regarding this is given in 'Marsden and Ratiu [15, Page 231]', but it includes an explanation that is too technical for me right now since it uses symplectic geomtery to explain this.
What exactly does this mean? What remains invariant here?
$F=ma$ is anyways coordinate invariant so surely that is not what is being referred to here or is it?
$\vec{F} = m\vec{a}$ is Newton's second law and it is a vector equation. That means, in 3D space, there are three components for which you have to solve for. Let us consider a simpler case, of a free particle in 2 dimensions.
Now we begin to see what coordinate invariance means. If you study the system using Cartesian coordinates, the equations of motion then are $m \ddot{x} = 0$ and $m \ddot{y} = 0$. This is fairly simple.
However, if you solve the same system in plane polar coordinates the equations of motion are NOT $m \ddot{r} = 0$ and $m \ddot{\theta} = 0$. Further, the two coordinates in the plane polar system $r, \theta$ do not even have the same dimensions. So simply replacing the cartesian coordinates with the polar coordinates in the equations of motion cannot work.
The Euler Lagrange equations of motion are invariant in the sense that as long as I have a set of coordinates {${q}$} and velocities {${\dot q}$}, you can always have equations of motion of the form:
$ \frac{\partial L}{\partial q} - \frac{d}{dt} \frac{\partial L}{\partial \dot q} = 0 $
where L is the Lagrangian of the system. So the q here could refer to coordinates in the Cartesian system or in the polar system. It is in this sense that they are coordinate invariant.
The equations that are being generalized to arbitrary coordinates are the $$\frac{dq^a}{dt} = v^a,\label{KE}\tag{Kinematic Equation}$$ $$\frac{dp_a}{dt} = f_a,\label{DE}\tag{Dynamic Equation}$$ for a system with $N$ degrees of freedom, with $0 ≤ a < N$. This doesn't just apply to case of a single body with $N = 3$ spatial coordinates, but to more general systems where the $q^a$ are its configuration coordinates. For instance, the mechanical degrees of freedom for a typical human arm, from collar on down, has at least $N = 24$ degrees of freedom, and the configuration variables would be the orientations of the hinges. If you count the spatial position and orientation of the body it's connected to, then there would be $6$ more degrees of freedom: $3$ for position, $3$ for orientation.
The two sets of equations ($\ref{KE}$) and ($\ref{DE}$), by themselves, do not complete the description of a system, but provide a framework for a description. You also need a set of constitutive relations that relate the dynamic variables $p_a$, $f_a$ to the kinematic variables $q^a$ and $v^a$. For the example you posed - a single body - you have the following equations and constitutive relations: $$v^x = \frac{dx}{dt},\quad f_x = \frac{dp_x}{dt},\quad p_x = m v^x,\quad f_x = -\frac{∂U}{∂x},$$ with similar equations for the for the $y$ and $z$ coordinates ... assuming that the forces are derivable from a potential $U(x,y,z)$.
The question, then, is: what do the $p_a$ (and $f_a$) become, when the configuration coordinates $q^a$ (and $v^a$) are transformed - where the transform may consist of general functions of the $q^a$? And what is the invariant condition or conditions that determine this transform?
For a conservative system, the $p_a$ are derivable from a function $T(v)$ as: $$p_a = \frac{∂T}{∂v^a},$$ so that they satisfy the following condition: $$\frac{∂p_a}{∂v^b} = m_{ab} = \frac{∂p_b}{∂v^a}.\label{Rpp}\tag{Reciprocal P}$$ The function $T(v)$ plays the role of kinetic energy, and the matrix of coefficients $m_{ab}$ comprise the coefficients of inertia for the system.
The $f_a$ for a conservative system are derivable from a function $U(q)$ as: $$f_a = -\frac{∂U}{∂q^a},$$ and a similar relation is satisfied: $$\frac{∂f_a}{∂q^b} = -k_{ab} = \frac{∂f_b}{∂q^a}.\label{Rff}\tag{Reciprocal F}$$ These play a role analogous to spring coefficients.
If the coefficients of inertia and "spring" coefficients were constant, then you would have a coupled harmonic oscillator ... or a generalization of one, if the $m_{ab}$ or $k_{ab}$ matrices are singular.
You can generalize further by considering the cases where the components $p_a$ of the momentum and $f_a$ of the force are functions of both the configuration coordinates $q^a$ and their corresponding velocities $v^a$. The archetypical example is the description, where $N = 3$, of a charged body with position $𝐫 = (x, y, z)$, velocity $𝐯 = \left(v^x, v^y, v^z\right)$, momentum $𝐩 = \left(p_x, p_y, p_z\right)$ and charge $e$, where the potential $U$ is velocity-dependent: $$𝐟 = -\frac{∂U}{∂𝐫},\quad U(t,𝐫,𝐯) = e\left(φ(t,𝐫) - 𝐯·𝐀(t,𝐫)\right),$$ and the momentum includes both a kinetic part proportional to $𝐯$ and a potential part proportional to $𝐀 = \left(A_x, A_y, A_z\right)$: $$𝐩 = \frac{m𝐯}{\sqrt{1 - α|𝐯|^2}} + e𝐀,$$ where $α = 0$ for the non-relativistic version and $α = 1/c^2$ for the relativistic version, where $c$ is in-vacuo light speed. The corresponding $T$ function has $𝐫$-dependence, as well as $𝐯$-dependence: $$𝐩 = \frac{∂T}{∂𝐯},\quad T(t,𝐫,𝐯) = \frac{m|𝐯|^2}{1 + \sqrt{1 - α|𝐯|^2}} + e𝐀(t,𝐫)·𝐯.$$ For this example, you also have reciprocal relations between the components of the force and momentum, e.g. $$\frac{∂p_x}{∂y} = e\frac{∂A_x}{∂y} = \frac{∂f_y}{∂v^x}.$$
Generalizing on this leads to a third set of reciprocal relations: $$\frac{∂p_a}{∂q^b} = \frac{∂f_b}{∂v^a}.\tag{Reciprocal PF}\label{Rpf}$$
Lagrange didn't have the example above on hand, since it's equivalent to the Lorentz force law, which came out of the end of the 19th century. I'm not even sure what other non-trivial example was available or considered at the time, where the force and momentum components are dependent on both configuration coordinates and velocities. So, I don't actually know where the the idea for the third set of relations ($\ref{Rpf}$) came from. Combined with ($\ref{Rpp}$) and ($\ref{Rff}$), this means that the differential: $$\sum_{0≤a<N} \left(f_a dq^a + p_a dv^a\right) = dL,$$ is (mostly, apart form subtleties relating to topology) integrable to a function $L(t,q,v)$ - the Lagrangian, itself - if the $q^a$ and $v^a$ are the independent variables. To remove the subtleties relating to topology, you can generalize this further by using Grassmann algebra on the differentials, e.g. $dx ∧ dy = -dy ∧ dx$, to write the three sets of reciprocal relations more succinctly as $$\sum_{0≤a<N} \left(df_a ∧ dq^a + dp_a ∧ dv^a\right) = 0.$$ Noting the ($\ref{KE}$) and ($\ref{DE}$), this can be treated as the application of $d/dt$ on $$ω = \sum_{0≤a<N} dp_a ∧ dq^a\\⇒\\\frac{d}{dt}ω = \sum_{0≤a<N} \left(d\left(\frac{dp_a}{dt}\right) ∧ dq^a + dp_a ∧ d\left(\frac{dq^a}{dt}\right)\right) = \sum_{0≤a<N} \left(df_a ∧ dq^a + dp_a ∧ dv^a\right),$$ while the total differential is the result of applying $d/dt$ on: $$θ = \sum_{0≤a<N} p_a dq^a\\⇒\\\frac{d}{dt}θ = \sum_{0≤a<N} \left(\frac{dp_a}{dt} dq^a + p_a d\left(\frac{dq^a}{dt}\right)\right) = \sum_{0≤a<N} \left(f_a dq^a + p_a dv^a\right).$$
Those are the basic objects of symplectic geometry, and that's how they arise. These are what are invariant. In particular, the effect of the transform of the coordinates $q^a$ to coordinates $Q^a(q)$ on the $p_a$ is a transform to components $P_a$ that satisfy the differential equation: $$\sum_{0≤a<N} P_a dQ^a = θ = \sum_{0≤a<N} p_a dq^a.$$
For a single body, that means $𝐩·d𝐫$ is the invariant that is to be preserved, if transforming from Cartesian coordinates $𝐫 = (x, y, z)$ to other coordinates, such as cylindrical or spherical. If the transform is done so as to satisfy this condition, then ($\ref{KE}$) can be used to determine what the components $V^a$ are, ($\ref{DE}$) what the components $F_a$ are, and the relations ($\ref{Rpp}$), ($\ref{Rff}$) and ($\ref{Rpf}$) will hold for the transformed quantities.
For example, if $$(x,y) = (r \cos θ, r \sin θ),$$ then $$\left(v^x, v^y\right) = \left(\cos θ v^r - r \sin θ v^θ, \sin θ v^r + r \cos θ v^θ\right),$$ and $$\begin{align} p_x dx + p_y dy &= p_x (\cos θ dr - r \sin θ dθ) + p_y (\sin θ dr + r \cos θ dθ)\\ &= \left(p_x \cos θ + p_y \sin θ\right) dr + (r \cos θ p_y - r \sin θ p_x) dθ, \end{align}$$ so that $$\left(p_r, p_θ\right) = \left(p_x \cos θ + p_y \sin θ, p_y r \cos θ - p_x r \sin θ\right) = \left(\frac{x p_x + y p_y}{r}, x p_y - y p_x\right), $$ and $$\begin{align} f_r &= f_x \cos θ - p_x \sin θ v^θ + f_y \sin θ + p_y \cos θ v^θ\\ &= \left(f_x \cos θ + f_y \sin θ\right) + \left(p_y \cos θ - p_x \sin θ\right) v^θ,\\ f_θ &= f_y r \cos θ + p_y v^r \cos θ - p_y r \sin θ v^θ - f_x r \sin θ - p_x v^r \sin θ - p_x r \cos θ v^θ\\ &= r \left(f_y \cos θ - f_x \sin θ\right) + v^r \left(p_y \cos θ - p_x \sin θ\right) - r v^θ \left(p_x \cos θ + p_y \sin θ\right). \end{align}$$
For a conservative system, the constitutive relations $$p_x = m v^x,\quad p_y = m v^y,\quad f_x = -\frac{∂U}{∂x},\quad f_y = -\frac{∂U}{∂y},$$ produce the corresponding differential equation for $L$: $$\begin{align} dL &= p_x dv^x + p_y dv^y + f_x dx + f_y dy\\ &= m\left(v^x dv^x + v^y dv^y\right) - \left(\frac{∂U}{∂x} dx + \frac{∂U}{∂y} dy\right)\\ &= d\left(m\frac{(v^x)^2 + (v^y)^2}2 - U(x,y)\right) \end{align}$$ transforms directly into $$dL = d\left(m\frac{(v^r)^2 + r^2(v^θ)^2}2 - U(r,θ)\right),$$ yielding the following transformed constitutive relations: $$\begin{align} p_r = \frac{∂L}{∂v^r} = m v^r,\quad p_θ = \frac{∂L}{∂v^θ} = m r^2 v^θ,\\ f_r = \frac{∂L}{∂r} = m r (v^θ)^2 - \frac{∂U}{∂r},\quad f_θ = \frac{∂L}{∂θ} = -\frac{∂U}{∂θ}. \end{align}$$ Notice, by the way, how the coefficients of inertia transformed: $$m_{rr} = m,\quad m_{rθ} = 0 = m_{θr},\quad m_{θθ} = m r^2.$$ The transformed components need not be constants, even if the Cartesian components were.
Edit: Also, as I look more closely at this example, I notice that the transforms to polar form $p_r$ and $p_θ$ are standard and were known in Lagrange's time, the latter being the orbital angular momentum. But, now, you have the property: $$\frac{∂p_θ}{∂r} = 2 m r v^θ = \frac{∂f_r}{∂v^θ}.$$ So, even for conservative systems, when applying transforms already well-known in Lagrange's time to polar, cylindrical or spherical coordinates, a non-trivial instance of ($\ref{Rpf}$) arises. So, that's probably where the idea for it arose. The velocity-dependent term $mr(v^θ)^2$ in $f_r$ is, of course, what's more commonly known as the "centrifugal force".
