# Numerical computation of dark solitons of a nonlocal nonlinear Schrödinger equation

André de Laire, Guillaume Dujardin, Salvador López-Martínez

► To cite this version:

André de Laire, Guillaume Dujardin, Salvador López-Martínez. Numerical computation of dark solitons of a nonlocal nonlinear Schrödinger equation. Journal of Nonlinear Science, 2024, 34 (1), 10.1007/s00332-023-10001-7. http://dx.doi.org/10.1007/s00332-023-10057-9

HAL Id: hal-04058467

https://hal.science/hal-04058467v2

Submitted on 2 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
NUMERICAL COMPUTATION OF DARK SOLITONS OF A NONLOCAL
NONLINEAR SCHRÖDINGER EQUATION

ANDRÉ DE LAIRE, GUILLAUME DUJARDIN, AND SALVADOR LÓPEZ-MARTÍNEZ

ABSTRACT. The existence and decay properties of dark solitons for a large class of nonlinear nonlocal Gross-Pitaevskii equations with nonzero boundary conditions in dimension one has been established recently [10]. Mathematically, these solitons correspond to minimizers of the energy at fixed momentum and are orbitally stable. This paper provides a numerical method to compute approximations of such solitons for these types of equations, and provides actual numerical experiments for several types of physically relevant nonlocal potentials. These simulations allow us to obtain a variety of dark solitons, and to comment on their shapes in terms of the parameters of the nonlocal potential. In particular, they suggest that, given the dispersion relation, the speed of sound and the Landau speed are important values to understand the properties of these dark solitons. They also allow us to test the necessity of some sufficient conditions in the theoretical result proving existence of the dark solitons.

**AMS Classification.** 35Q55; 35C07; 35B35; 35C08; 37K40

**Keywords.** Nonlocal Schrödinger equation, Gross–Pitaevskii equation, numerical methods, numerical computations, traveling waves, dark solitons, nonzero conditions at infinity.

1. INTRODUCTION

We consider the one-dimensional nonlocal Gross–Pitaevskii equation for $\Psi: \mathbb{R} \times \mathbb{R} \to \mathbb{C}$

$$ (1) \qquad i\partial_t \Psi = \partial_x^2 \Psi + (W \star (1 - |\Psi|^2))\Psi, $$

with nonvanishing boundary conditions at infinity, i.e.

$$ (2) \qquad \lim_{|x| \to +\infty} |\Psi(x, t)| = 1, \quad \text{for all } t, $$

and appropriate initial conditions, where $\star$ denotes the convolution in space with a tempered distribution $W$. This equation appears, for instance, as a model for the evolution of a one-dimensional optical beam of intensity $|\Psi|^2$ in a self-defocusing nonlocal Kerr-like medium, where $W$ characterizes the nonlocal response of the medium [18, 15]. In this case, condition (2) is natural when studying dark optical solitons. In all the physical situations, $W$ is assumed to be a real-valued and symmetric distribution, so that (1) is a Hamiltonian equation, and the energy and momentum

$$ (3) \qquad E(\Psi(t)) = \frac{1}{2} \int_{\mathbb{R}} |\partial_x \Psi(t)|^2 dx + \frac{1}{4} \int_{\mathbb{R}} (W \star (1 - |\Psi(t)|^2))(1 - |\Psi(t)|^2) dx, $$

$$ (4) \qquad P(\Psi(t)) = \int_{\mathbb{R}} \langle i \partial_x \Psi(t), \Psi(t) \rangle_{\mathbb{C}} \left( 1 - \frac{1}{|\Psi(t)|^2} \right) dx, $$

are formally conserved, where $\langle z_1, z_2 \rangle_{\mathbb{C}} = \Re(z_1 \bar{z}_2)$, for $z_1, z_2 \in \mathbb{C}$.

In the most typical first approximation, $W$ is considered as a Dirac delta function $\delta_0$, which, from (1), leads to the standard Gross–Pitaevskii equation

$$ (5) \qquad i\partial_t \Psi = \partial_x^2 \Psi + (1 - |\Psi|^2)\Psi, $$

with nonvanishing conditions at infinity. We refer to [11, 10, 9] for more details about the theory and applications of equations (1) and (5).
Dark solitons play a fundamental role in the dynamics of finite-energy solutions to (1) propagating in a constant background. Thus, these solutions have brought theoretical and experimental attention in the last decades (see e.g. [12, 3, 7]). Precisely, these solitons correspond to smooth traveling wave solutions of the form $\Psi(x, t) = u(x - ct)$, for some speed $c \in \mathbb{R}$, with localized derivative. Therefore, we focus on finding functions $u : \mathbb{R} \to \mathbb{C}$ solving

$$ (TW_c) \qquad icu' + u'' + (W \star (1 - |u|^2))u = 0, $$

and belonging to the energy space

