## 3 Dark Solitons  
In this section, we use multiple scales perturbation method to describe the evolution of the pressure using an effective lossy NLS equation. From the latter equation, we can obtain the approximate analytical dark soliton solutions of our lattice system and study numerically their dynamics. We start our analysis by introducing the slow variables,  
\[\chi_{n}=\epsilon^{n}\chi,\quad\tau_{n}=\epsilon^{n}\tau,\quad n=0,1,2,\ldots, \tag{12}\]  
and express \(P\) as an asymptotic series in \(\epsilon\),  
\[P=p_{0}+\epsilon p_{1}+\epsilon^{2}p_{2}+\ldots. \tag{13}\]  
Then, substituting Equations (12) and (13) into Equation (6), we obtain a hierarchy of equations at various orders in \(\epsilon\) (see more details in Appendix C). The leading order \(O(\epsilon^{0})\), Equation (10) in Appendix C, possesses a linear plane wave solution of the form  
\[p_{0}=A(\chi_{1},\chi_{2},\cdots,\tau_{1},\tau_{2},\cdots)\exp(i\theta)+ \text{c.c.}, \tag{14}\]  
where \(A\) is an unknown envelop function, \(\theta=k_{r}\chi_{0}-\omega\tau_{0}\) with the wave number \(k_{r}\) and the frequency \(\omega\) satisfying the lossless dispersion relation, Equation (9), and c.c. denotes a complex conjugate. Next, at the order \(O(\epsilon^{1})\), the solvable condition is that the secular part (i.e., the term \(\propto\exp(i\theta)\)) vanishes, which yields,  
\[\left(k^{\prime}\frac{\partial}{\partial\tau_{1}}+\frac{\partial}{\partial \chi_{1}}\right)A(\chi_{1},\chi_{2},\cdots,\tau_{1},\tau_{2},\cdots)=0, \tag{15}\]  
where  
\[k^{\prime}=\frac{\partial k}{\partial\omega}=\frac{1}{v_{g}}=\frac{\omega}{k_ {r}-2\zeta k_{r}^{3}}, \tag{16}\]  
is the inverse group velocity. Equation (15) is satisfied as long as \(A\) depends on the variables \(\chi_{1}\) and \(\tau_{1}\) through the traveling-wave coordinate \(\tilde{\tau}_{1}=\tau_{1}-k^{\prime}\chi_{1}\), namely \(A(\chi_{1},\tau_{1},\chi_{2},\tau_{2},\cdots)=A(\tilde{\tau}_{1},\chi_{2},\tau _{2},\cdots)\). At the same order, we could obtain the form of the field \(p_{1}\),  
\[p_{1}=8\beta_{0}\frac{i\omega^{3}}{D(2\omega,2k_{r})}A^{2}e^{2i\theta}+Be^{i \theta}+c.c., \tag{17}\]  
where \(B\) is an unknown function that can be found at a higher-order approximation. Finally, following the same process as above, the nonsecularity condition at \(O(\epsilon^{2})\) yields a lossy NLS equation for the envelop function \(A\),  
\[i\frac{\partial A}{\partial\chi_{2}}-\frac{1}{2}k_{r}^{\prime\prime}\frac{ \partial^{2}A}{\partial\tau_{1}^{2}}-q\left|A\right|^{2}A=-i\Lambda A, \tag{18}\]  
where the dispersion, nonlinearity and dissipation coefficients are respectively given by,  
\[k_{r}^{\prime\prime}\equiv\frac{\partial^{2}k_{r}}{\partial\omega^{2}}=\frac {1-k_{r}^{\prime 2}+6\zeta k_{r}^{2}k_{r}^{\prime 2}}{k_{r}-2\zeta k_{r}^{3}}, \tag{19}\]  
\[q(\omega,k_{r})=\beta_{0}^{2}\frac{4\omega^{4}}{3(m^{2}+4\zeta k_{r}^{4})(k_{ r}-2\zeta k_{r}^{3})}, \tag{20}\]  
\[\Lambda=\frac{m^{2}\gamma_{H}^{\prime}}{2\omega(k_{r}-2\zeta k_{r}^{3})}. \tag{21}\]  
The sign of the product \(\sigma\equiv\text{sgn}(qk^{\prime\prime})\) determines the nature of the NLS equation [39; 40]. A focusing NLS equation with \(\sigma=+1\) supports bright solitons, localized waves with vanishing tails towards infinity. A defocusing NLS equation with \(\sigma=-1\) supports dark solitons, which are density dips, with a phase jump across the density minimum, on top of a non-vanishing continuous wave background. Figure 2c shows that the product \(qk^{\prime\prime}\) is always negative, i.e., \(\sigma=-1\) in our case, which means that our system only supports dark solitons. In the absence of losses (\(\Lambda=0\)), the analytical dark soliton solution for the envelope function \(A\) is of the form,  
\[ \begin{split} A=\eta_{0}\left\{& b_{0}\text{tanh} \left[\eta_{0}b_{0}\left(\sqrt{\left|\frac{q}{k_{r}^{\prime}}\right|}\tilde{ \tau}_{1}-a_{0}q\eta_{0}\chi_{2}\right)\right]\right.\\ &\left.+ia_{0}\right\}e^{-iq\eta_{0}^{2}\chi_{2}},\end{split}