MAIN SECTION : The Lagrangian
Let express the equations of motion and the Euler-Lagrange equations with zero right hand sides
$$ \bbox[#FFFF88,12px]{\ddot{Q}_{k}+\omega^{2}_{k}Q_{k}-2\dfrac{\dot{q}}{q}\sum_{j}g_{kj}\dot{Q}_{j}-\dfrac{\ddot{q}q-\dot{q}^{2}}{q^{2}}\sum_{j}g_{kj}Q_{j}-\dfrac{\dot{q}^{2}}{q^{2}}\sum_{j\ell}g_{jk}g_{j\ell}Q_{\ell}=0} \tag{01a} $$
\begin{equation} \bbox[#FFFF88,12px]{m\ddot{q}+\dfrac{\partial V(q)}{\partial q}-\dfrac{1}{q}\sum_{k,j}(-1)^{k+j}\omega_{k}\omega_{j}Q_{k}Q_{j}=0} \tag{01b} \end{equation}
\begin{equation} \bbox[#E1FFFF,12px]{\dfrac{d}{dt}\left( \dfrac{\partial L}{\partial \dot{Q}_{k}}\right)-\dfrac{\partial L}{\partial Q_{k}}=0} \tag{02a} \end{equation}
\begin{equation} \bbox[#E1FFFF,12px]{\dfrac{d}{dt}\left( \dfrac{\partial L}{\partial \dot{q}}\right)-\dfrac{\partial L}{\partial q}=0} \tag{02b} \end{equation}
where $\:L\left(q,\dot{q},Q_{k},\dot{Q}_{k}\right)\:$ the Lagrangian.
We proceed to the following definitions in order to handle the large amount of variables and indices by means of compressed simplified expressions : \begin{equation} \mathbf{Q}\stackrel{\text{def}}{\equiv} \begin{bmatrix} Q_{1}\\ Q_{2}\\ \vdots\\ Q_{k}\\ \vdots\\ \end{bmatrix} \qquad \mathbf{\dot{Q}}\stackrel{\text{def}}{\equiv} \begin{bmatrix} \dot{Q}_{1}\\ \dot{Q}_{2}\\ \vdots\\ \dot{Q}_{k}\\ \vdots\\ \end{bmatrix} \qquad \mathbf{\ddot{Q}}\stackrel{\text{def}}{\equiv} \begin{bmatrix} \ddot{Q}_{1}\\ \ddot{Q}_{2}\\ \vdots\\ \ddot{Q}_{k}\\ \vdots\\ \end{bmatrix} \tag{03} \end{equation}
\begin{equation} \mathrm{G} \stackrel{\text{def}}{\equiv} \begin{bmatrix} 0& g_{12} & g_{13} & \cdots & g_{1k} & \cdots \\ g_{21} & 0 & g_{23} & \cdots & g_{2k} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ g_{k1} & g_{k2} & g_{k3} & \cdots & 0 & \cdots\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \end{bmatrix} =-\mathrm{G}^{\rm{T}} \tag{04} \end{equation}
\begin{equation} \Omega\left(q\right)\stackrel{\text{def}}{\equiv} \begin{bmatrix} \omega_{1} & 0 & \cdots & 0 & \cdots \\ 0 & \omega_{2} & \cdots & 0 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & \cdots & \omega_{k}& \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \end{bmatrix} = \dfrac{\pi}{q} \begin{bmatrix} 1 & 0 & \cdots & 0 & \cdots \\ 0 & 2 & \cdots & 0 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & \cdots & k & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \end{bmatrix} =\Omega^{\rm{T}}\left(q\right) \tag{05} \end{equation}
\begin{equation} \phi\left(q,\dot{q}\right)\stackrel{\text{def}}{\equiv}\dfrac{\dot{q}}{q} \tag{06} \end{equation}
We define also the real scalar below, something like the inner product of real vectors
\begin{equation} \boldsymbol{<}\mathbf{Q},\mathbf{P}\boldsymbol{>}\stackrel{\text{def}}{\equiv} \sum_{k}Q_{k}P_{k} \tag{07} \end{equation}
Under these definitions and using equations (A-01), see AUXILIARY SECTION, we have the following expressions (08) in place of the equations of motion (01) and (09) in place of (02):
\begin{equation} \mathbf{\ddot{Q}}+\Omega^{2}\left(q\right)\mathbf{Q}-2\phi\left(q,\dot{q}\right)\mathrm{G}\mathbf{\dot{Q}}-\dot{\phi}\left(q,\dot{q}\right)\mathrm{G}\mathbf{Q}+\phi^{2}\left(q,\dot{q}\right)\rm{G}^{2}\mathbf{Q}=\mathbf{0} \tag{08a} \end{equation}
\begin{equation} m\ddot{q}+\dfrac{\partial V(q)}{\partial q}-\dfrac{1}{q}\boldsymbol{<}\Omega^{2}\left(q\right)\mathbf{Q},\mathbf{Q}\boldsymbol{>}+\dfrac{1}{q}\boldsymbol{<}\Omega^{2}\left(q\right)\mathbf{Q},\mathrm{G}\mathbf{Q}\boldsymbol{>}=0 \tag{08b} \end{equation} \begin{equation} \dfrac{d}{dt}\left( \dfrac{\partial L}{\partial \mathbf{\dot{Q}}}\right)-\dfrac{\partial L}{\partial \mathbf{Q}}=\mathbf{0} \tag{09a} \end{equation} \begin{equation} \dfrac{d}{dt}\left( \dfrac{\partial L}{\partial \dot{q}}\right)-\dfrac{\partial L}{\partial q}=0 \tag{09b} \end{equation} while the Lagrangian of the system, see equation (3.1) in the question, using equations (A-02) is expressed as \begin{equation} \begin{split} & L\left(q,\dot{q},\mathbf{Q},\mathbf{\dot{Q}}\right)=\\ &\dfrac{1}{2}\boldsymbol{<}\mathbf{\dot{Q}}, \mathbf{\dot{Q}}\boldsymbol{>}-\dfrac{1}{2}\boldsymbol{<}\Omega^{2}\left(q\right)\mathbf{Q},\mathbf{Q}\boldsymbol{>}+\dfrac{1}{2}m\dot{q}^{2}-V(q)-\phi\boldsymbol{<}\mathrm{G}\mathbf{Q}, \mathbf{\dot{Q}}\boldsymbol{>}\\&-\dfrac{1}{2}\phi^{2}\boldsymbol{<}\mathrm{G}^{2}\mathbf{Q},\mathbf{Q}\boldsymbol{>} \end{split} \nonumber \end{equation} \begin{equation} \text{-----------------------------------------------------------} \tag{10} \end{equation} and in even more compact form \begin{equation} L\left(q,\dot{q},\mathbf{Q},\mathbf{\dot{Q}}\right)= \dfrac{1}{2}\left(\left\Vert \phi\mathrm{G}\mathbf{Q}-\mathbf{\dot{Q}}\right\Vert^{2}-\left\Vert\Omega\mathbf{Q}\right\Vert^{2} \right)+\dfrac{1}{2}m\dot{q}^{2}-V \tag{10$\;^{\boldsymbol{\prime}}$} \end{equation} We'll try to build the Lagrangian step by step by a trial and error procedure.
So, we expect the 1st term of equation (08a) to come from a Lagrangian part $\:L_{1}\left(\mathbf{\dot{Q}}\right)\:$ such that by (09a) \begin{equation} \dfrac{d}{dt}\left( \dfrac{\partial L_{1}}{\partial \mathbf{\dot{Q}}}\right)= \mathbf{\ddot{Q}}\Longrightarrow \dfrac{\partial L_{1}}{\partial \mathbf{\dot{Q}}}=\mathbf{\dot{Q}} \tag{11} \end{equation} From the rule (A-3.d), see AUXILIARY SECTION, $\:L_{1}\:$ is \begin{equation} L_{1}\left(\mathbf{\dot{Q}}\right)=\dfrac{1}{2}\boldsymbol{<}\mathbf{\dot{Q}}, \mathbf{\dot{Q}}\boldsymbol{>} \tag{12} \end{equation} since \begin{equation} \dfrac{\partial\left(\boldsymbol{<}\mathbf{\dot{Q}}, \mathbf{\dot{Q}}\boldsymbol{>}\right)}{\partial\mathbf{\dot{Q}}}=2\mathbf{\dot{Q}} \tag{13} \end{equation}
For the 2nd term of equation (08a) we expect a Lagrangian part $\:L_{2}\left(q,\mathbf{Q}\right)\:$ such that by (09a) \begin{equation} -\dfrac{\partial L_{2}}{\partial \mathbf{Q}}= \Omega^{2}\left(q\right)\mathbf{Q} \tag{14} \end{equation} so \begin{equation} L_{2}\left(q,\mathbf{Q}\right)=-\dfrac{1}{2}\boldsymbol{<}\Omega^{2}\left(q\right)\mathbf{Q},\mathbf{Q}\boldsymbol{>} \tag{15} \end{equation} since from the rule (A-3.c) and the symmetric (more exactly : diagonal) matrix $\:\Omega^{2}=\left(\Omega^{2}\right)^{\rm{T}}\:$ \begin{equation} \dfrac{\partial\left(\boldsymbol{<}\Omega^{2}\mathbf{Q}, \mathbf{Q}\boldsymbol{>}\right)}{\partial\mathbf{Q}}=\left[\Omega^{2}+\left(\Omega^{2}\right)^{\rm{T}}\right]\mathbf{Q}=2\Omega^{2}\mathbf{Q} \tag{16} \end{equation} But as the Lagrangian part $\:L_{2}\left(q,\mathbf{Q}\right)\:$ is a function of $\:q\:$ also, it produces items in the equations of motion if inserted to the 2nd term of (09b) : \begin{equation} -\dfrac{\partial L_{2}}{\partial q}=+\dfrac{1}{2}\dfrac{\partial\left(\boldsymbol{<}\Omega^{2}\mathbf{Q}, \mathbf{Q}\boldsymbol{>}\right)}{\partial\mathbf{Q}}=+\boldsymbol{<}\Omega \dfrac{\partial \Omega}{\partial q} \mathbf{Q},\mathbf{Q}\boldsymbol{>}=-\dfrac{1}{q}\boldsymbol{<}\Omega^{2}\left(q\right)\mathbf{Q},\mathbf{Q}\boldsymbol{>} \tag{17} \end{equation} that is exactly the 3rd term in equation (08b).
On the other hand the first two terms of (08b) are those of a particle moving in a potential, so they come from a Lagrangian part $\:L_{3}\left(q,\dot{q}\right)\:$ : \begin{equation} L_{3}\left(q,\dot{q}\right)=\dfrac{1}{2}m\dot{q}^{2}-V(q) \tag{18} \end{equation} This part $\:L_{3}\left(q,\dot{q}\right)\:$ if inserted to (9a) produces nothing (no term in equations of motion). Now, in (08a) half of the 3rd term and the 4th term give \begin{equation} -\phi\left(q,\dot{q}\right)\mathrm{G}\mathbf{\dot{Q}}-\dot{\phi}\left(q,\dot{q}\right)\mathrm{G}\mathbf{Q}=\dfrac{d}{dt}\left(-\phi\mathrm{G}\mathbf{Q}\right) \tag{19} \end{equation} so we expect a Lagrangian part $\:L_{4}\left(q,\dot{q},\mathbf{Q},\mathbf{\dot{Q}}\right)\:$ such that by (09a) \begin{equation} \dfrac{\partial L_{4}}{\partial \mathbf{\dot{Q}}}= -\phi\left(q,\dot{q}\right)\mathrm{G}\mathbf{Q} \tag{20} \end{equation} that is \begin{equation} L_{4}\left(q,\dot{q},\mathbf{Q},\mathbf{\dot{Q}}\right)=-\phi\left(q,\dot{q}\right)\boldsymbol{<}\mathrm{G}\mathbf{Q}, \mathbf{\dot{Q}}\boldsymbol{>} \tag{21} \end{equation} But, because of the antisymmetry of $\:\mathrm{G}\;$, this part may be expressed also as \begin{equation} L_{4}\left(q,\dot{q},\mathbf{Q},\mathbf{\dot{Q}}\right)=+\phi\left(q,\dot{q}\right)\boldsymbol{<}\mathrm{G} \mathbf{\dot{Q}},\mathbf{Q}\boldsymbol{>} \tag{22} \end{equation} so inserting this in the 2nd term of (09a) \begin{equation} -\dfrac{\partial L_{4}}{\partial \mathbf{Q}}= -\phi\left(q,\dot{q}\right)\mathrm{G}\mathbf{\dot{Q}} \tag{23} \end{equation} which is the other half of the 3rd term in (08a). This means that $\:L_{4}\:$, if inserted in (09a), produces the 3rd and 4th terms of (08a) \begin{equation} \dfrac{d}{dt}\left( \dfrac{\partial L_{4}}{\partial \mathbf{\dot{Q}}}\right)-\dfrac{\partial L_{4}}{\partial \mathbf{Q}}=-2\phi\left(q,\dot{q}\right)\mathrm{G}\mathbf{\dot{Q}}-\dot{\phi}\left(q,\dot{q}\right)\mathrm{G}\mathbf{Q} \tag{24} \end{equation} The output of the insertion of $\:L_{4}\:$ in (09b) would be examined later together with $\:L_{5}\:$. The 5th term of (08a) may be come from a Lagrangian part $\:L_{5}\left(q,\dot{q},\mathbf{Q},\mathbf{\dot{Q}}\right)\:$ such that by (09a)
\begin{equation} -\dfrac{\partial L_{5}}{\partial \mathbf{Q}}=\phi^{2}\left(q,\dot{q}\right)\mathrm{G}^{2}\mathbf{Q} \tag{25} \end{equation} so \begin{equation} L_{5}\left(q,\dot{q},\mathbf{Q},\mathbf{\dot{Q}}\right)=-\dfrac{1}{2}\phi^{2}\boldsymbol{<}\mathrm{G}^{2}\mathbf{Q},\mathbf{Q}\boldsymbol{>} \tag{26} \end{equation} since from (A-03.c) and the symmetry of $\:\mathrm{G}^{2}\:$ \begin{equation} \dfrac{\partial\left(\boldsymbol{<}\mathrm{G}^{2}\mathbf{Q},\mathbf{Q}\boldsymbol{>}\right)}{\partial \mathbf{Q}}=\left(\mathrm{G}^{2}+\left(\mathrm{G}^{2}\right)^{\rm{T}}\right)\mathbf{Q}=2\mathrm{G}^{2}\mathbf{Q} \tag{27} \end{equation} It can be proved, see A PROOF SECTION, that the sum $\:L_{45}=L_{4}+L_{5}\:$ \begin{equation} L_{45}\left(q,\dot{q},\mathbf{Q},\mathbf{\dot{Q}}\right)=L_{4}+L_{5}=-\dfrac{1}{2}\phi^{2}\boldsymbol{<}\mathrm{G}^{2}\mathbf{Q},\mathbf{Q}\boldsymbol{>}-\phi\left(q,\dot{q}\right)\boldsymbol{<}\mathrm{G}\mathbf{Q}, \mathbf{\dot{Q}}\boldsymbol{>} \tag{28} \end{equation} if inserted in (09b) produces the 4th term of (08b) \begin{equation} \dfrac{d}{dt}\left( \dfrac{\partial L_{45}}{\partial \dot{q}}\right)-\dfrac{\partial L_{45}}{\partial q}=+\dfrac{1}{q}\boldsymbol{<}\Omega^{2}\left(q\right)\mathbf{Q},\mathrm{G}\mathbf{Q}\boldsymbol{>} \tag{29} \end{equation}
In equation (30) below we sum up the found Lagrangian parts and the final Lagrangian is \begin{equation} \begin{split} & L\left(q,\dot{q},\mathbf{Q},\mathbf{\dot{Q}}\right)=\\ &\underbrace{\dfrac{1}{2}\boldsymbol{<}\mathbf{\dot{Q}}, \mathbf{\dot{Q}}\boldsymbol{>}}_{L_{1}}\underbrace{-\dfrac{1}{2}\boldsymbol{<}\Omega^{2}\left(q\right)\mathbf{Q},\mathbf{Q}\boldsymbol{>}}_{L_{2}}\underbrace{+\dfrac{1}{2}m\dot{q}^{2}-V(q)}_{L_{3}}\underbrace{-\phi\boldsymbol{<}\mathrm{G}\mathbf{Q}, \mathbf{\dot{Q}}\boldsymbol{>}}_{L_{4}}\\&\underbrace{-\dfrac{1}{2}\phi^{2}\boldsymbol{<}\mathrm{G}^{2}\mathbf{Q},\mathbf{Q}\boldsymbol{>}}_{L_{5}} \end{split} \nonumber \end{equation} \begin{equation} \text{-----------------------------------------------------------} \tag{30} \end{equation} identical to that given in the paper, equation (10).
Equations (31) are the equations of motion (08) with braces under items indicated from which Lagrangian terms $\:L_{m}\:$ these items come from :
\begin{equation} \underbrace{\mathbf{\ddot{Q}}}_{L_{1}}\underbrace{+\Omega^{2}\left(q\right)\mathbf{Q}}_{L_{2}}\underbrace{-2\phi\left(q,\dot{q}\right)\mathrm{G}\mathbf{\dot{Q}}-\dot{\phi}\left(q,\dot{q}\right)\mathrm{G}\mathbf{Q}}_{L_{4}}\underbrace{+\phi^{2}\left(q,\dot{q}\right)\rm{G}^{2}\mathbf{Q}}_{L_{5}}=\mathbf{0} \tag{31a} \end{equation} \begin{equation} \underbrace{m\ddot{q}+\dfrac{\partial V(q)}{\partial q}}_{L_{3}}\underbrace{-\dfrac{1}{q}\boldsymbol{<}\Omega^{2}\left(q\right)\mathbf{Q},\mathbf{Q}\boldsymbol{>}}_{L_{2}}\underbrace{+\dfrac{1}{q}\boldsymbol{<}\Omega^{2}\left(q\right)\mathbf{Q},\mathrm{G}\mathbf{Q}\boldsymbol{>}}_{L_{4}+L_{5}}=0 \tag{31b} \end{equation}
Note that the canonical momenta $\:\mathbf{P},p\:$ conjugate to $\:\mathbf{Q},q\:$ respectively are
\begin{align} \mathbf{P}&=\dfrac{\partial L}{\partial \mathbf{\dot{Q}}}=\mathbf{\dot{Q}}-\frac{\dot{q}}{q}\mathrm{G}\mathbf{Q} \tag{32a}\\ p&=\dfrac{\partial L}{\partial \dot{q}}=m\dot{q}-\dfrac{1}{q}\boldsymbol{<}\mathrm{G}\mathbf{Q},\mathbf{P}\boldsymbol{>} \tag{32b} \end{align} where for the proof of (32b) \begin{align} p&=\dfrac{\partial L}{\partial \dot{q}}=m\dot{q}-\dfrac{1}{q}\boldsymbol{<}\mathrm{G}\mathbf{Q}, \mathbf{\dot{Q}}\boldsymbol{>}-\frac{\dot{q}}{q^{2}}\boldsymbol{<}\mathrm{G}^{2}\mathbf{Q},\mathbf{Q}\boldsymbol{>} \nonumber\\ &=m\dot{q}-\dfrac{1}{q}\boldsymbol{<}\mathrm{G}\mathbf{Q}, \mathbf{\dot{Q}}\boldsymbol{>}+\frac{\dot{q}}{q^{2}}\boldsymbol{<}\mathrm{G}\mathbf{Q},\mathrm{G}\mathbf{Q}\boldsymbol{>} \nonumber\\ &=m\dot{q}-\dfrac{1}{q}\boldsymbol{<}\mathrm{G}\mathbf{Q}, \underbrace{\mathbf{\dot{Q}}-\frac{\dot{q}}{q}\mathrm{G}\mathbf{Q}}_{\mathbf{P}}\boldsymbol{>}=m\dot{q}-\dfrac{1}{q}\boldsymbol{<}\mathrm{G}\mathbf{Q},\mathbf{P}\boldsymbol{>} \tag{32b$\;^{\boldsymbol{\prime}}$} \end{align} Equations (32a) and (32b) are identical to (3.3) and (3.4) of the paper respectively, given below \begin{align} P_{k}&=\dot{Q}_{k}-\frac{\dot{q}}{q}\sum_{j}g_{kj}Q_{j} \tag{3.3}\\ p&=m\dot{q}-\dfrac{1}{q}\sum_{jk}g_{kj}P_{k}Q_{j} \tag{3.4} \end{align}