$$ \mathcal{E}(\mathbb{R}) = \{v \in H_{\text{loc}}^1(\mathbb{R}) : 1 - |v|^2 \in L^2(\mathbb{R}), v' \in L^2(\mathbb{R})\}. $$

Notice that this is justified by assuming that the Fourier transform of $W$ is bounded, i.e. that $\hat{W} \in L^\infty(\mathbb{R})$ (see [10]). Note that, in this paper, we use the convention that for $f \in L^1(\mathbb{R})$,

$$ \hat{f}(\xi) = \int_{\mathbb{R}} e^{-ix\xi} f(x) dx, \quad \xi \in \mathbb{R}, $$

so that $\hat{\delta}_0 = 1$. Moreover, we only consider distributions $W$ such that $\hat{W}$ is (bounded and) continuous, so that we can, and will, also assume the normalization condition $\hat{W}(0) = 1$.

Notice that any constant function $u$ of modulus 1 is a solution to $(TW_c)$ with zero energy. We refer to such functions as trivial solutions. Observe that there are also oscillating solutions to $(TW_c)$ with infinite energy such as

$$ u_{r,c}^{\pm}(x) = r \exp \left( i x \left( \frac{-c \pm \sqrt{c^2 + 4(1-r^2)}}{2} \right) \right), \quad \text{for all } r \in (0,1) \text{ and } c \in \mathbb{R}. $$

Another remark is that, by taking the complex conjugate of $u$ in equation $(TW_c)$, we only need to consider $c \ge 0$.

In the sequel, when we refer to a solution or a soliton to $(TW_c)$, we implicitly assume that it is a nontrivial finite-energy solution to $(TW_c)$. Let us recall the following result concerning the properties of these solutions.

**Lemma 1** (Decay properties of dark solitons (see [10]). Let $W$ be a real-valued even tempered distribution such that $\hat{W} \in L^\infty(\mathbb{R})$. Assume that $c \ge 0$ and that $u \in \mathcal{E}(\mathbb{R})$ is a solution to $(TW_c)$. Then $u$ is bounded and of class $C^\infty(\mathbb{R})$. Moreover, if $c > 0$, then $u$ does not vanish on $\mathbb{R}$ and there exists a smooth lifting of $u$. More precisely, there is real-valued function $\theta \in C^\infty(\mathbb{R})$ and another $C^\infty$ function $\eta$ over $\mathbb{R}$ with values in $(-\infty, 1)$ such that $u = \sqrt{1 - \eta e^{i\theta}}$ on $\mathbb{R}$, with for all $k \in \mathbb{N}$, $\theta', \eta \in H^k(\mathbb{R})$, and
$$ (6) \qquad |u|(\pm\infty) = 1, \quad D^j u(\pm\infty) = D^j \theta(\pm\infty) = D^j \eta(\pm\infty) = 0, \quad \text{for all } j \ge 1. $$

In addition, the existence of such solutions is also proved in [10], and we can rewrite this result as follows. For the sake of completeness, we provide a proof in Appendix.

**Theorem 2.** Let $W$ be a real-valued even tempered distribution such that $\hat{W} \in L^\infty(\mathbb{R})$. Assume that there is $\sigma \in (0, 1]$ such that

$$ (7) \qquad \inf_{\mathbb{R}} (\hat{W}(\xi) + \xi^2/2) = \sigma. $$

Then, for almost every $c \in (0, \sqrt{2\sigma})$, there exists a nontrivial solution $u \in \mathcal{E}(\mathbb{R})$ to $(TW_c)$.

Furthermore, by imposing some extra conditions on $W$, it is proved in [10] that the existence holds true for every $c \in (0, \sqrt{2\sigma})$. These kinds of conditions are satisfied, for instance, by the Dirac delta function $W = \delta_0$, so that $\hat{W} = 1$. This agrees with the results in [5], where the authors show that for $c \in [0, \sqrt{2})$, the unique solutions, up translation and multiplication by a complex constant of modulus one, are explicitly given by

$$ (8) \qquad u_c(x) = \sqrt{\frac{2-c^2}{2}} \tanh \left( \frac{\sqrt{2-c^2}}{2} x \right) - i \frac{c}{\sqrt{2}}. $$
In addition, they show that if $c \ge \sqrt{2}$, the only solutions to (2) are the trivial ones.

Concerning the nonlocal equation (1), as explained in [11], the Bogoliubov dispersion relation is given by

$$ (9) \qquad \omega(\xi) = \sqrt{\xi^4 + 2\hat{W}(\xi)\xi^2}. $$

Thus, assuming that $\hat{W}$ is continuous and non-negative in 0, we get $\omega(\xi) \approx (2\hat{W}(0))^{1/2}|\xi|$, for $\xi \approx 0$. We infer that the value of the sonic speed is

$$ (10) \qquad c_s(W) := \lim_{\text{def } \xi \to 0} \frac{\omega(\xi)}{\xi} = (2\hat{W}(0))^{1/2}. $$

With our normalization $\hat{W}(0) = 1$, the sonic speed is hence set to $c_s = \sqrt{2}$.

In [10], the authors proved that if $\hat{W}$ is smooth in a neighborhood of the origin, with $\hat{W} \ge 0$ a.e. on $\mathbb{R}$, $\hat{W}(0) = 1$ and $(\hat{W})''(0) \ne -1$, then $(TW_c)$ admits no nontrivial solution for the sonic speed $c_s = \sqrt{2}$.

Another important concept in the study of superfluids is the Landau's critical speed, that explains the existence of excitations in a superfluid called rotons. Mathematically, the Landau speed is defined as

$$ (11) \qquad c_L(W) := \inf_{\text{def } \mathbb{R}} \frac{\omega(\xi)}{|\xi|}, $$

and the Landau's criterion is that there exists a point $\xi_{\text{rot}} > 0$, that we call roton minimum, such that $c_L(W) = \omega(\xi_{\text{rot}})/|\xi_{\text{rot}}|$. Thus, if $\omega$ is differentiable at $\xi_{\text{rot}} > 0$, a roton minimum $\xi_{\text{rot}}$ is a point where the group velocity $\omega'(\xi_{\text{rot}})$ and the phase velocity $\omega(\xi_{\text{rot}})/\xi_{\text{rot}}$ are equal. In classical weakly interacting BEC, the Landau's speed and sonic speed are equal $c_L(W) = c_s(W)$, while in some superfluids such as $^4$He, the Landau's speed is smaller than the speed of sound [2]. In the physical literature, it is expected that when $c_L(W) < c_s(W)$, solitons solution should include some additional density oscillations around the vortex, due to the presence of a roton minimum. Moreover, Berloff and Roberts provided some formal arguments in higher dimensions, indicating that the condition $c < c_L(W)$ should be necessary to obtain nontrivial solitons [4].

To give an idea of the existence of a roton minimum, let us consider the Gaussian potential $\hat{W}_{\lambda}(\xi) = e^{-\lambda^2\xi^2}$. If $\lambda \in [0, \sqrt{2}/2)$, the dispersion curve $\xi \mapsto \omega(\xi)$ is convex, we have $c_L(W_{\lambda}) = c_s = \sqrt{2}$, and the aspect of this curve is similar to the one depicted in the left panel of Figure 1. In contrast, if $\lambda > \sqrt{2}/2$, $\omega$ is concave near the origin and there is a roton minimum $\xi_{\text{rot}} > 0$ as depicted in the right panel in Figure 1, where the slope of the green line corresponds to the Landau speed.

In this context, we can recast Theorem 2 as follows

**Corollary 3.** Let $W$ as in Theorem 2. Then the Landau's critical speed is given by

$$ (12) \qquad c_L(W) = \sqrt{2\sigma}, $$

and, for almost every $c \in (0, c_L(W))$, there exists a nontrivial solution $u \in \mathcal{E}(\mathbb{R})$ to $(TW_c)$.

As seen in Lemma 1, if $c > 0$, then solutions to $(TW_c)$ do not vanish, i.e. $|u(x)| > 0$, for all $x \in \mathbb{R}$. Moreover, as shown in [10], for $c > 0$, the change of variables, $\eta = 1 - |u|^2$ allows us to recast $(TW_c)$ as a single real-valued equation

$$ (13) \qquad J_c(\eta, \lambda) = 0, $$

where

$$ (14) \qquad J_c(\eta, \lambda) := \eta'' - 2W_{\lambda} \star \eta + c^2\eta + \frac{c^2\eta^2}{2(1-\eta)} + \frac{(\eta')^2}{2(1-\eta)} + 2(W_{\lambda} \star \eta)\eta. $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Dispersion curve $\xi \mapsto \omega_{\lambda}(\xi)$ (red), sound speed curve $\xi \mapsto \sqrt{2} \times \xi$ (blue), and Landau speed curve $\xi \mapsto c_L \times \xi$ (green), for the Gaussian potential $\tilde{W}(\xi) = e^{-\lambda^2 \xi^2}$ with $\lambda = 1/2$ (left panel), and $\lambda = 5$ (right panel). In the left panel, the sonic speed $c_s = \sqrt{2}$ is equal to the Landau speed $c_L$.</figcaption></figure>

Notice that in (14), instead of a single potential $W$, we have introduced a family of potentials $W_{\lambda}$ indexed by a real parameter $\lambda$. This will allow us to relate the nonlocal equation (13) to the local one. In fact, in several examples below, we will choose the family $W_{\lambda}$ such that $W_0 = \delta_0$.

By Lemma 1, if $u$ is a solution to (TW$_c$) with $c > 0$, then $\eta < 1$ on $\mathbb{R}$, so that $J_c$ is well-defined. Also, we can recover $u$ from $\eta$ by setting

$$ (15) \qquad u(x) = \sqrt{1 - \eta(x)} e^{i\theta(x)}, \quad \text{where} \quad \theta(x) = \frac{c}{2} \int_{0}^{x} \frac{\eta(s)}{1 - \eta(s)} ds. $$

Moreover, if $\eta$ solves (13), then the momentum (4) and the energy (3) can be computed as follows (see [10] for details)

$$ (16) \qquad E(\eta) = \frac{c^2}{8} \int_{\mathbb{R}} \frac{\eta(x)^2}{1-\eta(x)} dx + \frac{1}{8} \int_{\mathbb{R}} \frac{\eta'(x)^2}{1-\eta(x)} dx + \frac{1}{4} \int_{\mathbb{R}} (W_{\lambda} \star \eta)(x) \eta(x) dx, $$

$$ (17) \qquad P(\eta) = \frac{c}{4} \int_{\mathbb{R}} \frac{\eta(x)^2}{1-\eta(x)} dx. $$

For instance, in these new variables, the soliton $u_c$ in (8), associated with the solution to (TW$_c$) when $W = \delta_0$, reads, for $c \in (0, \sqrt{2})$,

$$ (18) \qquad \eta_c(x) = \frac{2-c^2}{2} \operatorname{sech}^2 \left( \frac{\sqrt{2-c^2}}{2} x \right), \quad \text{and} \quad \theta_c(x) = \arctan \left( \frac{\sqrt{2-c^2}}{c} \tanh \left( \frac{\sqrt{2-c^2}}{2} x \right) \right). $$

The main purpose of this article is to provide a numerical scheme to obtain approximate solutions to (TW$_c$) for some $c > 0$ and to implement this method to provide actual numerical experiments providing us with families of such dark solitons for several physically relevant potentials $W$. To do so, we introduce in Section 2 six families of physically relevant potentials $W_{\lambda}$. We also analyze the corresponding dispersion relations, sonic speeds and Landau speeds, and also discuss the existence of roton minima. Section 3 is devoted to a stability criterion which ensures orbital stability of dark solitons. We describe the numerical method we propose in Section 4 for the computation of approximation of
dark solitons (i.e. solutions to (TW$_c$)). This method uses the equivalent equation (13) in $\eta$, which is first reduced to a bounded interval with homogeneous Dirichlet boundary conditions thanks to the decay properties (Lemma 1), and then discretized by finite differences. The corresponding discrete analog to (13) is then solved numerically by minimizing a residue by gradient descent. This allows us to obtain numerical results in Section 5 for the six families $W_\lambda$ of nonlocal potentials introduced in Section 2. In particular, this gives insight into the profiles of dark solitons to (TW$_c$), and allows for an investigation of the role of the Landau speed, and the stability of solitons (in connection with Section 3).

## 2. EXAMPLES OF INTERACTION POTENTIALS

We consider six examples of interacting potentials $W_\lambda$ studied in [10]. They fulfill the hypotheses of Theorem 2 and Corollary 3. More precisely, they are even tempered distributions, with $\tilde{W} \in L^\infty(\mathbb{R})$ and $\tilde{W}$ is of class $C^2$ in a neighborhood of the origin, satisfying the normalization condition

$$ (19) \qquad \tilde{W}(0) = 1. $$

In addition, in all our examples, except the last one, the potentials are parametrized by $\lambda \in \Lambda \subset \mathbb{R}$, such that $W_0 = \delta_0$. Thus, for $\lambda = 0$, some solutions to (TW$_c$) are explicitly given by the soliton (18), for speed $c \in (0, \sqrt{2})$, that we will use later to initiate our numerical method.

We use the dispersion curve $\omega$ defined in (9), and the Landau speed in (11). When the potential is parametrized by $\lambda \in \Lambda$, we denote for simplicity $\omega_\lambda$ the corresponding dispersion curve and

$$ c_L(\lambda) = c_L(W_\lambda). $$

We remark that the sign of $\omega''(\xi)$, for small positive $\xi$ corresponds to the sign of $1 + (\tilde{W})''(\xi)$. Therefore, if $(\tilde{W})''(0) > -1$, then the curve $\omega$ is convex close to the origin, and lies above the tangent $\sqrt{2}\xi$, for small $\xi > 0$. If $(\tilde{W})''(0) < -1$, then the curve $\omega$ is concave close to the origin, and lies below the tangent $\sqrt{2}\xi$, for small $\xi > 0$. See for instance the potentials in Figure 1.

**Example 1.** Let $\beta > 0$ and $\lambda \in (-\infty, \beta/2)$, we consider

$$ (E1) \qquad W_\lambda(x) = \frac{\beta}{\beta - 2\lambda} \left( \delta_0 - \lambda e^{-\beta|x|} \right), \quad x \in \mathbb{R}, \quad \text{i.e.} \quad \tilde{W}_\lambda(\xi) = \frac{\beta}{\beta - 2\lambda} \left( 1 - \frac{2\lambda\beta}{\xi^2 + \beta^2} \right), \quad \xi \in \mathbb{R}. $$

This kind of potential has been used in [18] for the study of dark solitons in a self-defocusing nonlocal Kerr-like medium. If $\lambda > 0$, the potential $W_\lambda$ represents a strong repulsive interaction between particles that coincide in space, while the interaction becomes attractive otherwise, being this attraction more significant at short distances. In contrast, for $\lambda < 0$, the potential $W_\lambda$ is purely repulsive. According to Theorem 2 and Corollary 3, for this potential, we distinguish the following cases:

(1) If $\beta \ge \sqrt{2}$, then there is a soliton for a.e. $c \in (0, \sqrt{2})$ and for every $\lambda \in (-\infty, \beta/2)$.

(2) If $\beta \in (0, \sqrt{2})$, then there is a soliton for a.e. $c \in (0, c_L(\lambda))$ and for every $\lambda \in (-\infty, \beta/2)$, where $c_L(\lambda) = \sqrt{2\sigma(\lambda)}$ and

$$ (20) \qquad \sigma(\lambda) = \begin{cases} 1 & \text{if } -\frac{\beta^3}{2(2-\beta^2)} < \lambda < \frac{\beta}{2}, \\ \frac{\beta(1 + \sqrt{-\lambda(\beta-2\lambda)})}{\beta - 2\lambda} + \beta\sqrt{\frac{-\lambda}{\beta - 2\lambda}} - \frac{\beta^2}{2} & \text{if } \lambda \le -\frac{\beta^3}{2(2-\beta^2)}. \end{cases} $$

In particular, in this second case, as a function of $\lambda$, $c_L$ is continuous, positive and $\lim_{\lambda \to -\infty} c_L(\lambda) = (\beta(2\sqrt{2}-\beta))^{1/2}$. Notice that the critical value $-\beta^3/(2(2-\beta^2))$, corresponds to the zero of the function $1 + (\tilde{W}_\lambda)''$.

In Figure 2, we depict the dispersion curve $\omega_\lambda$ for two values of $(\beta, \lambda)$. In the left panel, we have $\beta = 0.15$ and $\lambda = 0.05$, so that $\tilde{W}_\lambda''(0) > -1$, the function $\omega$ is convex and its graph lies above that of the tangent at the origin $\xi \mapsto \sqrt{2}\xi$. There is no roton minimum and the Landau speed is $c_L = \sqrt{2}$. In the center panel, we have $\beta = 0.5$ and $\lambda = -1$, so that $\tilde{W}_\lambda''(0) < -1$, and $\omega_\lambda$ is concave near the
origin, and is below the tangent $\sqrt{2}\xi$, for $\xi \lesssim 1$. In this case, the Landau speed $c_L \sim 1.19$ is less than $\sqrt{2}$ and there is a roton for $\omega_\lambda$. For the sake of clarity, the right panel shows the curve $\omega_\lambda(\xi)/\xi$ for $\beta = 0.5$ and $\lambda = -1$.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Dispersion curves $\xi \mapsto \omega_\lambda(\xi)$ for potential (E1). Left: $\beta = 0.15$ and $\lambda = 0.05$. Center and right: $\beta = 0.5$ and $\lambda = -1$, so that $c_L \sim 1.19$.</figcaption></figure>

**Example 2.** Another interesting example proposed in [17] is the Gaussian function, for $\lambda \in \mathbb{R}\setminus\{0\}$,

$$ (E2) \qquad W_{\lambda}(x) = \frac{1}{2|\lambda|\sqrt{\pi}} e^{-\frac{x^2}{4\lambda^2}}, x \in \mathbb{R}, \quad \text{i.e.} \quad \widehat{W}_{\lambda}(\xi) = e^{-\lambda^2 \xi^2}, \xi \in \mathbb{R}, $$

with $W_0 = \delta_0$ and $\widehat{W}_0 \equiv 1$. This potential is a classical model for a smooth approximation of the Dirac delta function, as $\lambda$ goes to zero. The function $1 + (\widehat{W}_\lambda)''$ vanishes at $\lambda^* = \sqrt{2}/2$. According to Corollary 3, for $\lambda \in [0, 1/\sqrt{2})$, we have existence for almost every $c \in (0, \sqrt{2})$. Also, for $\lambda \ge \sqrt{2}/2$, we have existence of solitons for almost every $c \in (0, c_L(\lambda))$, where

$$ (21) \qquad c_L(\lambda) = \frac{1}{\lambda} \sqrt{1 + \ln(2\lambda^2)}. $$

In Figure 3, we depict the dispersion curve $\omega_\lambda$ for two values of $\lambda$. In the left and center panels, we have $\lambda = 1$, there is a roton minimum and $c_L \sim 1.3$, so the Landau speed is very close to the speed of sound. In the right panel, we have $\lambda = 3$, there is also a roton minimum and $c_L \sim 0.66$.

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Dispersion curves $\xi \mapsto \omega_\lambda(\xi)$ for potential (E2). Left and center: $\lambda = 1$ so that $c_L \sim 1.3$. Right: $\lambda = 3$ so that $c_L \sim 0.66$.</figcaption></figure>

**Example 3.** The next potential, so-called rectangular potential, was used in [1] to study supersolids, and also in [14] as a model in nonlocal materials. In addition, it can be seen as a nonsmooth
approximation of the Dirac delta when $\lambda > 0$ is small,

$$ (E3) \qquad W_{\lambda}(x) = \begin{cases} \frac{1}{2|\lambda|} & \text{if } |x| \le |\lambda| \\ 0 & \text{otherwise,} \end{cases} \quad x \in \mathbb{R}, \quad \text{i.e.} \quad \hat{W}_{\lambda}(\xi) = \text{sinc}(\lambda\xi) = \frac{\sin(\lambda\xi)}{\lambda\xi}, \quad \xi \in \mathbb{R}. $$

In this case, as a function of $\lambda$, the function $1 + (\hat{W}_{\lambda})''(0)$ vanishes at $\lambda^* = \sqrt{3}$. Moreover, the Landau speed $c_L(\lambda)$ is less than $\sqrt{2}$ for $\lambda > \sqrt{3}$. Note that $\lambda = 2$ yields a Landau speed of $c_L(2) \sim 1.374$, and $\lambda = 4.5$, yields a Landau speed of $c_L(4.5) \sim 0.624$. In Figure 4, we show the respective dispersion curves: In the left and center panels, we see that for $\lambda = 2$, there is a roton minimum. In the right panel, we see a roton minimum for $\lambda = 4.5$.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Dispersion curves $\xi \mapsto \omega_{\lambda}(\xi)$ for potential (E3). Left and center: $\lambda = 2$, so that $c_L \sim 1.374$. Right: $\lambda = 4.5$, so that $c_L \sim 0.624$.</figcaption></figure>

**Example 4.** The following potential was proposed in [20] as a simple model for interactions in a Bose-Einstein condensate. It is given by a contact interaction $\delta_0$ and two Dirac delta functions centered at $\pm \lambda$, as

$$ (E4) \qquad W_{\lambda} = 2\delta_{0} - \frac{1}{2}(\delta_{-\lambda} + \delta_{\lambda}), \quad \text{i.e.} \quad \hat{W}_{\lambda}(\xi) = 2 - \cos(\lambda\xi), \quad \xi \in \mathbb{R}. $$

For any $\lambda \ge 0$, we have $c_L(\lambda) = \sqrt{2}$ and Corollary 3 provides the existence of solitons for almost every speed $c \in (0, \sqrt{2})$. Notice the dispersion curve $\omega_{\lambda}$ can have inflection points as seen in Figure 5. However, the dispersion curve has no roton minimum, for any $\lambda \ge 0$.

**Example 5.** As pointed out in [15] in the context of solitons in nonlocal media, when the response function is narrow compared to the extent of the beam, it is possible to use the Taylor expansion

$$ \hat{W}(\xi) \approx \hat{W}(0) + \xi(\hat{W}(0))' + \frac{\xi^2}{2}(\hat{W}(0))'' = 1 - \lambda \frac{\xi^2}{2}, \quad \text{with } \lambda = \int_{\mathbb{R}} x^2 W(x) dx. $$

Here we used (19), and the fact that $W$ is even. In this case, $\hat{W}$ is not bounded and it leads to a quasilinear equation. In order to have a truly nonlocal model, we modify it by considering the the Bochner-Riesz potential given by

$$ (E5) \qquad \hat{W}_{\lambda}(\xi) = \left(1 - \lambda \frac{\xi^2}{2}\right)^{+}, \quad \xi \in \mathbb{R}. $$

From a mathematical point of view, the Bochner-Riesz potential is also interesting because it is not smooth in the Fourier variable. Even though we do not need to use the potential in the physical variable for our simulations, we give its expression to show the slow decay and oscillations at infinity:

$$ W_{\lambda}(x) = \frac{\sqrt{\lambda}}{\pi x^2} \left( \sqrt{\lambda} \frac{\sin(x\sqrt{2/\lambda})}{x} - \sqrt{2} \cos(x\sqrt{2/\lambda}) \right). $$
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Dispersion curves $\xi \mapsto \omega_{\lambda}(\xi)$ for potential (E4). Left: $\lambda = 2$. Right: $\lambda = 10$.</figcaption></figure>

In Section 5, we perform numerical simulation for $\lambda = 1$ and $\lambda = 4$. For $\lambda = 1$, we get $c_L(1) = \sqrt{2}$. This is a critical value, since the dispersion curve coincides with tangent $\sqrt{2}\xi$, for $\xi \in [0, \sqrt{2}]$, as seen in the left panel of Figure 5. For $\lambda = 4$, the Landau speed is $c_L = \sqrt{2}/2 \sim 0.7$, and there is a roton minimum, as seen in the right panel of Figure 5.

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Dispersion curves $\xi \mapsto \omega_{\lambda}(\xi)$ for (E5). Left: $\lambda = 1$. Right: $\lambda = 4$, so that $c_L = \sqrt{2}/2$.</figcaption></figure>

**Example 6.** For our last example, we consider the potential defined by its Fourier transform

$$ (E6) \qquad \hat{W}(\xi) = (1 + a\xi^2 + b\xi^4)e^{-c\xi^2}. $$

This potential was proposed in [4, 19] to describe a quantum fluid exhibiting a roton-maxon spectrum such as Helium 4. In [11], some numerical simulations were done for $a = -36$, $b = 2687$, $c = 30$, and a branch of solitons was found with speeds in $(0, \sqrt{2})$. Therefore, we will consider the same values of $(a, b, c)$, to compare our results with theirs. Also, the value of $c_L$ is approximately equal to 0.596, so that Corollary 3 gives the existence of solitons for almost every speed $c < c_L$.
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Dispersion curves for (E6), with $c_L \#0.596$.</figcaption></figure>

### 3. STABILITY OF DARK SOLITONS

In the case $W = \delta_0$, where the dark solitons are explicitly given by the formulas in (18), Lin showed in [16] their orbital stability for $c \in (0, \sqrt{2})$ (see also [8]). His proof relied on the general Grillakis-Shatah-Strauss theory [13], that could still be applied to the nonlocal evolution equation (1), to determine if the dark solitons are orbitally stable. We now recall briefly the idea of the method in our context. Let us fix $W$ and assume that there exists $c > 0$ such that for any $c \in (0, c)$, there is a solution $u(c)$ to $(TW_c)$, and that the local branch $c \in (0, c) \to u(c) \in \mathcal{E}(\mathbb{R})$ is $C^1$, so that $c \in (0, c) \to \eta(c) \in H^1(\mathbb{R})$ is $C^1$, where $\eta(c) = 1 - |u(c)|^2$.

Using (16) and (17), we can check that the functions $\tilde{E}(c) = E(\eta(c))$ and $\tilde{P}(c) = P(\eta(c))$ satisfy Hamilton group relation

$$ (22) \qquad \tilde{E}'(c) = c \tilde{P}'(c). $$

Let us assume also that the spectral assumptions in the Grillakis-Shatah-Strauss theory hold. This is probably true because of the mountain-pass theorem used to prove Theorem 2. Then, the second derivative of the function $d(c) = \tilde{E}(c) - c\tilde{P}(c)$ determines the stability. More precisely, if $d''(c) < 0$, the soliton is stable. In view of (22), we have $d''(c) = -\tilde{P}(c)$, so that the stability holds if $c \mapsto \tilde{P}(c)$ has negative derivative.

In the literature, the energy-momentum $(E, P)$ diagram is preferred to depict the local branch of solitons, and the stability is guaranteed by the strict concavity of the map $P \to E$, assuming it is of class $C^2$, since

$$ (23) \qquad \frac{dE}{dP} = c, \quad \text{and} \quad \frac{d^2E}{dP^2} = \frac{d}{dP} \left( \frac{dE}{dP} \right) = \frac{dc}{dP}. $$

An alternative method to establish the orbital stability was given in [11] by using a minimization approach, assuming that $W$ is an even tempered distribution with $\tilde{W} \in L^\infty(\mathbb{R}) \cap C^3(\mathbb{R})$, $\tilde{W} \ge 0$ on $\mathbb{R}$, with $\tilde{W}(0) = 1$ such that $\tilde{W}(\xi) \ge 1 - \xi^2/2$, for all $|\xi| < 1/\sqrt{2}$. For $q \ge 0$, they consider the minimization curve

$$ (24) \qquad E_{\min}(q) := \inf\{E(v) : v \in \mathcal{E}(\mathbb{R}), p(v) = q\}, $$

and set

$$ (25) \qquad q_* = \sup\{q > 0 \mid \forall v \in \mathcal{E}(\mathbb{R}), E(v) \le E_{\min}(q) \Rightarrow \inf_{\mathbb{R}} |v| > 0\}. $$
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Curve $E_{\min}$ and solitons in the case $W = \delta_0$.</figcaption></figure>

Theorem 4 in [11] establishes that if $E_{\min}$ is concave on $\mathbb{R}^{+}$, then the set of solutions of the minimization problem (24) is orbitally stable, for all $\mathbf{q} \in (0, \mathbf{q}_{*})$. Their proof is based on the Lions-Cazenave argument. When $W = \delta_0$, using (8), it is possible to check directly that $E_{\min}$ is smooth on $\mathbb{R}^{+}$, strictly concave on $[0, \pi/2]$, and constant on $(\pi/2, \infty)$, as depicted in Figure 8. Moreover, in view of (23), the slope of the tangent to the curve is exactly $c$, for $\mathbf{q} \in (0, \pi/2)$, and the curve remains below the critical line $\sqrt{2}\mathbf{q}$.

Therefore, both criteria explained above can be used to prove that the dark solitons are stable for every $c \in (0, \sqrt{2})$. The case $c = 0$ corresponding to the black soliton is more involved and is beyond the scope of this work. One problem is that the momentum (4) is ill-defined for perturbations of the black soliton, thus it is necessary to introduce some renormalized momentum (see [6]).

Concerning the examples in Section 2, de Laire and Mennuni were able to prove in [11] that the minimization curve $E_{\min}$ is indeed concave for the potential (E1), which implies the stability of the solitons. They conjectured that in all these examples we should have $\mathbf{q}_{*} = \pi/2$, that $E_{\min}$ should be strictly concave on $[0, \pi/2]$, and constant on $(\pi/2, \infty)$. This was corroborated by some numerical simulations in [11].

In this paper, we adopt a different approach. Indeed, we do not minimize the energy at fixed momentum, to end up with dark solitons at some speed that need to be computed afterward. In contrast, we take advantage of the equivalent form (13) of (TW$_c$), which uses implicitly parts of the results of Lemma 1, and we solve it by a simple gradient descent method on its residue, after discretization. This method is described in Section 4 and implemented on the six potentials of Section 2 in Section 5. It has at least two advantages. First, it allows to us fix the speed of the soliton *prior* to the computation. Second, it allows us to work in a vector space, since the change of unknown from $u$ in (TW$_c$) to $\eta$ in (13) took away the inhomogeneous boundary condition.

## 4. THE DISCRETIZED PROBLEM

4.1. **The discretization of** (13). For the discretization of $J_c(\cdot, \lambda)$ defined in (14), we fix some $L > 0$ and we work on the bounded closed interval $[-L/2, L/2]$. We introduce a number $N \ge 1$ of interior points, and set $\delta x = L/(N + 1)$ for $k \in \{0, \dots, N+1\}$, $x_k = k\delta x - L/2$. We
define the matrix corresponding to the centered first derivative with homogeneous Dirichlet boundary conditions as

$$ (26) \qquad D = \frac{1}{2\delta x} \begin{pmatrix} 0 & 1 & 0 & \dots & 0 \\ -1 & 0 & 1 & \dots & \vdots \\ 0 & -1 & 0 & \dots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & 0 & -1 & 0 \end{pmatrix}, $$

so that $D \in \mathcal{M}_N(\mathbb{R})$ is skew-symmetric. We also define the matrix corresponding to the classical second order derivative with homogeneous Dirichlet boundary conditions as

$$ (27) \qquad A = \frac{1}{\delta x^2} \begin{pmatrix} -2 & 1 & 0 & \dots & 0 \\ 1 & -2 & 1 & \dots & \vdots \\ 0 & 1 & -2 & \dots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & 0 & 1 & -2 \end{pmatrix}. $$

Of course, $A \in \mathcal{M}_N(\mathbb{R})$ is symmetric. For $\lambda \in \Lambda$, we define the interacting potential vector $V_\lambda \in \mathbb{R}^N$ with the entries $(V_\lambda)_k = W_\lambda(x_k)$ for $k \in \{1, \dots, N\}$, provided $W_\lambda$ is a smooth enough function.

Eventually, we define the discretized analogue of $J_c$ defined in (14) for $\eta \in \mathbb{R}^N$ and $\lambda \in \Lambda$ as

$$ (28) \quad J_c^{\delta x}(\eta, \lambda) = A\eta - 2V_\lambda * \eta + c^2\eta + c^2\eta.\eta / (2(1-\eta)) + (D\eta.D\eta) / (2(1-\eta)) + 2(V_\lambda * \eta).\eta, $$

where $\cdot$ and $./$ stand for componentwise multiplication and division respectively, and $*$ stands for the convolution of vectors. This discrete convolution corresponds to the restriction to $\{1, \dots, N\}$ of the usual convolution of the two infinite sequences consisting in its two arguments completed by zeros. This is consistent with the continuous context, since we expect the first and last components of both arguments $(V_\lambda$ and $\eta$) to be small, provided $L > 0$ has been chosen large enough. We define the numerical momentum of $\eta$ a solution to $J_c^{\delta x}(\eta, \lambda) = 0_{\mathbb{R}^N}$, by the discrete analogues to (17) and (16)

$$ (29) \qquad P^{\delta x} = \frac{c}{4} \delta x \sum_{k=1}^{N} \frac{\eta_k^2}{1 - \eta_k}, $$

$$ (30) \qquad E^{\delta x} = \frac{c^2}{8} \delta x \sum_{k=1}^{N} \frac{\eta_k^2}{1 - \eta_k} + \frac{1}{8} \delta x \sum_{k=1}^{N} \frac{(D\eta)_k^2}{1 - \eta_k} + \frac{1}{4} \delta x \sum_{k=1}^{N} (V_\lambda * \eta)_k \eta_k. $$

4.2. **The minimization algorithm.** Since we are interested in nontrivial zeros of the function $J_c$ at fixed $c \in \mathbb{R}$, we adopt the following strategy: we minimize at fixed $\lambda \in \Lambda$ and $c \in \mathbb{R}$ the squared euclidean norm of the vector-valued function $\eta \mapsto J_c^{\delta x}(\eta, \lambda)$, starting from a discretization of the known nontrivial solution (18) (for $\lambda = 0$ and for the given $c$), using a gradient method with adaptive step. This requires the introduction of the function

$$ (31) \qquad F_{\lambda}^{\delta x}(\eta) = \delta x \| J_{c}^{\delta x}(\eta, \lambda) \|^2. $$
Observe that $\eta \mapsto J_c^{\delta x}(\eta, \lambda)$ is a smooth function on some neighborhood of 0 in $\mathbb{R}^N$, and so is $F_{\lambda}^{\delta x}$. Moreover, its partial differential with respect to $\eta$ at some point $(\eta, \lambda)$ is given by

$$ (32) \quad \forall \omega \in \mathbb{R}^N, \quad \partial_{\eta} J_c^{\delta x}(\eta, \lambda)(\omega) = A\omega - 2V_{\lambda} * \omega + c^2\omega \\ \qquad + \left[ \frac{(2c^2\eta.(1-\eta/2) + D\eta.D\eta)}{(2(1-\eta).(1-\eta))} \right] \omega \\ \qquad + \left( \frac{D\eta.}{(1-\eta)} \right) \cdot (D\omega) + 2(V_{\lambda} * \eta) \cdot \omega + 2\eta \cdot (V_{\lambda} * \omega). $$

From this, we infer that

$$ (33) \quad \nabla F_{\lambda}^{\delta x}(\eta) = 2\delta x \left[ \partial_{\eta} J_{c}^{\delta x}(\eta, \lambda) \right]^{t} J_{c}^{\delta x}(\eta, \lambda), $$

where $[\partial_{\eta} J_c^{\delta x}(\eta, \lambda)]^t$ is the adjoint of the linear mapping $\omega \mapsto \partial_{\eta} J_c^{\delta x}(\eta, \lambda)(\omega)$ from $\mathbb{R}^N$ to itself defined in (32).

With the notations above, the gradient iteration starting from a given $\eta_n \in \mathbb{R}^N$ reads

$$ (34) \quad \eta_{n+1} = \eta_n - h \nabla F_{\lambda}^{\delta x}(\eta_n), $$

for some $h > 0$. After this step, the number $F_{\lambda}^{\delta x}(\eta_{n+1})$ is computed and compared to $F_{\lambda}^{\delta x}(\eta_n)$. If it is greater, then the step $h$ is divided by 2 and step (34) is performed again with the new $h$ and the same $n$. If not, then one changes $n$ in $n+1$ and performs (34) again. The method stops if (at least) one of the three criteria below is fulfilled:

(i) $F_{\lambda}^{\delta x}(\eta_n)$ is smaller than $\varepsilon^2$ for some given $\varepsilon > 0$. In this case, we consider that the method has converged.

(ii) $h$ is below some given $h_{\min} > 0$. In this case, we consider that the method has not converged.

(iii) $n$ has become greater than some given integer $n_{\max}$. In this case, we consider that the method has not converged.

Of course, if the method (34) converges, then its limit is a critical point of the function $F_{\lambda}^{\delta x}$. This point may not be a zero of the function $F_{\lambda}^{\delta x}$. This is a reason for the check of the size of $F_{\lambda}^{\delta x}$ in point (i) above. Moreover, the convergence of the method is only local in $\mathbb{R}^N$. In order to avoid convergence to the trivial critical point of $F_{\lambda}^{\delta x}$ (that is $\eta \equiv 0$), we initialize the method (34) with a projection of the continuous soliton $\eta_c$ defined in (18) on the grid $x_1, \dots, x_N$. Of course, this soliton $\eta_c$ is not a soliton for a fixed $W_{\lambda}$ in general (in particular when $\lambda \neq 0$), but it proves itself far away enough from the trivial soliton to obtain convergence to a nontrivial soliton in the numerical examples presented in the next section.

Instead of the explicit gradient method (34) to minimize (31), one could have chosen several other options, such as implicit gradient methods, or Gauss-Newton methods, for example. In this paper, we focus on (34) because it is explicit (in contrast to the other aforementioned methods, it does not require to solve a linear or nonlinear system at each step) and it proves itself stable enough in this context. Indeed, one would expect to have a CFL-type condition on $h/\delta x^3$ to ensure stability of (34). However, for the values of $\delta x$ used in the numerical experiments of Section 5, this condition appears less restrictive than the condition on $h$ implied by the condition on the sign of the difference of the values of $F_{\lambda}^{\delta x}$ between two successive iterations (i.e. $F_{\lambda}^{\delta x}(\eta_{n+1}) - F_{\lambda}^{\delta x}(\eta_n) < 0$).

## 5. NUMERICAL RESULTS

In this section, we consider the six potentials described in Section 2 and compute numerically nontrivial minimizers $\eta \in \mathbb{R}^N$ of the function $F_{\lambda}^{\delta x}$ defined in (31) for several values of $\lambda$
and $c \in (0, \sqrt{2})$, by using the algorithm described in Section 4.2. In the sequel, the number of interior points $N$ is odd and $M = (N+1)/2$ is therefore an integer corresponding to the index $k$ of the middle of the vector $\eta$ (it satisfies $x_M = 0$ with the notations of Section 4.1). Using (15), we recover from $\eta$ an approximate solution $u$ of the solution $u$ to (TW$_c$) by setting

$$ \theta_k = \begin{cases} -\frac{c}{2}\delta x \sum_{p=k}^{M-1} \frac{1}{2} \left( \frac{\eta_p}{1-\eta_p} + \frac{\eta_{p+1}}{1-\eta_{p+1}} \right) & \text{if } k \le M-1, \\ 0 & \text{if } k = M, \\ \frac{c}{2}\delta x \sum_{p=M+1}^{k} \frac{1}{2} \left( \frac{\eta_{p-1}}{1-\eta_{p-1}} + \frac{\eta_p}{1-\eta_p} \right) & \text{if } k \ge M+1, \end{cases} $$

and

$$ u_k = \sqrt{1 - \eta_k} e^{i\theta_k}, $$

for $k \in \{1, \dots, N\}$.

Recall that we initialize the method with $\eta_c$, which corresponds to the solution to (TW$_c$) for $W = \delta_0$ (see (18)). One of the difficulties is that, for values of $c$ close to $\sqrt{2}$, the size of the numerical support of the soliton $\eta$ tends to infinity.

### 5.1. Solitons with potential (E1).
We consider the potential $W_\lambda$ defined in (E1) with $\beta = 1.0$ and $\lambda = 0.4$, so that $\beta/(\beta - 2\lambda) = 5.0$. In particular, according to (20), we have $c_L(\lambda) = \sqrt{2}$, $\hat{W}''(0) > -1$, and the aspect of $\omega$ is the same as the one in the left panel of Figure 2. The numerical experiments are carried out with $L = 50$, $N = 999$, so that $\delta x = 0.05$, and $\varepsilon^2 = \delta x/4$. The parameter $h_{\min}$ is so small and the parameter $n_{\max}$ is so big that the method has always converged in these numerical experiments. The results displayed in Figure 9 allow us to see the shape of the numerical minimizers $\eta$ as functions of $x$, and the phases $\theta$, for several values of $c$. These results confirm that the support of the minimizer $\eta$ tends to spread as the speed $c$ increases. Moreover, the maximum of $\eta$ (corresponding to the minimum of $u$) tends to decrease as $c$ increases. This is similar to the case of the Dirac potential $\lambda = 0$ in this nonlocal case ($\lambda = 0.4$). The results displayed in Figure 10 show the evolution of the momentum $P^{\delta x}$ and the energy $E^{\delta x}$ (see (29) and (30)) of the numerical minimizers as functions of the speed $c$, together with the evolution of $E^{\delta x}$ as a function of $P^{\delta x}$.

Concerning the stability criteria and conjectures explained in Section 3, we see that everything is fulfilled numerically: the mapping $c \mapsto P^{\delta x}$ is strictly decreasing, the curve $P^{\delta x} \mapsto E^{\delta x}$ is below and tangent to the line $\sqrt{2}P^{\delta x}$, and is strictly concave on $(0, \pi/2)$. Therefore, these facts confirm numerically the stability of the solitons for $c \in (0, \sqrt{2})$.

Next, we consider the same potential (E1) with $\beta = 0.15$ and $\lambda = 0.05$, so that $\beta/(\beta - 2\lambda) = 3.0$ and $c_L = \sqrt{2}$. The numerical experiments are carried out with $L = 400$ and $N = 6399$, so that $\delta x = 6.25e-2$ and $\varepsilon^2 = 0.01$. The parameter $h_{\min}$ is so small and the parameter $n_{\max}$ is so big that the method has always converged in these numerical experiments. The numerical results are displayed in Figure 11. Note that, even though $c_L = \sqrt{2}$ in this case, the numerical method converged to 0 when starting from the soliton with speed $c = 1.25$, which is not displayed in the numerical results. For the other speeds, the numerical method converged to nonzero solitons, which are displayed in Figure 11. The energy $P^{\delta x}$ as a function of the momentum $P^{\delta x}$ is displayed on the left panel of Figure 12. The dispersion curve is depicted in the left panel in Figure 2, where there is no roton minimum and the computed
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Numerically computed solitons for potential (E1) with $\beta = 1.0$ and $\lambda = 0.4$ showing $\eta_k$ (left panel) and $\theta_k$ (right panel) as functions of $x_k$. Numerical parameters are indicated in the text.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 10. Energy and momentum for numerically computed solitons for potential (E1) with $\beta = 1.0$ and $\lambda = 0.4$. The first two panels show $P^{\delta x}$ and $E^{\delta x}$ as functions of $c$, while the last one depicts $E^{\delta x}$ as a function of $P^{\delta x}$. On the last panel, the dashed line corresponds to $P^{\delta x} \mapsto P^{\delta x} \sqrt{2}$.</figcaption></figure>

solitons are strictly decreasing on $\mathbb{R}^+$, as well as the curves of momentum $P^{\delta x}$ and energy $E^{\delta x}$.

We examine the potential (E1) with $\beta = 0.5$ and $\lambda = -1.0$, so that $\beta/(\beta - 2\lambda) = 0.2$, whose dispersion curve is depicted in the right panel in Figure 2. Notice that $\omega_\lambda$ has a roton minimum and that $c_L(\lambda) \sim 1.19$. The numerical experiments are carried out with $L = 60$ and $N = 2399$, so that $\delta x = 2.5e-2$ and $\varepsilon^2 = 6.25e-3$. The parameter $h_{\min}$ is so small and the parameter $n_{\max}$ is so big that the method has always converged in these numerical experiments. The numerical results are displayed in Figure 13. Note that, despite the fact that the Landau speed is $c_L \# 1.19$, the numerical method converged to a nonzero soliton for $c = 1.25$, which is displayed in the numerical results. Observe that, in contrast to the two previous experiments for potential (E1), the minimizers $\eta$ are *not* strictly decreasing on $\mathbb{R}^+$ except for $c = 1.25$, and take negative values, which means that $|u| = \sqrt{1-\eta}$ has values above 1.0. This leads thinking that the existence of a roton minimum of the dispersion curve
<figure><img src="image_12.png" /><figcaption>FIGURE 11. Numerically computed solitons for potential (E1) with $\beta = 0.15$ and $\lambda = 0.05$, showing $\eta_k$ (left panel) and $\theta_k$ (right panel) as function of $x_k$. First line: from $-L/2$ to $L/2$. Second line: from $-20$ to $+20$.</figcaption></figure>

is related to the oscillation of solitons. Remarkably, the only minimizer without oscillations corresponds to $c = 1.25$, which is above the Landau speed.

The energy $E^{\delta x}$ as a function of the momentum $P^{\delta x}$ is displayed in the right panel of Figure 12 for $(\beta, \lambda) = (0.15, 0.05)$ and $(\beta, \lambda) = (0.5, -1.0)$. We see that the existence of a roton minimum does not change the aspect of the curve $(E^{\delta x}, P^{\delta x})$, that seems to be strictly concave, so that the computed solitons should be stable.

We end this subsection by studying the effect of the parameter $\lambda$, for the same potential (E1), on the minimum value of $|u|$, and the $L^2$ and $H^1$-norms of the solitons. For this purpose, we fix $\beta = 0.5$ and $c = 0.1$, and we let $\lambda$ vary in the interval $[-10, 0.2]$. For the numerical computations, we set $L = 60$ and $N = 1201$, so that $\delta x \# 4.99e-2$ and $\varepsilon^2 \# 1.25e-2$. Numerical solitons $\eta$ and their phases $\theta$ are displayed in Figure 14. The behavior of the minimum of $u$ the soliton width, and the $L^2$ and $H^1$ norms of $\eta$ as functions of $\lambda$ are displayed in Figure 15. Let us recall that the critical value for the presence of a roton minimum in $\omega$ is $\lambda^* = -\beta^3/(2(2-\beta^2))$, that is $\lambda^* \# -0.0357$ when $\beta = 0.5$. In Figure 14, we see that the solitons exhibit a negative bump for $\lambda \le \lambda^*$, and that they are strictly decreasing for $\lambda > \lambda^*$,
<figure><img src="image_13.png" /><figcaption>FIGURE 12. $E^{\delta x}$ as a function of $P^{\delta x}$ for several speeds. Left: $\beta = 0.15$ and $\lambda = 0.05$. Right: $\beta = 0.5$ and $\lambda = -1.0$. On both panels, the dashed line corresponds to $P^{\delta x} \mapsto P^{\delta x}\sqrt{2}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 13. Numerically computed solitons for potential (E1) with $\beta = 0.5$ and $\lambda = -1.0$, showing $\eta_k$ (left panel) and $\theta_k$ (right panel) as function of $x_k$.</figcaption></figure>

supporting the link between a roton minimum and the oscillations. In addition, we remark a monotonic behavior on the minimum value of $|u|$, on the width of $\eta$, and on the $L^2$ and $H^1$-norms of the solitons, independently of the presence of a roton minimum.

5.2. **Solitons with potential (E2).** We consider now the Gaussian potential $W_\lambda$ defined in (E2) in Example 2. We take first $\lambda = 1.0$ and compute numerically nontrivial minimizers of the function $F_\lambda^{\delta x}$ defined in (31) with $\lambda = 1.0$ and several values of $c$ using the algorithm described in Subsection 4.2. The numerical experiments are carried out with $L = 50$, $N = 999$, so that $\delta x = 0.05$, and $\varepsilon^2 = \delta x/4$. In this case, $\delta\xi = 2\pi/L\#0.13$. The results are displayed in Figures 16 and 17. As explained in Example 2, we expect solutions for every $c \in (0, c_L(1))$, where $c_L(1) = \#1.30$, in agreement with Figures 16 and 17. Moreover, the presence of a roton minimum (see middle panel in Figure 3) should imply the presence of oscillations on
<figure><img src="image_15.png" /><figcaption>FIGURE 14. Numerically computed solitons for potential (E1) with $\beta = 0.5$ and $c = 0.1$ for several values of $\lambda$. Top panel: $\eta_k$ (left) and $\theta_k$ (right) as function of $x_k$. Bottom panel is the legend of the top panel.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 15. Minimum of $u$ (left panel), width of the soliton (middle panel), and $L^2$ and $H^1$ norms (right panel) as a function of $\lambda$ for potential (E1) with $\beta = 0.5$ and $c = 0.1$. Numerical parameters are indicated in the text.</figcaption></figure>

$\mathbb{R}^+$, which is indeed seen in these figures. Notice that for $c = 1.25$, $\eta_k$ shows no oscillation on $\mathbb{R}^+$, but there is a bump in its phase $\theta_k$, so there is an oscillation of $u_k$.

Next, we run another experiment with a larger $\lambda$ to make the potential more long range. We chose $\lambda = 3$, and run another series of computation of numerical minimizers for several values of $c$ as before. As seen in Example 2, there is a roton minimum and the Landau speed is $c_L(3)\#0.66$. To take into account the larger supports of the minimizers, we take $L = 100$ and $N = 1999$ so that $\delta x = 0.05$ and $\epsilon^2 = 1.25e-2$. The results are displayed in Figures 18
<figure><img src="image_17.png" /><figcaption>FIGURE 16. Numerically computed solitons for potential (E2) with $\lambda = 1.0$, showing $\eta_k$ (left panel) and $\theta_k$ (right panel) as function of $x_k$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 17. Numerically computed solitons for potential (E2) with $\lambda = 1.0$. On the left: $P^{\delta x}$ as a function of $c$. In the middle: $E^{\delta x}$ as a function of $c$. On the right: $E^{\delta x}$ as a function of $P^{\delta x}$.</figcaption></figure>

and 19. We have thus obtained soliton for speeds above the Landau speed, which seem to be stable, due to the aspect of the curves in Figure 19. This is an unexpected result, that shows that the physical conjecture of the role of the Landau speed is not valid in this case. However, the fact that there is bigger gap between $c_L$ and $\sqrt{2}$ (than for $\lambda = 1$), could be an indication of the more oscillating behavior of the solitons in Figure 18.

5.3. **Solitons with potential (E3).** We study now the rectangular potential $W_\lambda$ defined in (E3) in Example 3. We start by computing numerically nontrivial minimizers of the function $F_\lambda^{\delta x}$ in (31) for $\lambda = 2.0$ and several values of $c$ as described in Subsection 4.2. The numerical experiments are carried out with $L = 50$, $N = 999$, so that $\delta x = 0.05$, and $\varepsilon^2 = \delta x/4$. In this case, $\delta\xi = 2\pi/L \sim 0.13$. The results are displayed in Figures 20 and 21. As seen in Figure 4, there is a roton minimum and the Landau speed $c_L$#1.374 is very close to the speed of sound. This agrees of the small bump on $\eta_k$, except for the soliton with speed $c = 1.25$, whose profile $\eta_k$ and phase $\eta_k$ seem to be strictly monotone on $\mathbb{R}^+$. As in all previous examples, the monotonicity of $c \mapsto P^\delta$ indicates that these are stable solitons.
<figure><img src="image_19.png" /><figcaption>FIGURE 18. Numerically computed solitons for potential (E2) with $\lambda = 3.0$, showing $\eta_k$ (left panel) and $\theta_k$ (right panel) as function of $x_k$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 19. Numerically computed solitons for potential (E2) with $\lambda = 3.0$, showing $P^{\delta x}$ (left) and $E^{\delta x}$ (center) as a function of $c$. On the right, $E^{\delta x}$ as a function of $P^{\delta x}$.</figcaption></figure>

We perform another experiment for $\lambda = 4.5$, with $L = 80$, $N = 1599$ and a tolerance of $\varepsilon^1 = \delta x/4 = 1.25\varepsilon - 2$. The results are displayed in Figures 22 and 23.

In this case, $c_L \# .624$ and the dispersion curve is depicted in Figure 4. We found a change on the behavior of soliton similar to one detected in Subsection 5.2, i.e. a more oscillating behavior of the solitons. Concerning the stability, the curve $P^{\delta x}$ is strictly decreasing for $c \le 075$, but it is not clear the exact behavior of the curve for $c > 075$.

**5.4. Solitons with potential (E4).** In this subsection, we consider the potential $W_\lambda$ in (E4) in Example 4, given by three Dirac delta functions, so that the Landau speed is $c_L(\lambda) = \sqrt{2}$, for any $\lambda \ge 0$, and there is no roton minimum. We will investigate numerical soliton for $\lambda = 2$ and $\lambda = 10$, whose dispersion curve are given in Figure 5.

We first compute numerically nontrivial minimizers of $F_\lambda^{\delta x}$ in (31) for $\lambda = 2.0$ using the algorithm described in Subsection 4.2. The numerical experiments are carried out with $L = 50$, $N = 999$, so that $\delta x = 0.05$, and $\varepsilon^2 = \delta x/4$. In this case, $\delta\xi = 2\pi/L \sim 0.13$. The results are displayed in Figures 24 and 25. As expected, the aspect of the dispersion curve in the left
<figure><img src="image_21.png" /><figcaption>FIGURE 20. Numerically computed solitons for potential (E3) with $\lambda = 2.0$, showing $\eta_k$ (left panel) and $\theta_k$ (right panel) as a function of $x_k$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>FIGURE 21. Numerically computed solitons for potential (E3) with $\lambda = 2.0$, showing $P^{\delta x}$ (left) and $E^{\delta x}$ (center) as a function of $c$. On the right, $E^{\delta x}$ as a function of $P^{\delta x}$.</figcaption></figure>

panel of Figure 5 agrees with the monotonicity of the numerical solitons in Figure 24, which should be stable in view of the results in Figure 25.

We perform another experiment for $\lambda = 10.0$ with $L = 300$, $N = 5999$, so that $\delta x = 0.05$ and $\varepsilon^2/\delta x/4$. In this case, $\delta\xi = 2\pi/L \sim 0.021$. The numerical results are displayed in Figures 26 and 27. Even though there is no roton minimum and the Landau speed coincides with the speed of sound, we see an oscillating behavior of the solutions. In view of the dispersion curve in Figure 5, we conjecture that the oscillation are related to the existence of inflection points of the dispersion curve. Once again, these solitons should be stable in view of the curves in Figure 27.

Another point is that, despite the oscillations in $\eta$ in Figure 26, the minimum values are of order *minus* a few $10^{-3}$ for all speed. This may indicate that the corresponding continuous solitons $\eta$ are indeed non-negative, so that the corresponding function $u$ has modulus bounded by 1.

5.5. **Solitons with potential (E5).** In this subsection, we consider the Bochner-Riesz potential $W_\lambda$ defined in (E5), for $\lambda = 2$ and $\lambda = 4$, as explained in Example 5.
<figure><img src="image_23.png" /><figcaption>FIGURE 22. Numerically computed solitons for potential (E3) with $\lambda = 4.5$, showing $\eta_k$ (left) and $\theta_k$ (right) as a function of $x_k$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIGURE 23. Numerically computed solitons for potential (E3) with $\lambda = 4.5$, showing $P^{\delta x}$ (left) and $E^{\delta x}$ (center) as a function of $c$. On the right, $E^{\delta x}$ as a function of $P^{\delta x}$.</figcaption></figure>

We compute first nontrivial minimizers of the function $F_{\lambda}^{\delta x}$ for $\lambda = 1.0$, as before, taking $L = 200$ and $N = 3999$, so that $\delta x = 0.05$, and $\varepsilon^2 = \delta x/4$. In this case, $\delta\xi = 2\pi/L \sim 0.13$. The results are displayed in Figures 28 and 29. Clearly, $\lambda = 1$ is a critical case for the dispersion curve, since it is given by a straight line in interval $[0, \sqrt{2}]$, as seen in Figure 6. Even though $\omega_{\lambda}$ has no inflection point, the non-differentiability at $\xi = \sqrt{2}$ could be related to the oscillations of solitons in Figure 28. These solitons should be stable in view of Figure 29.

We run another experiment with $\lambda = 4.0$, $L = 300$, $N = 5999$, so that $\delta x = 0.05$. The results are displayed in Figures 30 and 31. Even though in this case $c_L\#0.7$, we obtain again solitons for every subsonic speed, that seem stable since $c \mapsto P^{\delta x}$ is decreasing. The more oscillating behavior of the soliton, could be explained by the gap between $c_L$ and $\sqrt{2}$, as in Subsection 5.3.

5.6. **Solitons with potential (E6).** We consider the potential $W$ defined in (E6). The minimizers are computed with $L = 800$ and $N = 7999$ so that $\delta x = 0.1$. The results are displayed in Figures 32 and 33. As seen in Figure 7, there is a roton minimum and $c_L\#0.596$. This is no obstacle for the computation of soliton for every subsonic speed. Since $P^{\delta x}$ is a
<figure><img src="image_25.png" /><figcaption>FIGURE 24. Numerically computed solitons for potential (E4) with $\lambda = 2.0$, showing $\eta_k$ (left) and $\theta_k$ (right) as a function of $x_k$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>FIGURE 25. Numerically computed solitons for potential (E4) with $\lambda = 2.0$, showing $P^{\delta x}$ (left) and $E^{\delta x}$ (center) as a function of $c$. On the right, $E^{\delta x}$ as a function of $P^{\delta x}$.</figcaption></figure>

strictly decreasing function, these solitons should be stable. The highly oscillating behavior could be explained again by the gap between the Landau speed and the speed of sound.

## 6. CONCLUSION

This paper introduces a numerical method to compute dark solitons to the nonlocal Gross-Pitaevskii equation with nonlocal realistic potentials, with a prescribed speed. Using an *ad hoc* formulation of the problem (see (13)-(14), and the discrete analogue (28)), one has to solve a nonlinear nonlocal problem ($J_c^{\delta x}(\eta, \lambda) = 0_{\mathbb{R}^N}$ for a fixed parameter value $\lambda \in \mathbb{R}$) in the neighborhood of a guess, which is linked with the known continuous soliton with same speed for $\lambda = 0$. This problem is solved numerically by a gradient descent method minimizing the residue, and allows for numerical simulations for several physically realistic nonlocal interaction potentials.

These numerical simulations provide insight into the behavior of solitons and allow us to discuss the influence of the nonlocal interactions on the shape of the dark solitons, as well as on their stability. Moreover, they suggest that the speed of sound and the Landau speed,
<figure><img src="image_27.png" /><figcaption>FIGURE 26. Numerically computed solitons for potential (E4) with $\lambda = 10.0$, showing $\eta_k$ (left) and $\theta_k$ (right) as a function of $x_k$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>FIGURE 27. Numerically computed solitons for potential (E4) with $\lambda = 10.0$, showing $P^{\delta x}$ (left) and $E^{\delta x}$ (right) as a function of $c$. On the right, $E^{\delta x}$ as a function of $P^{\delta x}$.</figcaption></figure>

given by the dispersion relation, are important values for understanding the properties of these dark solitons.

In future research, we hope to find more precise criteria for determining the behavior of dark solitons. For instance, we aim to obtain conditions on the interaction potential to ensure that the corresponding dark solitons exhibit either monotonic or oscillating behavior.

*Acknowledgments.* The authors acknowledge support from the Labex CEMPI (ANR-11-LABX-0007-01). A. de Laire was also supported by the ANR project ODA (ANR-18-CE40-0020-01). S. López-Martínez was supported by the Madrid Government (Comunidad de Madrid – Spain) under the multiannual Agreement with UAM in the line for the Excellence of the University Research Staff in the context of the V PRICIT (Regional Program of Research and Technological Innovation).
<figure><img src="image_29.png" /><figcaption>FIGURE 28. Numerically computed solitons for potential (E5) with $\lambda = 1.0$, showing $\eta_k$ (left) and $\theta_k$ (right) as a function of $x_k$.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>FIGURE 29. Numerically computed solitons for potential (E5) with $\lambda = 1.0$, showing $P^{\delta x}$ (left) and $E^{\delta x}$ (center) as a function of $c$. On the right, $E^{\delta x}$ as a function of $P^{\delta x}$.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>FIGURE 30. Numerically computed solitons for potential (E5) with $\lambda = 4.0$, showing $\eta_k$ (left) and $\theta_k$ (right) as a function of $x_k$.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>FIGURE 31. Numerically computed solitons for potential (E5) with $\lambda = 4.0$, showing $P^{\delta x}$ (left) and $E^{\delta x}$ (center) as a function of $c$. On the right, $E^{\delta x}$ as a function of $P^{\delta x}$.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>FIGURE 32. Numerically computed solitons for potential (E6), showing $\eta_k$ (left) and $\theta_k$ (right) as a function of $x_k$.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>FIGURE 33. Numerically computed solitons for potential (E6), showing $P^{\delta x}$ (left) and $E^{\delta x}$ (center) as a function of $c$. On the right, $E^{\delta x}$ as a function of $P^{\delta x}$.</figcaption></figure>
APPENDIX

Proof of Theorem 2. Let $m = \inf_{\xi \in \mathbb{R}} (\tilde{W}(\xi)^{-})$, where $s^{-} = \min\{s, 0\}$, and let $R = \sqrt{2(\sigma - m)}$. For any $\tilde{\sigma} \in (0, \sigma)$, we choose $k = (\tilde{\sigma} - m)/R^2 = (\tilde{\sigma} - m)/(2(\sigma - m))$. Clearly, $k \in (0, 1/2)$. On the one hand, for a.e. $|\xi| \ge R$, one trivially has

$$ \tilde{W}(\xi) \ge m = \tilde{\sigma} - kR^2 \ge \tilde{\sigma} - k\xi^2. $$

On the other hand, observe that $\sigma - \xi^2/2 \ge \tilde{\sigma} - k\xi^2$ if, and only if, $|\xi| \le R$. Therefore, for a.e. $|\xi| \le R$, we get

$$ \tilde{W}(\xi) \ge \sigma - \frac{\xi^2}{2} \ge \tilde{\sigma} - k\xi^2. $$

We may now apply Theorem 1.1 in [10] and get a solution to $(TW_c)$ for a.e. $c \in (0, \sqrt{2\tilde{\sigma}})$. Taking $\tilde{\sigma}$ arbitrarily close to $\sigma$ we obtain a solution for a.e. $c \in (0, \sqrt{2\sigma})$, which completes the proof. □

REFERENCES

[1] A. Aftalion, X. Blanc, and R. L. Jerrard. Mathematical issues in the modelling of supersolids. *Nonlinearity*, 22(7):1589–1614, 2009.

[2] C. F. Barenghi, R. J. Donnelly, and W. Vinen. *Quantized vortex dynamics and superfluid turbulence*, volume 571. Springer Science & Business Media, 2001.

[3] C. Becker, S. Stellmer, P. Soltan-Panahi, S. Dörscher, M. Baumert, E.-M. Richter, J. Kronjäger, K. Bongs, and K. Sengstock. Oscillations and interactions of dark and dark–bright solitons in Bose–Einstein condensates. *Nature Physics*, 4(6):496–501, 2008.

[4] N. G. Berloff and P. H. Roberts. Motions in a Bose condensate VI. Vortices in a nonlocal model. *J. Phys. A*, 32(30):5611–5625, 1999.

[5] F. Béthuel, P. Gravejat, and J.-C. Saut. Existence and properties of travelling waves for the Gross-Pitaevskii equation. In *Stationary and time dependent Gross-Pitaevskii equations*, volume 473 of *Contemp. Math.*, pages 55–103. Amer. Math. Soc., Providence, RI, 2008.

[6] F. Béthuel, P. Gravejat, J.-C. Saut, and D. Smets. Orbital stability of the black soliton for the Gross-Pitaevskii equation. *Indiana Univ. Math. J.*, 57(6):2611–2642, 2008.

[7] F. Bethuel, P. Gravejat, and D. Smets. Asymptotic stability in the energy space for dark solitons of the Gross-Pitaevskii equation. *Ann. Sci. Éc. Norm. Supér. (4)*, 48(6):1327–1381, 2015.

[8] D. Chiron. Stability and instability for subsonic traveling waves of the nonlinear Schrödinger equation in dimension one. *Anal. PDE*, 6(6):1327–1420, 2013.

[9] A. de Laire. Global well-posedness for a nonlocal Gross-Pitaevskii equation with non-zero condition at infinity. *Comm. Partial Differential Equations*, 35(11):2021–2058, 2010.

[10] A. de Laire and S. López-Martínez. Existence and decay of traveling waves for the nonlocal Gross-Pitaevskii equation. *Comm. Partial Differential Equations*, 47(9):1732–1794, 2022.

[11] A. de Laire and P. Mennuni. Traveling waves for some nonlocal 1D Gross-Pitaevskii equations with nonzero conditions at infinity. *Discrete Contin. Dyn. Syst.*, 40(1):635–682, 2020.

[12] D. J. Frantzeskakis. Dark solitons in atomic Bose–Einstein condensates: from theory to experiments. *Journal of Physics A: Mathematical and Theoretical*, 43(21):213001, 2010.

[13] M. Grillakis, J. Shatah, and W. Strauss. Stability theory of solitary waves in the presence of symmetry. *I. J. Funct. Anal.*, 74(1):160–197, 1987.

[14] Q. Kong, Q. Wang, O. Bang, and W. Krolikowski. Analytical theory of dark nonlocal solitons. *Opt. Lett.*, 35(13):2152–2154, 2010.

[15] W. Królikowski and O. Bang. Solitons in nonlocal nonlinear media: Exact solutions. *Physical Review E*, 63(1):016610, 2000.

[16] Z. Lin. Stability and instability of traveling solitonic bubbles. *Adv. Differential Equations*, 7(8):897–918, 2002.

[17] S. Lopez-Aguayo, A. S. Desyatnikov, Y. S. Kivshar, S. Skupin, W. Krolikowski, and O. Bang. Stable rotating dipole solitons in nonlocal optical media. *Opt. Lett.*, 31(8):1100–1102, 2006.
[18] N. I. Nikolov, D. Neshev, W. Królikowski, O. Bang, J. J. Rasmussen, and P. L. Christiansen. Attraction of nonlocal dark optical solitons. *Optics letters*, 29(3):286-288, 2004.

[19] J. Reneuve, J. Salort, and L. Chevillard. Structure, dynamics, and reconnection of vortices in a nonlocal model of superfluids. *Physical Review Fluids*, 3(11):114602, 2018.

[20] H. Veksler, S. Fishman, and W. Ketterle. Simple model for interactions and corrections to the Gross-Pitaevskii equation. *Phys. Rev. A*, 90(2):023620, 2014.