amplitude functions $d^{ml}_j(t)$ as $$V_1 \! =\!\!\!\!
\sum_{m,l,k=-\infty}^{\infty}\sum_{j,p=0}^{\infty}\frac
{4\pi^2}{D_j(m)}\Lambda^{mk}_{jp}\Psi^{ml}_j(\textbf{\textit{J}})
d^{mk}_p(t)
e^{{\rm i}\left(m\theta_\phi+l\theta_R \right)}.
\label{eq:expansion-V1-vs-dmlj}$$ It would be computationally favorable to collect $d^{ml}_j(t)$ in a single vector $\textbf{\textit{z}}(t)=\{z_n(t)\}$ by defining a map $(m,l,j)\rightarrow n$. In practice the infinite sums in (\[eq:expansion-f1\]) are truncated and approximated by finite sums so that $-l_{\rm max}\le l \le l_{\rm max}$, $-m_{\rm max}\le m\le m_{\rm max}$ and $0\le j\le
j_{\rm max}$. For $1\le n\le n_{\rm max}$, a simple map between indices will be $$\begin{aligned}
n &=& \left ( m+m_{\rm max}\right )
\left ( 2l_{\rm max}+1\right )(j_{\rm max}+1)
\nonumber \\
&{}& + \left (l+l_{\rm max}\right )(j_{\rm max}+1)+j+1,
\label{eq:map-mlj-to-i} \\
n_{\rm max} &=& \left ( 2m_{\rm max}+1\right )
      \left ( 2l_{\rm max}+1\right )
      (j_{\rm max}+1). \label{eq:maximum-n}\end{aligned}$$ One can now use (\[eq:expansion-f1\]) and (\[eq:expansion-V1-vs-dmlj\]) in (\[eq:CBE-perturbed\]) and apply the Petrov-Galerkin method to construct the weighted residual form of the CBE. That is to multiply (\[eq:CBE-perturbed\]) by some weighting functions $W^{ml}_j(\Theta,\textbf{\textit{J}})$ and to integrate the identity over the action-angle space. The outcome is the following system of nonlinear ODEs $${\rm i} { d z_p\over  d t}\!=\! \sum_{q=1}^{n_{\rm max}}
\! A_{pq} z_q \!+ \!\! \sum_{q,r=1}^{n_{\rm max}}
\! B_{pqr} z_q z_r,~~p=1,2,\cdots,n_{\rm max},
\label{eq:nonlinear-ODE}$$ for the amplitude functions $z_n(t) \equiv d^{ml}_j(t)$. The elements 