Firstly, we assume that the deterministic form of the chemical Langevin
equation, Eq.(REF ), has a unique stable point \(\mathbf {x}_{s}\) 
with \(\mathbf {F}(\mathbf {x}_{s})\equiv 0\) , which loses stability at the supercritical
HB \(\mu =\mu _{c}\) , where \(\mu \)  is the control parameter. Based on
the Hopf theorem [1]}, the Jacobian matrix \(\mathbf {J}\) ,
whose components \(J_{ij}=\left(\partial f_{i}/\partial x_{j}\right)|_{\mathbf {x}=\mathbf {x}_{s}}\) ,
has a pair of conjugate eigenvalues \(\lambda _{\pm }=\alpha (\mu )\pm i\omega \) 
with \(\alpha (\mu _{c})=0\) . The other \(N-2\)  eigenvalues of \(\mathbf {J}\) ,
\(-\lambda _{j(\ge 3)}\) , all have negative real parts with absolute
values considerably larger than 0. Performing the variable transformation
\(\mathbf {u}=\mathbf {T}^{-1}(\mathbf {x}-\mathbf {x}_{s})\) , the linear part of Eq.(REF )
can be transformed to Jordan form as \(\dot{\mathbf {u}}=\mathbf {\Lambda \mathbf {u}}+O(\mathbf {u}^{2})+\frac{1}{\sqrt{V}}\mathbf {\eta }\left(t\right),\) where
\(\mathbf {\Lambda }=\left(\begin{array}{cc}\alpha & -\omega \\\omega & \alpha \end{array}\right)\oplus \text{diag$\left(-\lambda _{1},\dots ,-\lambda _{N}\right)$}\)  and \(\mathbf {\eta }=\mathbf {T}^{-1}\zeta (\mathbf {x}_{s},t)\)  with \(\zeta (\mathbf {x},t)=\sum _{\rho }v_{\rho }^{j}\sqrt{w_{\rho }(\mathbf {x})}\xi _{\rho }(t)\) .
The variances of \(\mathbf {\eta }\)  are \(\left\langle \eta _{i}(t)\eta _{j}(s)\right\rangle =2D_{ij}\delta (t-s)\) 
with \(\mathbf {D}=\mathbf {T}^{-1}\mathbf {G}(\mathbf {T}^{-1})^{\text{T}}\) . When the
system locates near the HB (\(\left|\alpha \right|\ll 1\) ), the emotion
of the oscillatory mode related to \((u_{1},u_{2})\)  is much slower
than the other \(N-2\)  stable modes due to the time-scale separation.
Hence, the system's dynamics will be dominated by the slow motion
on a 2D center manifold spanned by the eigenvectors of \(\lambda _{\pm }\) .
The oscillatory mode are ruled by a normal form equation involving
the time evolution of a complex variable \(Z=u_{1}+iu_{2}\) , or a pair
of coupled equations for the oscillation amplitude \(r\)  and phase
\(\theta \)  via \(Z=re^{i\theta }\) .We follow the standard procedure to
get the normal form,
\(\frac{dZ}{dt}=\left(\alpha +i\omega \right)Z+\left(C_{r}+iC_{i}\right)\left|Z\right|^{2}Z+\frac{1}{\sqrt{V}}\sum _{\rho }\left(v_{1\rho }^{\prime }+iv_{2j\rho }^{\prime }\right)\sqrt{w_{\rho }}\xi _{\rho },\) 
