## 3 Dark Solitons
### Gray Solitons  
Apart from the stationary black soliton (characterized by a zero density minimum), there exists the moving gray soliton solution, with \(\phi\
eq 0\) (with a nonzero density minimum). The gray soliton solution of Equation (6) reads:  
\[ \begin{split} P(\chi,\tau)\approx& 2\eta_{0}b_{0}\tanh\left[ \epsilon\eta_{0}b_{0}\sqrt{\left|\frac{q}{k_{r}^{\prime}}\right|}\left(\tau-k_{ r}^{\prime}\chi\right)-\epsilon^{2}a_{0}b_{0}q\eta_{0}^{2}\chi\right]\\ &\times\cos\left(\omega\tau-k_{r}\chi-\epsilon^{2}q\eta_{0}^{2} \chi\right)\\ &-2\eta_{0}a_{0}\sin\left(\omega\tau-k_{r}\chi-\epsilon^{2}q\eta_ {0}^{2}\chi\right)\end{split}  \tag{34}\]  
and is a function of parameters \(\chi\) and \(\tau\). In terms of the original coordinates \(x\) and \(t\), the approximate gray soliton solution for the pressure \(p(x,t)\) is given by:  
\[ \begin{split}& p(x,t)\approx 2\epsilon\eta_{0}P_{0}b_{0}\\ &\times\tanh\left[\epsilon\eta_{0}b_{0}\sqrt{\left|\frac{q}{k_{r }^{\prime}}\right|}\omega_{B}\left(t-\frac{k_{r}^{\prime}}{c_{0}}x\right)- \epsilon^{2}a_{0}b_{0}q\eta_{0}^{2}\frac{\omega_{B}}{c_{0}}x\right]\\ &\times\cos\left(\omega\omega_{B}t-\frac{k_{r}\omega_{B}}{c_{0}}x -q\epsilon^{2}\eta_{0}^{2}\frac{\omega_{B}}{c_{0}}x\right)\\ &-2\epsilon\eta_{0}P_{0}a_{0}\sin\left(\omega\omega_{B}t-\frac{k_ {r}\omega_{B}}{c_{0}}x-q\epsilon^{2}\eta_{0}^{2}\frac{\omega_{B}}{c_{0}}x \right).\end{split}  \tag{35}\]  
Figure 5: (color online) Effects of radiation and viscothermal losses on black soliton propagation in an water-filled acoustic metamaterial with \(d=0.05\) m, \(r=0.025\) m, \(r_{H}=0.002\) m, \(l_{H}=0.002\) m, calculated at (**a**) \(t=0.56\) s; (**b**) \(t=0.58\) s and (**c**) \(t=0.6\) s, respectively. A numerical dissipative black soliton (blue line) is obtained by integrating the lattice model Equation (3) with the driver corresponding to \(\epsilon=0.01\) (\(2\epsilon\eta P_{0}=4.39\times 10^{7}\) Pa) and \(f=1300\) Hz. Green dashed line is the corresponding analytical envelope results, Equation (33). The light green thick line denotes the nonlinear length \(L_{NL}\) and dispersion length \(L_{D}\), with \(L_{NL}=L_{D}=36.54\) m.  
Here, the velocity \(v_{g}\) is given by:  
\[v_{g}=\left(\frac{k_{r}^{\prime}}{c_{0}}+\frac{\epsilon a_{0}\eta_{0}\sqrt{|qk_{r }^{\prime\prime}|}}{c_{0}}\right)^{-1} \tag{36}\]  
and is the sum of the background velocity, \(v_{b}=c_{0}/k_{r}^{\prime}\), controlled by the carrier frequency, and the gray soliton's velocity, \(v_{s}=c_{0}/\left(\epsilon a_{0}\eta_{0}\sqrt{|qk_{r}^{\prime}|}\right)\), which depends on \(a_{0}=\sin\phi\) and \(\eta_{0}\). We now focus on the numerical study of the evolution of the gray soliton in the lattice. We start with the lossless case, by integrating Equation (3) with a driver of the form given by Equation (35) at \(x=0\), with \(\epsilon=0.01\) (\(2\epsilon\eta P_{0}=4.39\times 10^{7}\) Pa) and \(f=1300\) Hz. The phase angle in the simulations is chosen as \(\phi=\pi/4\), \(\phi=\pi/8\) and \(\phi=-\pi/8\)--see, respectively, Figure 6a-f. We observe that the gray solitons propagate at a constant velocity with a constant amplitude, width and minimum intensity, as shown in the contour plot in Figure 6a,c,e. The numerical spatial profiles calculated at \(t=0.35\) s (blue lines) are in good agreement with the corresponding analytical ones, Equation (35) (green dashed lines), as shown in Figure 6b,d,f.