## 4 Single DKS Stability Analysis  
The selection of \(|S|^{2}\) and \((\Delta)\) values is essential for achieving a stable DKS. Conducting a stability analysis is crucial, as it identifies the specific ranges of \(|S|^{2}\) and \((\Delta)\) within which the DKS state can be maintained, thereby avoiding unintended frequency comb patterns or the absence of comb generation. To analyze the stability region of the single DKS in the DCM, we propose a new method based on previous work [8,15]. Our proposed method first solves Equation (10) while allowing \(S\) or \(\Delta\) to vary. To solve Equation (10), we employ the split-step Fourier evolution method, which yields a highly stable and stationary solution. To obtain this stable solution, we use the paths proposed above, which sweeps the detuning from \(\Delta_{1}\) to \(\Delta_{2}\) while keeping \(S\) constant. Subsequently, we jump to a final point in the range of \(0<|S|^{2}<2.5\) with a resolution of 0.025 in both parameters. At each point, we wait for 3 \(\upmu\)s to evaluate the stability of DKS. \[ \begin{split} 0=-i\Bigg{[}&\frac{8\pi^{3}R}{t_{R}^{2}} \frac{\beta_{2}}{2\alpha}\Bigg{]}\frac{\partial^{2}\psi}{\partial\theta^{2}} +i|\psi|^{2}\psi-(i\Delta+1)\psi+S\\ &+i\frac{1}{\alpha}\sum_{\mu}\left[\frac{a/2}{\mu-b}\widetilde{ \psi}_{\mu}(T)e^{j\mu\theta}\right]+\left[\frac{t_{s}}{\alpha}\frac{2\pi}{T_{ R}}\right]\frac{\partial\psi}{\partial\theta}\end{split}  \tag{10}\]  
In parallel with solving Equation (10) to find the wave solutions, we determine their stability. We write  
\[\psi=\psi_{0}+\Delta\psi\qquad\quad\vec{\psi}=\psi_{0}^{*}+\vec{\Delta}\psi \tag{11}\]  
where \(\Delta\psi\) and \(\vec{\Delta}\psi\) are perturbation of the solution \(\psi_{0}\) and its complex conjugate \(\psi_{0}^{*}\). The perturbations \(\Delta\psi\) and \(\vec{\Delta}\psi\) obey the linearized equation  
\[\frac{\partial\Delta\psi}{\partial T}=\mathcal{L}\Delta\psi \tag{12}\]  
where  
\[\Delta\psi= \begin{bmatrix}\Delta\psi\\ \Delta\psi\end{bmatrix}  \tag{13}\]  
and  
\[\mathcal{L}= \begin{bmatrix}\mathcal{L}_{11}&\mathcal{L}_{12}\\ \mathcal{L}_{21}&\mathcal{L}_{22}\end{bmatrix}  \tag{14}\]  
\[\mathcal{L}_{11}\!=\!-i\!\left[\frac{8\pi^{3}\!R\beta_{2}}{2\pi t_{R}{}^{2}} \!\right]\!\frac{\partial^{2}}{\partial\theta^{2}}\!+\!2i|\psi_{0}|^{2}\!-i \Delta\!-\!1\!+\!\frac{i}{\alpha}F_{t}^{-}\!\!\left[\!\frac{a/2}{\mu-b}\! \right]\!\!F_{t}\!+\!\frac{t_{s}2\pi}{\alpha t_{R}}\frac{\partial}{\partial\theta} \tag{15}\]  
\[\mathcal{L}_{12}=i\psi_{0}^{2}\qquad\mathcal{L}_{21}=-i(\psi_{0}{}^{*})^{2} \tag{16}\]  
\[\mathcal{L}_{22}=i\!\left[\frac{8\pi^{3}\!R\beta_{2}}{2\pi t_{R}{}^{2}}\! \right]\!\frac{\partial^{2}}{\partial\theta^{2}}\!-\!2i|\psi_{0}|^{2}\!+\!i \Delta\!-\!1\!-\!\frac{i}{\alpha}F_{t}^{-}\!\!\left[\!\frac{a/2}{\mu-b}\! \right]\!\!F_{t}\!+\!\frac{t_{s}2\pi}{\alpha t_{R}}\frac{\partial}{\partial\theta} \tag{17}\]  
The terms \(F_{t}\) and \(F_{t}^{-1}\) correspond to the Fourier transform and its inverse function. Those functions enable mathematical modeling, as these transforms have a matrix representation, and the term \(\frac{a/2}{\mu-b}\) can be assigned to a diagonal matrix, resulting in a series of matrix operations of dimension \(N\times N\) that operate on the solution vector \(\psi(N\times 1)\). Next, We convert Equation (12) into an eigenvalue equation \((\mathcal{L}-\lambda I)\Delta\psi=0\), where \(I\) is the identity operator, and we discretize this equation. We then proceed to computationally find the set of all eigenvalues \(\lambda_{j}\) of this equation. If any eigenvalues are positive in its real component, the solution \(\frac{\partial\psi}{\partial\theta^{2}}\) is unstable for \(\Delta\) and \(S\) values. \(t_{s}\) is a drift term introduced in the model because the avoided crossing introduces a slight drift in the solution, causing it to become non-stationary. To find \(t_{s}\), we calculate the soliton centroid in two consecutive round trips, determining the \(\theta\) difference between both to obtain this value.