## 2 Classical Hamiltonian
### Lagrangian  
The associated Lagrangian which leads to the above set of Euler equations is given by  
\[\begin{array}{l}\mathcal{L}=\frac{1}{2}m\dot{q}^{2}-V(q)+\frac{1}{2}\sum \biggl{(}\dot{Q_{k}}^{2}-\omega_{k}^{2}Q_{k}{}^{2}\biggr{)}-\frac{\dot{q}}{q} \sum g_{kj}Q_{j}\dot{Q}_{k}\\ +\frac{\dot{q}^{2}}{2q^{2}}\sum\Bigl{(}h_{kj}-2g_{kj}+r_{k}\delta_{kj}\Bigr{)} Q_{k}Q_{j}.\end{array} \tag{13}\]  
It should be again noticed that \(h_{kk}=0\) by (5), which together with the antisymmetry of coefficients \(g_{kj}\) from (5), allows further simplification to reach  
\[\mathcal{L}=\frac{1}{2}m\dot{q}^{2}-V(q)+\frac{1}{2}\sum\biggl{(}\dot{Q_{k}}^{ 2}-\omega_{k}^{2}Q_{k}{}^{2}\biggr{)}+\frac{\dot{q}^{2}}{2q^{2}}\sum d_{kj}Q_ {k}Q_{j}-\frac{\dot{q}}{q}\sum g_{kj}Q_{j}\dot{Q}_{k}, \tag{14}\]  
where \(d_{kj}=\frac{1}{2}\left(h_{kj}+h_{jk}\right)+r_{k}\delta_{kj}\) is related to the symmetric part of \(h_{kj}\), given by  
\[d_{kj}=d_{jk}=\begin{cases}r_{k,}&k=j,\\ 4\frac{\left(-1\right)^{k+j}kj\left(k^{2}+\dot{q}^{2}\right)}{\left(k^{2}-\dot {q}^{2}\right)^{2}},&k\
eq j.\end{cases} \tag{15}\]  
This is to be compared with the quite different expression by [PERSON] [5] given by  
\[\mathcal{L}=\frac{1}{2}m\dot{q}^{2}-V(q)+\frac{1}{2}\sum\biggl{(}\dot{Q_{k}}^ {2}-\omega_{k}^{2}Q_{k}{}^{2}\biggr{)}+\frac{\dot{q}^{2}}{2q^{2}}\sum g_{ik}g _{lj}Q_{k}Q_{j}-\frac{\dot{q}}{q}\sum g_{kj}Q_{j}\dot{Q}_{k}, \tag{16}\]  
which indeed consistently satisfies the [PERSON]'s pair of equations.  
Now, it is a matter of speculation whether (9) is equivalent to (11) or not. Firstly, the diagonal terms are equal if  
\[\sum_{k}g_{kj}g_{kj}=r_{j}=d_{jj}, \tag{17}\]  
which holds true because of the identity  
\[\sum_{j=1,j\
eq k}^{\infty}\frac{4k^{2}\dot{q}^{2}}{\left(\dot{j}^{2}-k^{2} \right)^{2}}=k^{2}\frac{\pi^{2}}{3}+\frac{1}{4}. \tag{18}\]  
Now, for a single-mode system with only one radiation mode, we may easily notice that (9) and (11) become readily identical. This becomes more clear by noticing that \(g_{kk}=h_{kk}=0\), and hence for a system with only one electromagnetic radiation mode, (9) becomes \(\ddot{Q}=-\omega^{2}Q+\left(rq^{-2}\dot{q}^{2}\right)Q\), with \(r=r_{1}\). This is while [PERSON]'s expression (11) is \(\ddot{Q}=-\omega^{2}Q+q^{-2}\dot{q}^{2}\bigl{(}\sum g_{1j}g_{1j}\bigr{)}Q\). However, (17) requires that \(r=\sum g_{1j}g_{1j}\), which confirms the equivalency for single-mode systems.  
As for a multimode system, and by comparing (13) and (16) one would need  
\[\sum g_{kl}g_{jl}=d_{kj}=r_{k}\delta_{kj}-2g_{kj}+h_{kj}, \tag{19}\]  
in order to (9) and (11) be identical. This can be put to numerical tests (here done by the author by coding a simple Mathematica program), and is in fact accurately satisfied. Hence, the single-mode Lagrangian in the non-relativistic limit can be written in the form  
\[\mathcal{L}=\frac{1}{2}m\dot{q}^{2}-V(q)+\frac{1}{2}\bigg{[}\dot{Q}^{2}- \omega^{2}(q)Q^{2}\bigg{]}+\frac{\dot{q}^{2}}{2q^{2}}rQ^{2}. \tag{20}\]  
The last term has been usually ignored so far in the literature, and will result in momentum-field coupling. In the remainder of the paper, we focus on nonlinear terms arising from this interaction, and then also add up the relativistic corrections in the end.