††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Topologically nontrivial $1/3$-magnetization plateau state in a spin-1/2
trimer chain

Y. Y. Han1 B. C. Yu3 Z. Du1 L. S. Ling1 L. Zhang1<EMAIL_ADDRESS>W. Tong1
C. Y. Xi1 J. L. Zhang1 T. Shang3,4<EMAIL_ADDRESS>Li Pi1,2 Long Ma1
<EMAIL_ADDRESS>1Anhui Key Laboratory of Low-Energy Quantum Materials and
Devices, High Magnetic Field Laboratory, Chinese Academy of Sciences, Hefei
230031, China
2 Hefei National Research Center for Physical Sciences at the Microscale,
University of Science and Technology of China, Hefei 230026, China
3 Key Laboratory of Polar Materials and Devices (MOE), School of Physics and
Electronic Science, East China Normal University, Shanghai 200241, China
4 Chongqing Key Laboratory of Precision Optics, Chongqing Institute of East
China Normal University, Chongqing 401120, China

###### Abstract

Topologically nontrivial Haldane phase is theoretically proposed to be
realized in the 1/3-magnetization ($M$) plateau of spin-1/2 trimer systems.
However, the spin excitation gap, typical characteristic of Haldane phase, is
not yet experimentally verified. Here, we report the nuclear magnetic
resonance investigations into the low-energy spin dynamics in the $S=1/2$
spin-trimer antiferromagnetic chain compound Na2Cu3Ge4-xSixO12 ($x=0,0.1\sim
1.5$). In the parent compound ($x=0$), the spin-lattice relaxation rate
(1/$T_{1}$) shows significantly different temperature dependence when the
external magnetic field is increased above the critical field of
$\mu_{0}$$H_{c}$ = 29 T. The spin excitation gap is evidenced from the
thermally activated behavior of $1/T_{1}(T)$ in the 1/3-$M$ plateau state. By
substituting Ge4+ with Si4+, the critical field for the 1/3-$M$ plateau
significantly decreases, e.g. $\mu_{0}H_{c}=17$ T in $x=1.0$ samples, which
results from the suppressed inter-trimer coupling $J_{2}$. The gapped spin
excitation is confirmed again above 17 T, whose size shows temperature-
dependent behavior for $\mu_{0}H\geq 25.72$ T. These observations provide
further insights into the Haldane physics.

Topology has become an important concept in physics since the discovery of
quantum spin Hall effect and topological insulatorsHasan and Kane (2010). It
supplies a new perspective on states of matter beyond the local order
parameters in conventional Ginzburg-Landau theory. From the application
standpoint, as a robust physical property against disturbance of circumstance,
topological nontrivial states serve as a robust physical resistant to the
circumstance disturbance, supplying the hardware basis for quantum computation
with low error ratesBartlett et al. (2010). In the fruitful field of
topological insulator, the research starts from the band theory based on
single particle approximation. While, in strongly-correlated electron systems,
exploring topological nontrivial states is still challengingWang et al.
(2014); Chen et al. (2011); Pollmann et al. (2012).

Due to the enhanced spin fluctuations as well as reduced dimensionality,
quantum spin chains provide a valuable playground for studying exotic quantum
phases, transitions, and emergent spin excitations. As indicated by the Bethe
ansatz solutionBethe (1931), the ground state of one-dimensional(1D) spin-1/2
Heisenberg spin chain is the quantum critical Tomonaga-Luttinger liquid. The
spin excitation spectrum is dominated by the gapless spinon continuum with a
fractional quantum number. While for the integer spin case, a finite spin
excitation gap exists, as conjectured by HaldaneHaldane (1983a, b).
Subsequently, Affleck, Kennedy, Lieb, and Tasaki (AKLT) propose a model
Hamiltonian for spin-1 by artificially introducing the quadratic term of spin
exchange interactionsAffleck et al. (1987), which can be solved exactly. A
novel topological string order and edge states have been found to exist in
this AKLT stateAffleck et al. (1987), which is believed to capture the nature
of Haldane phase. Further theoretical work reveals that Haldane chains with
odd integer spin host typical symmetry protected topological orderGu and Wen
(2009); Pollmann et al. (2012). The Haldane gap and edge states are verified
by extensive experimental study of nickel-based materials with spin-1 magnetic
sites.

The Haldane phase also can be realized in the magnetization ($M$) plateau
state of spin-1/2 trimer chains. Theoretically, the precondition for the
existence of plateau in the $M(H)$ curve for quantum spin chains is the
Oshikawa-Yamanaka-Affleck ruleOshikawa et al. (1997),$S\cdot
n(1-M/M_{\rm{sat}})=\rm{integer}$. Here, $S$ and $n$ denote the spin value of
the magnetic ion and number of magnetic ions per unit cell. The
$M/M_{\rm{sat}}$ represents the ratio between the magnetization and the
saturated magnetization. The 1/3-$M$ plateau state of spin-1/2 trimmer
antiferromagnetic chains is proposed theoretically to be topological
nontrivial with a finite Haldane gapHu et al. (2014).

In the spin trimmer chain compound Cu3(P2O6OH)2 ($J_{1}-J_{2}-J_{2}$), a
1/3-$M$ plateau was observed above the critical field of 12 THase et al.
(2006). A singlet-triplet-like excitation gap responsible for the plateau was
observed from the inelastic neutron scattering under zero magnetic fieldHase
et al. (2007). Recently, emergent composite excitations of doublons and
quartons have been reported in another spin-trimmer antiferromagnet
Na2Cu3Ge4O12, where the 1/3-$M$ plateau appears above $\mu_{0}H_{c}\sim 28$
TBera et al. (2022). However, to the best of our knowledge, the experimental
evidence supporting the existence of a finite Haldane gap in the magnetization
plateau state is still missing in the spin-1/2 trimmer chains, clearly
requires further investigations.

In this paper, we use nuclear magnetic resonance (NMR) technique to explore
the low-energy spin excitations of the magnetization plateau state in
Na2Cu3Ge4-xSixO12 ($x=0,0.1\sim 1.5$). In the parent compound ($x=0$), the
incommensurate magnetic ordering is evidenced by the spectral intensity
modulation below $T_{N}=2$ K in the powder samples aligned by magnetic field.
The 1/3-$M$ plateau state appears above $\mu_{0}H_{c}=29$ T, confirmed by
high-field magnetization measurements. A thermally activated behavior is
observed from the temperature dependent nuclear spin-lattice relaxation rates.
This clearly indicates the gap opening in the spin excitation spectrum, which
is believed to be the topologically nontrivial Haldane gap. We further tune
the magnetic coupling between the spin trimmers via replacing Ge4+ with
smaller Si4+. Upon doping, the magnetic transition is quickly suppressed and
the critical field above which the spin gap emerges shifts to the lower field
side. This effect is most likely attributed to the reduced inter-trimer
magnetic coupling. For $x=1.0$ samples under magnetic fields ranging from
25.72 T to 32.5 T, the spin gap size obtained from the spin-lattice relaxation
rates varies for different temperature region, which can be attributed to the
dominant spin excitations at different wave vector $\bf{q}$ in the reciprocal
space.

Polycrystalline samples of Na2Cu3Ge4-xSixO12 were synthesized using a solid-
state reaction methodBera et al. (2022). The NMR measurements were performed
on 23Na nuclei ($\gamma_{n}=11.262$ MHz/T, $I=3/2$) with a phase-coherent NMR
spectrometer. The spectrum was obtained by summing up the spectral intensity
with the frequency swept. The spin-lattice relaxation rate ($1/T_{1}$) was
measured with the inversion-recovery sequence, and obtained by fitting the
time dependence of nuclear magnetization to ($M(t)$) to
$M(t)/M(\infty)=1-0.1\exp[-(t/T_{1})^{\beta}]-0.9\exp[-(6t/T_{1})^{\beta}]$.

Figure 1: (color online) (a) 23Na NMR spectra at different temperatures under
a 10 T magnetic field. The powder sample was finely grounded, loosely put into
a capsule and aligned with the applied NMR magnetic field. The red lines are
simulations to the spectra. (b) Temperature dependence of 23Na Knight shifts
under 10 T and $dc$-susceptibility. The $dc$-susceptibility was measured under
a weak magnetic field of 0.1 T. The crystal structure of Na2Cu3Ge4O12 is
sketched in the inset.

The Na2Cu3Ge4-xSixO12 samples crystallize in a triclinic structure with the
$P\text{-}1$ space groupBera et al. (2022), as shown in the inset of Fig.1(b).
Three edge-sharing CuO4 planes form Cu3O8 spin trimers with a intra-trimer
coupling $J_{1}$ (shown by yellow lines). The Cu3O8-trimers interact with each
other via O-Ge-O super-superexchange interactions, whose coupling strength is
denoted as $J_{2}$ (shown by red lines). For each trimer, the Cu2+ on both
sides couple via oxygen-oxygen ions, resulting in the magnetic coupling of
$J_{3}$. The magnetic transition at $T_{N}=2$ K results from the weak inter-
chain coupling $J_{4}$, bridged by the Na+ and Ge4+ ions.

Microscopic information about the magnetic order can be obtained from the NMR
spectrum of the aligned sample. Above $T_{N}$, 23Na spectrum is composed by
one central peak and double satellite peaks. We simulate the spectrum at
$T=140$ K with three lorentz peaks and one gaussian background counting for
the deviation from complete alignment. This is consistent with the single
Wyckoff position of Na+ ions carrying nuclear spin $I=3/2$. For temperatures
below $T_{N}$, each peak broadens and exhibits a distinctive double-horn
shape, expected for the incommensurate magnetic orderZeng et al. (2022); Higa
et al. (2017). The resonance line shape for the incommensurate modulation is
simulated as shown by the red line at the bottom of Fig.1(a), consistent with
the previous NMR studyYasui et al. (2014). The Knight shifts, measuring the
spin susceptibility, exhibit a highly similar behavior to the temperature-
dependent $dc$-susceptibility (See Fig.1(b)). The broad maximum at $T^{*}\sim
11$ K indicate the setup of short-range order upon cooling, which is commonly
observed in low-dimensional magnetsBonner and Fisher (1964).

Figure 2: (color online) (a) Field dependence of the magnetization at $T=1.8$
K and $T=4$ K in Na2Cu3Ge4O12. The sketched spin configurations of the low
magnetic state and plateau state are shown as insets. (b) Spin-lattice
relaxation rates as a function of temperature under magnetic fields ranged
from 10 T to 32.85 T. Inset: Semi-log plot of the $1/T_{1}$ versus the inverse
of temperature. Both the red lines in the main figure and inset denote the
thermally activated behavior of $1/T_{1}(T)$.

Fig.2(a) presents the $dc$-magnetization versus external magnetic field
measured at $T=1.8$ K and 4 K. The magnetization first shows a linear field
dependence at $\mu_{0}H<22$ T, then increases rapidly up to $\mu_{0}H_{c}=29$
T, followed by a level-off behavior under higher field. The slope change point
is defined as the critical field for the 1/3-$M$ plateau state, which is close
to the $\mu_{0}H_{c}=28$ T obtained by previous magnetization measurements
under pulsed magnetic field. The corresponding spin states in different field
range are depict in the inset of Fig.2(a). The circled dumbbells represent the
spin singlets formed by two adjacent Cu2+. The moments of the remaining Cu2+
ion in the spin trimmer is denoted by black arrows, which are polarized in the
plateau.

To verify the existence of spin gap in the magnetization plateau state, we
measured the temperature dependence of the spin-lattice relaxation rate under
various magnetic fields (see Fig.2(b)). For $\mu_{0}H=10$ T, the $1/T_{1}$
first decreases with cooling the sample, then increases at $T_{N}<T<5$ K, and
finally drops steeply below $T_{N}$. The peak centered at $T_{N}$ results from
the critical slowing down of electron spin dynamics near the transition
temperature. At $\mu_{0}H=27.92$ T just below $\mu_{0}H_{c}$, the upturn
behavior persists down to the lowest temperature, indicating strong spin
fluctuations near the magnetic instability. When the field intensity increases
to 30.3 T, just 1.3 T higher than $\mu_{0}H_{c}$, a thermally activated
behavior dominates the low temperature $1/T_{1}(T)$ instead, in sharp contrast
with that below $\mu_{0}H_{c}$. For higher $\mu_{0}H=32.85$ T, the
$1/T_{1}(T)$ drops more steeply upon cooling, demonstrating a larger spin
excitation gap. In the inset, we present the semi-log plot of $1/T_{1}$ versus
inverse temperature. The gapped behavior is shown more clearly by the straight
lines, whose slopes are directly proportional to the gap size.

Figure 3: (color online) (a) The magnetization at $T=1.8$ K as a function of
field intensity in Na2Cu3Ge4-xSixO12 with $x=$0, 0.5, 1, and 1.5. The red
lines represent the magnetization measured by the VSM option of PPMS below
$\mu_{0}H=14$ T. (b) Temperature dependence of the 23Na spin-lattice
relaxation rate under $\mu_{0}H=10$ T in samples with $x=$0, 0.1, 0.2, 0.5,
0.7, 1, and 1.5.

The critical field for the magnetization plateau state in spin trimmer systems
is closely related to the ratio of inter-trimer to intra-trimer magnetic
coupling strengthBera et al. (2022). We take a step forward to tune the
magnetic coupling and thus, the critical field via chemical doping. The
interaction between spin trimmers is bridged by the oxygen, germanium and
oxygen. By substituting the Ge with Si, we expect effective modification to
the inter-trimer coupling $J_{2}$. As shown in Fig.3(a), low temperature
field-dependent magnetization was measured for Na2Cu3Ge4-xSixO12 with various
Si contents. The predominant characteristic of the $M(H)$ curve is maintained
upon doping except two facts. First, for higher doping levels, the
magnetization below 20 T enhances more rapidly with increasing magnetic field.
The level-off trend appears at much lower critical fields upon doping. Second,
for $x=1$ sample, the $M(H)$ curve bends over at a moderate field
$\mu_{0}H_{c}=17$ T before its flattening at $\mu_{0}H=25$ T, which are marked
by the up and down arrows. The level-off behavior becomes gradual in the
$x=1.5$ sample, which makes the determination of $H_{c}$ difficult. The
decreasing $H_{c}$ in doped samples is consistent with the reduced $J_{2}$
according to the theoretical phase diagramBera et al. (2022).

The temperature dependence of spin-lattice relaxation rates in doped samples
is measured under a constant magnetic field of 10 T (see Fig.3(b)). The
relaxations of 23Na nuclei overall get faster with the doping level
increasing. The $1/T_{1}$ shares a similar linear temperature dependence at
$T>60$ K. When the sample is further cooled down, a low-temperature broad
maximum in $1/T_{1}(T)$ gradually emerges as the doping level increases. For
$x=0\sim 0.2$, the sharp peak, corresponding to the critical slowing down of
spin fluctuations near the magnetic transition, is suppressed to a lower
temperature of $T_{N}=1.67$ K for $x=0.1$, and finally disappears in the
$x=0.2$ sample for the present temperature range.

The linear temperature dependence of $1/T_{1}$ at $T\geq 60$ K is contributed
by the uniform spin susceptibility as indicated by the scaling theory study of
$S=1/2$ Heisenberg antiferromagnetic spin chain systemsSachdev (1994). This
relation is verified experimentally in various $S=1/2$ spin chains, e.g.
Sr2CuO3Thurber et al. (2001). The suppression of the magnetic transition
indicates the fragility of the magnetic order against chemical substitution in
Na2Cu3Ge4O12, which may relate with the reduced inter-chain coupling $J_{4}$
in the doped samples. As indicated by the low temperature $M(H)$ curve shown
in Fig.3(a), no obvious change of slope is observed near $\mu_{0}H=10$ T,
where the magnetization is far from 1/3 of the saturated one. Thus, the spin
state at $\mu_{0}H=10$ T cannot be regarded as that in the plateau state.
Based on this consideration, the low temperature broad maximum in the samples
of $x=1.0$ and $1.5$ is most likely due to the emergence of certain type of
short-range order with typical timescale close to the inverse Larmor frequency
($\sim 112.5$ MHz).

Figure 4: (color online) (a) Spin-lattice relaxation rates as a function of
temperature in the $x=1.0$ sample with the field intensity ranging from 10 T
to 32.5 T. (b) The semi-log plot of $1/T_{1}$ versus the inverse temperature.
(c) Field dependence of the spin excitation gap size determined by
$1/T_{1}(T)$ in the $x=0$ and $1.0$ samples.

The suppressed critical field for the 1/3-$M$ plateau state in doped samples
provides a valuable opportunity to track the field dependence of the spin
excitation gap. In Fig.4(a), we show the temperature dependence of 23Na spin-
lattice relaxation rates under various magnetic fields up to 32.5 T. The low
temperature $1/T_{1}(T)$ drops steeply upon cooling when the field intensity
exceeds 17 T. Under higher fields, even more pronounced drop in $1/T_{1}(T)$
was observed at low temperatures. A kinked characteristic was observed for
fields larger than 25.72 T as demonstrated by the thermally activated function
shown by the red lines in Fig.4(a). The gapped behavior can be seen more
clearly by the straight lines in the semi-log plot of $1/T_{1}$ versus the
inverse temperature (see Fig.4(b)). The kinks reflect the change of gap size,
which equals to the slopes of straight lines. The field dependence of the spin
excitation gap as determined from $1/T_{1}(T)$ is shown in Fig.4(c) for the
parent sample and sample with $x=1.0$. The $\Delta_{r1}$ and $\Delta_{r2}$
denote the gap size determined from the $1/T_{1}(T)$ of the $x=1.0$ sample at
intermediate and low temperatures, respectively. For the field range in this
study, the gap size increases rapidly with magnetic field strength. The
$\Delta_{r1}/\Delta_{r2}$ equals to 1.54, 2.04, 1.86, and 1.2 for
$\mu_{0}H=25.72$ T, 28 T, 30.4 T and 32.5 T, respectively.

Information about the spin dynamics of the magnetization plateau state can be
inferred from nuclear relaxation measurements. The spin-lattice relaxation
results from the hyperfine coupling between a nuclear spin and electron spins.
The $1/T_{1}$ can be expressed as a weighted sum over momenta
$\overrightarrow{q}$ of the spin-correlation function
$S(\overrightarrow{q},\omega_{0})$Moriya (1963), i.e.,
$1/T_{1}=\frac{\gamma_{N}^{2}\hbar}{2\mu_{B}^{2}}\sum_{\overrightarrow{q}}\lvert
A_{\perp}(\overrightarrow{q})\rvert^{2}S(\overrightarrow{q},\omega_{0})$.
Here, $A_{\perp}(\overrightarrow{q})$ denotes the momentum-dependent hyperfine
coupling constant in the plane perpendicular to the applied magnetic field.
The $\omega_{0}$ is NMR frequency, generally treated as zero in the studied
correlated spin system.

For the Haldane spin chains, the $1/T_{1}(T)$ at low temperatures is expected
to be dominated by a thermally activated behavior with the function form,
i.e., $1/T_{1}\propto\exp(-\gamma\Delta_{H}/T)$, where $\Delta_{H}$ is the
Haldane gap size, and $\gamma$ is the ratio between the gap $\Delta_{r}$
determined from $1/T_{1}(T)$ and $\Delta_{H}$. The $\gamma$-value varies among
different Haldane spin systems. Experimentally, the temperature dependence of
$1/T_{1}$ gives $\gamma=1$ in the typical Haldane spin system
Ni(C2H8N2)2NO2ClO4 (NENP)Gaveau et al. (1990). However, $\gamma$ is determined
to be $\sim 1.2$ in Y2BaNiO5 and $\sim 1.5$ in AgVP2S6, respectivelyShimizu et
al. (1995); Takigawa et al. (1995, 1996). In the $S=1$ spin chain compound
SrNi2V2O8, the $1/T_{1}(T)$ in the temperature range of 15 K $<T<$ 80 K and 4
K $<T<$ 10 K give very different $\gamma$-valuesPahari et al. (2006). The
behavior sharply contrasts with the temperature dependence of Knight shifts,
which measures the spin susceptibility at $\overrightarrow{q}=0$, implying
prominent $\overrightarrow{q}$-dependence of the low-energy spin dynamics.
Theoretically, the $1/T_{1}$ can be calculated by assuming a uniform hyperfine
coupling constant in the reciprocal spaceJolicoeur and Golinelli (1994); Sagi
and Affleck (1996); Sachdev and Damle (1997); Konik (2003). Controversies
about the $\gamma$-value still exists depending on the different theoretical
approaches.

Recently, the NMR relaxation rate in the Haldane spin chains has been
calculated using numerical techniquesCapponi et al. (2019), by taking into
count the contributions from modes with momenta close to $q\approx 0$ and
$q\approx\pi$. For the low temperature regime ($T<0.41\Delta_{H}$), the
$1/T_{1}$ is mainly contributed by two-magnon Raman process at $q=0$ as a
result of $\omega_{0}\ll\Delta_{H}$ and the energy conservation law. The
$1/T_{1}(T)$ follows a simple thermally activated behavior, i.e.,
$1/T_{1}\propto\exp(-\Delta_{H}/T)$. For $0.41\Delta_{H}<T<\Delta_{H}$, the
spin excitations at $q=\pi$ dominate nuclear relaxations, leading to a steeper
drop of $1/T_{1}$ upon cooling,where the $\gamma$-value is larger than 1.5.

As shown in Fig.4, the temperature dependence of $1/T_{1}$ in the $x=1.0$
sample yields distinct spin gap sizes denoted as $\Delta_{r1}$ and
$\Delta_{r2}$ for the intermediate and low temperature regime. We understand
this phenomenon with the numerical results for $S=1$ Haldane chains as
illustrated above. The variation of $\Delta_{r1}/\Delta_{r2}$ is attributed to
the fine tuning of the weighted momentum contribution to $1/T_{1}$ by applied
fields.

For both $x=0$ and $x=1.0$ samples, the spin excitation gap size of the
1/3-$M$ plateau state increases with the increasing magnetic field as depicted
in Fig.4(c). For typical $S=1$ Haldane chains, the spin excitation gap
decreases with applying magnetic field, and finally closes above the critical
field, which results from Zeeman splitting of the three-fold degenerated
excited stateFujiwara et al. (1992, 1993). In this study, the magnetic field
as high as 32.5 T enables us to access the 1/3-$M$ plateau state. However, it
is insufficient to fully suppress the gapless spinon excitations observed at
zero fieldBera et al. (2022). When the spin-trimer system enters the plateau
state from low magnetic field side, the gradually suppressed gapless spinon
excitation leads to the increasing gap size observed here. We expect that the
spin excitation gap size will decrease under higher magnetic fields.

For isotropic $S=1$ Heisenberg spin chains, the Haldane gap $\Delta_{H}$
equals to $0.41J$ ($J$ is the exchange interaction constant)Nightingale and
Blote (1986). In Na2Cu3Ge4O12, the intra-trimer coupling $J_{1}$ is determined
to be $\sim 235$ K from numerical simulations of the magnetization dataBera et
al. (2022). The Cu3O8 spin trimers formed by edge-sharing CuO4 planes remain
unchanged upon chemical doping at Ge4+-sites. Thus, the Haldane gap for both
$x=0$ and $x=1.0$ samples is estimated to be $\sim 96$ K. However, the upper
limit for the Haldane gap observed from $1/T_{1}(T)$ under a 32.5 T field is
$\sim 26$ K, about one quarter of the calculated value. There are two possible
reasons for this discrepancy. One is the fore mentioned contribution of the
gapless spinon excitations to $1/T_{1}$, which prevents the determination of
intrinsic Haldane gap size. The other one is the existence of single ion
anisotropy, which reduces the Haldane gap size, and is not considered in the
proposed spin model.

In summary, we have performed high-field NMR study on the low-energy spin
dynamics in the $S=1/2$ spin-trimer chain antiferromagnet Na2Cu3Ge4-xSixO12
($x=0\sim 1.5$). A non-zero spin gap is verified from the temperature
dependence of $1/T_{1}$ for the 1/3-$M$ plateau state in the parent compound.
By suppressing the inter-trimer coupling via silicon doping, the critical
field decreases to 17 T for samples with $x=1.0$, which is accessible for
commercial superconducting magnets. The gapped spin excitation is confirmed
again above 17 T, whose size shows temperature-dependent behavior for
$\mu_{0}H\geq 25.72$ T. The temperature dependent gap size can be understood
within the numerical calculations of $1/T_{1}$ contributed by different
momenta. The variation of $\Delta_{r1}/\Delta_{r2}$ under different fields is
most likely attributed to the fine tuning of the weighted momentum
contribution to $1/T_{1}$ by applied fields. These observations provide direct
experimental evidence for the existence of topologically nontrivial Haldane
phase in the 1/3-$M$ plateau of $S=1/2$ spin-trimer chain antiferromagnets.
The doped samples with lower $H_{c}$ provide a better platform to perform
further experimental investigations on the 1/3-$M$ plateau state.

This research was supported by the National Natural Science Foundation of
China (Grants No. 11874057, 12374105 and 12374128), the Natural Science
Foundation of Shanghai (Grants No. 21ZR1420500 and 21JC1402300), Natural
Science Foundation of Chongqing (Grant No. CSTB-2022NSCQ-MSX1678),the Basic
Research Program of the Chinese Academy of Sciences Based on Major Scientific
Infrastructures (Grant No. JZHKYPT-2021-08) and Fundamental Research Funds for
the Central Universities. A portion of this work was supported by the High
Magnetic Field Laboratory of Anhui Province.

## References

  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Bartlett et al. (2010) S. D. Bartlett, G. K. Brennen, A. Miyake, and J. M. Renes, Phys. Rev. Lett. 105, 110502 (2010).
  * Wang et al. (2014) C. Wang, A. C. Potter, and T. Senthil, Science 343, 629 (2014).
  * Chen et al. (2011) X. Chen, Z.-C. Gu, and X.-G. Wen, Phys. Rev. B 83, 035107 (2011).
  * Pollmann et al. (2012) F. Pollmann, E. Berg, A. M. Turner, and M. Oshikawa, Phys. Rev. B 85, 075125 (2012).
  * Bethe (1931) H. Bethe, Z. Phys. 71, 205 (1931).
  * Haldane (1983a) F. Haldane, Phys. Lett. A 93, 464 (1983a).
  * Haldane (1983b) F. D. M. Haldane, Phys. Rev. Lett. 50, 1153 (1983b).
  * Affleck et al. (1987) I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Phys. Rev. Lett. 59, 799 (1987).
  * Gu and Wen (2009) Z. C. Gu and X. G. Wen, Phys. Rev. B 80, 155131 (2009).
  * Oshikawa et al. (1997) M. Oshikawa, M. Yamanaka, and I. Affleck, Phys. Rev. Lett. 78, 1984 (1997).
  * Hu et al. (2014) H. P. Hu, C. Cheng, Z. H. Xu, H. G. Luo, and S. Chen, Phys. Rev. B 90, 035150 (2014).
  * Hase et al. (2006) M. Hase, M. Kohno, H. Kitazawa, N. Tsujii, O. Suzuki, K. Ozawa, G. Kido, M. Imai, and X. Hu, Phys. Rev. B 73, 104419 (2006).
  * Hase et al. (2007) M. Hase, M. Matsuda, K. Kakurai, K. Ozawa, H. Kitazawa, N. Tsujii, A. Donni, M. Kohno, and X. Hu, Phys. Rev. B 76, 064431 (2007).
  * Bera et al. (2022) A. K. Bera, S. M. Yusuf, S. K. Saha, M. Kumar, D. Voneshen, Y. Skourski, and S. A. Zvyagin, Nat. Comm. 13, 6888 (2022).
  * Zeng et al. (2022) K. Y. Zeng, F. Y. Song, L. S. Ling, W. Tong, S. L. Li, Z. M. Tian, L. Ma, and L. Pi, Chin. Phys. Lett. 39, 107501 (2022).
  * Higa et al. (2017) N. Higa, Q. P. Ding, M. Yogi, N. S. Sangeetha, M. Hedo, T. Nakama, Y. Onuki, D. C. Johnston, and Y. Furukawa, Phys. Rev. B 96, 024405 (2017).
  * Yasui et al. (2014) Y. Yasui, Y. Kawamura, Y. Kobayashi, and M. Sato, J. Appl. Phys. 115, 17E125 (2014).
  * Bonner and Fisher (1964) J. C. Bonner and M. E. Fisher, Phys. Rev. 135, A640 (1964).
  * Sachdev (1994) S. Sachdev, Phys. Rev. B 50, 13006 (1994).
  * Thurber et al. (2001) K. R. Thurber, A. W. Hunt, T. Imai, and F. C. Chou, Phys. Rev. Lett. 87, 247202 (2001).
  * Moriya (1963) T. Moriya, J. Phys. Soc. Jpn. 18, 516 (1963).
  * Gaveau et al. (1990) P. Gaveau, J. P. Boucher, L. P. Regnault, and J. P. Renard, Europhys. Lett. 12, 647 (1990).
  * Shimizu et al. (1995) T. Shimizu, D. E. MacLaughlin, P. C. Hammel, J. D. Thompson, and S.-W. Cheong, Phys. Rev. B 52, R9835 (1995).
  * Takigawa et al. (1995) M. Takigawa, T. Asano, Y. Ajiro, and M. Mekata, Phys. Rev. B 52, R13087 (1995).
  * Takigawa et al. (1996) M. Takigawa, T. Asano, Y. Ajiro, M. Mekata, and Y. J. Uemura, Phys. Rev. Lett. 76, 2173 (1996).
  * Pahari et al. (2006) B. Pahari, K. Ghoshray, R. Sarkar, B. Bandyopadhyay, and A. Ghoshray, Phys. Rev. B 73, 012407 (2006).
  * Jolicoeur and Golinelli (1994) T. Jolicoeur and O. Golinelli, Phys. Rev. B 50, 9265 (1994).
  * Sagi and Affleck (1996) J. Sagi and I. Affleck, Phys. Rev. B 53, 9188 (1996).
  * Sachdev and Damle (1997) S. Sachdev and K. Damle, Phys. Rev. Lett. 78, 943 (1997).
  * Konik (2003) R. M. Konik, Phys. Rev. B 68, 104435 (2003).
  * Capponi et al. (2019) S. Capponi, M. Dupont, A. W. Sandvik, and P. Sengupta, Phys. Rev. B 100, 094411 (2019).
  * Fujiwara et al. (1992) N. Fujiwara, T. Goto, S. Maegawa, and T. Kohmoto, Phys. Rev. B 45, 7837 (1992).
  * Fujiwara et al. (1993) N. Fujiwara, T. Goto, S. Maegawa, and T. Kohmoto, Phys. Rev. B 47, 11860 (1993).
  * Nightingale and Blote (1986) M. P. Nightingale and H. W. J. Blote, Phys. Rev. B 33, 659 (1986).


# Number of solitons produced from a large initial pulse in the generalized
NLS dispersive hydrodynamics theory

L. F. Calazans de Brito  Higher School of Economics, Moscow, Russia  A. M.
Kamchatnov  Higher School of Economics, Moscow, Russia  Institute of
Spectroscopy, Russian Academy of Sciences, Troitsk, Moscow, 108840, Russia

###### Abstract

We show that the number of solitons produced from an arbitrary initial pulse
of the simple wave type can be calculated analytically if its evolution is
governed by a generalized nonlinear Schrödinger equation provided this number
is large enough. The final result generalizes the asymptotic formula derived
for completely integrable nonlinear wave equations like the standard NLS
equation with the use of the inverse scattering transform method.

###### pacs:

05.45.Yv, 42.65.Tg, 47.35.Fg

##  I Introduction

It is well known that if a nonlinear wave system supports propagation of
solitons with some polarity (bright or dark ones), then a large enough initial
pulse with the same polarity evolves at asymptotically large times to a
sequence of $N$ separate solitons with some amount of linear radiation. The
number $N$ of soliton depends on the initial data and for $N\gg 1$ the energy
contained in linear radiation is negligibly small compared with the solitons’
energy in the main approximation with respect to the small parameter $1/N\ll
1$. Consequently, the number of solitons $N$ is one of the most important
characteristics of the initial pulse. Besides that, it is relatively easy to
measure $N$ experimentally. Thus, the possibility to predict the number of
solitons from the initial data is an important task in the theory of nonlinear
waves.

If dynamics of the system under consideration is described by a completely
integrable wave equation, then the number of solitons is equal to the number
of eigenvalues of the associated linear spectral problem (see, e.g., [1, 2]).
Consequently, for large number $N\gg 1$ one can apply the quasiclassical WKB
method to the linear spectral problem and obtain the asymptotic formula for
$N$. For example, in case of the famous KdV equation

$u_{t}+6uu_{x}+u_{xxx}=0$ (1)

such a formula was first obtained by Karpman [3] and it reads

$N\approx\frac{1}{\pi}\int_{-\infty}^{\infty}\sqrt{u_{0}(x)}\,dx,\qquad N\gg
1,$ (2)

where $u_{0}(x)$ is the initial distribution of the amplitude $u$ and $N$ is
the number of eigenvalues of the Schrödinger spectral problem [4, 1, 2] with
the “potential” $u_{0}(x)$. A similar formula was derived in Ref. [5] for the
NLS equation

$i\psi_{t}+\frac{1}{2}\psi_{xx}-|\psi|^{2}\psi=0$ (3)

associated with the Zakharov-Shabat spectral problem [6]. This formula is
conveniently formulated in terms of the initial data for the hydrodynamics-
like form of Eq. (3) obtained by means of the Madelung substitution

$\psi(x,t)=\sqrt{\rho(x,t)}\exp(i\phi(x,t)),\qquad\phi_{x}=u(x,t),$ (4)

so that the NLS equation reduces to the system

$\begin{split}&\rho_{t}+(\rho u)_{x}=0,\\\
&u_{t}+uu_{x}+\rho_{x}+\left(\frac{\rho_{x}^{2}}{8\rho^{2}}-\frac{\rho_{xx}}{4\rho}\right)_{x}=0.\end{split}$
(5)

In the context of the Gross-Pitaevskii theory [7, 8] of Bose-Einstein
condensates of diluted gases, $\rho$ has the meaning of the gas density and
$u$ is its flow velocity. This interpretation becomes especially clear in the
limit of large pulses with their characteristic size $l\gg 1$. Then the space
derivative has the order of magnitude $\partial_{x}\sim l^{-1}\ll 1$ and the
last term in the second equation (5) can be neglected. Hence we arrive at
equations of Euler hydrodynamics

$\rho_{t}+(\rho u)_{x}=0,\qquad u_{t}+uu_{x}+\rho_{x}=0$ (6)

for a compressible fluid with the equation of state $P=\frac{1}{2}\rho^{2}$,
where $P$ plays the role of “pressure”.

Dynamics of such a gas is conveniently described with the use of the variables

$r^{\pm}=\frac{1}{2}u\pm\sqrt{\rho}$ (7)

called Riemann invariants, so that Eqs. (6) acquire a diagonal form

$\begin{split}&\frac{\partial r^{+}}{\partial
t}+\frac{1}{2}(3r^{+}+r^{-})\frac{\partial r^{+}}{\partial x}=0,\\\
&\frac{\partial r^{-}}{\partial t}+\frac{1}{2}(r^{+}+3r^{-})\frac{\partial
r^{-}}{\partial x}=0.\end{split}$ (8)

The Gross-Pitaevskii equation has dark soliton solutions [9] propagating in
the form of dips along a uniform background
$\rho=\overline{\rho},u=\overline{u}$. If the initial distributions
$\overline{\rho}_{0}(x),\overline{u}_{0}(x)$ correspond to a large dip in the
density, $0<\overline{\rho}_{0}(x)<\overline{\rho}$ and
$\overline{u}_{0}(x)\to 0$ as $|x|\to\infty$, then we can transform them to
distributions of the Riemann invariants
$r_{0}^{\pm}(x)=\overline{u}_{0}(x)/2\pm\sqrt{\overline{\rho}_{0}(x)}$. The
asymptotic formulas for a number of dark solitons moving in the positive or
negative direction are given, correspondingly, by the equations [5]

$N_{\pm}=\frac{1}{\pi}\int_{-\infty}^{\infty}\sqrt{(\pm\overline{r}-r_{0}^{+}(x))(\pm\overline{r}-r_{0}^{-}(x))}\,dx,$
(9)

where $\overline{r}=\sqrt{\overline{\rho}}$.

Less rigorous simple approach to this problem was suggested in Ref. [10] for
the whole AKNS hierarchy [11] associated with $2\times 2$-matrix spectral
problem written in an equivalent scalar form [12]. Naturally, in case of Eq.
(3) it reproduces Eq. (9) (see Ref. [13]). The approach of Refs. [10, 13] is
based on the supposition that the solution of the linear spectral problem (the
Baker-Akhiezer function) corresponding to the periodic nonlinear wave keeps
its form for slightly modulated waves and can be formally continued to the
initial state where it becomes a quasiclassical eigenfunction of the spectral
problem.

For not completely integrable equations the associated linear spectral problem
does not exist and the above approach becomes impossible. Since the process of
formation of solitons from an initially smooth pulse is universal from
physical point of view and, generally speaking, it does not depend on the fact
whether the wave equation is completely integrable or not, other ideas are
necessary. An alternative approach was suggested in Refs. [14, 15] and it was
based on the Gurevich-Pitaevskii theory of dispersive shock waves (DSWs).
According to this theory, transformation of an initially smooth pulse goes
through several stages. In the first stage, the pulse changes its form
remaining smooth up to the moment of wave breaking when the dispersionless
approximation breaks down due to appearance of infinite space derivatives
(“gradient catastrophe”). After the wave breaking moment, the second stage
starts during which the pulse in the simplest case consists of two parts—a
smooth part for which the dispersionless approximation remains correct, and
the DSW part represented by a modulated periodic solution of the nonlinear
wave equation under consideration where the modulation parameters obey the
Whitham modulation equations [17] (see also review articles [18, 19] for more
details). In Whitham approximation, the boundary between the two parts is
sharp and it is called the small-amplitude edge of the DSW. At last, in the
third asymptotic stage of the pulse’s evolution its smooth part becomes
negligibly small and the whole pulse evolves to a sequence of separate
solitons. The method of Refs. [14, 15] was based on the assumption that the
solution of the Whitham equations at the small-amplitude edge of the DSW can
be formally prolonged to the smooth part of the pulse. The necessary
parameters of the DSW at its small-amplitude edge can be obtained in the
important case of initial simple waves with one of the Riemann invariants
constant by El’s method [20]. Although the results of this approach were
confirmed in several particular cases by their comparison with the results of
numerical solutions, it hardly can be considered as obvious, because the
Whitham equations are obtained by averaging over fast nonlinear wave
oscillations, whereas there are no any oscillations in the smooth part of the
pulse and introduction of their wave vector here looks quite artificial.
Therefore other approach based on clearer physical ideas seems very desirable.

Recently a new approach has been suggested in Refs. [19, 21]. It is based on
an old remark of Gurevich and Pitaevskii [22] that at the mentioned above
second stage of evolution the number $N(t)$ of oscillations (wave crests)
inside the DSW part increases with time according to the equation

$\frac{dN}{dt}=\frac{1}{2\pi}k(v_{g}-v_{ph}),$ (10)

where $k$ is the wave number and $v_{g},v_{ph}$ are the group and the phase
velocities, correspondingly, of the wave at the small-amplitude edge of the
DSW. If we denote the background amplitude at this edge as $u$ then the
dependence $k=k(u)$ can be found by El’s method and the dependence $t=t(u)$
along the path of the small-amplitude edge can be found by the method of Ref.
[23]. As a result, Eq. (10) can be integrated along the path of the small-
amplitude edge and in the limit $t\to\infty$ we obtain the total number of
crests in the DSW which evolve eventually to separate solitons. As was shown
in Refs. [21, 24] for several simple particular cases, this method yields the
formulas for the number of solitons coinciding with those derived by the
method of Refs. [14, 15], so in these cases the formulas are justified by a
more reliable approach.

In this paper, we consider the problem of calculation of the number of
solitons produced from a simple-wave type of an initial pulse which evolves
according to the generalized NLS (gNLS) equation

$i\psi_{t}+\frac{1}{2}\psi_{xx}-\frac{1}{p}|\psi|^{2p}\psi=0.$ (11)

It reduces to the standard NLS equation (3) for $p=1$, but for $p\neq 1$ it is
not completely integrable. Various forms of nonlinearity in the gNLS equation
are used in nonlinear physics. In particular, the regime $2/3<p<1$ describes
the superfluid BEC–Bardeen–Cooper–Schrieffer transition in ultracold Fermi
gases [25] and the case $p=2/3$ of the so-called unitary limit has drawn
considerable attention (see, e.g., [26] and reference therein). Besides that,
Eq. (11) with different values of $p$ are very useful and convenient for the
study of various properties of dispersive Euler hydrodynamics [27] since it
models both integrable and non-integrable situations amenable to analytical
study. Here we show that the formula for the number of solitons can be derived
by the direct method of Refs. [19, 21] and this formula can be used in
applications and general nonlinear physics investigations.

##  II Basic formulas

We shall present here the basic formulas necessary for derivation of the
expression for the number $N$ of solitons.

###  II.1 Simple waves

We assume that at the dispersionless stage of evolution before the wave
breaking moment the pulse has a form of a simple wave. Actually, this is not a
serious restriction imposed on the form of the pulse, since for quite
arbitrary initial conditions the dispersionless evolution leads in a natural
way to splitting the pulse to two pulses propagating in opposite directions
and they both are of the simple wave type. Besides that, simple wave pulses
can be created by a proper choice of the initial conditions as it happens, for
example, in the flow caused by a unidirectional motion of a piston (see, e.g.,
[28]).

Equations of the dispersionless approximation for Eq. (11) can be obtained by
means of substitution (4) into it, separation of real and imaginary parts and
neglecting the dispersive terms with higher order space derivatives. As a
result, we obtain the Euler equations

$\rho_{t}+(\rho u)_{x}=0,\qquad u_{t}+uu_{x}+\rho^{p-1}\rho_{x}=0$ (12)

corresponding to the equation of state $P=\rho^{p+1}/(p+1)$. It is convenient
to define the “sound velocity” variable

$c=\sqrt{dP/d\rho}=\rho^{p/2}.$ (13)

Then for the Riemann invariants

$r_{\pm}=\frac{u}{2}\pm\frac{c}{p}$ (14)

Eqs. (12) transform to

$\partial_{t}r_{\pm}+(u\pm c)\partial_{x}r_{\pm}=0,$ (15)

where the characteristic velocities $v_{\pm}=u\pm c$ have simple physical
meaning—they correspond to waves propagating with the sound velocity
$c=c(\rho)$ downstream or upstream, correspondingly.

Figure 1:  (a) The initial distribution of the local sound velocity
$\overline{c}(x)$. (b) The inverse function $\overline{x}(c)$ has two branches
$x_{1}(c)$ and $x_{2}(c)$.

In a simple wave one of the Riemann invariants is constant and, to be
definite, we assume that $r_{-}=u/2-c/p=-\overline{c}_{0}/p$, where
$\overline{c}_{0}$ is the sound velocity far from the localized initial pulse.
Then the flow velocity $u(x,t)$ is expressed in terms of the distribution
$c(x,t)$ of the local sound velocity,

$u(x,t)=\frac{2}{p}(c(x,t)-\overline{c}_{0}),$ (16)

hence $r_{+}=(2c(x,t)-\overline{c}_{0})/p$ and Eq. (15) for $r_{+}$ is
transformed to the equation

$\frac{\partial c}{\partial
t}+\frac{1}{p}[(2+p)c-2\overline{c}_{0}]\frac{\partial c}{\partial x}=0.$ (17)

This is the well known Hopf equation (see, e.g., [29]) and its solution reads

$x-\frac{1}{p}[(2+p)c-2\overline{c}_{0}]t=\overline{x}(c),$ (18)

where $\overline{x}(c)$ is the function inverse to the initial distribution
$c=\overline{c}(x)$ of the sound velocity (13). We are interested in the
initial distributions in the form of localized dips in the uniform state with
$c=\overline{c}_{0}$ (see Fig. 1(a)). Hence, the inverse function has two
branches $x=\overline{x}_{1,2}(c)$ (see Fig. 1(b)), and each branch gives the
solution (18) which defines in implicit form the function $c=c(x,t)$.
Evolution of the profile $c=c(x,t)$ leads to wave breaking and, to simplify
the notation, we assume that the wave breaking moment corresponds to $t=0$.

###  II.2 Dispersion relation and El’s equation

After the wave breaking moment a DSW is generated and its small amplitude edge
starts its motion along the smooth part of the pulse. Actually, this edge is
represented by a linear wave packet propagating with the group velocity of
linear waves. The corresponding dispersion relation can be easily found after
linearization of Eq. (11) and it can be written in the form (see, e.g., [27])

$\omega(k,c)=k\left(u\pm\sqrt{c^{2}+k^{2}/4}\right).$ (19)

The background distribution $c=c(x,t)$ is changing according to the Hopf
equation (17) which should be combined with the Hamilton equations

$\frac{dx}{dt}=\frac{\partial\omega}{\partial
k},\qquad\frac{dk}{dt}=-\frac{\partial\omega}{\partial x}$ (20)

for the packet’s motion. As a simple consequence of these three equations we
easily obtain [21] the equation

$\frac{dk}{dc}=\frac{(\partial\omega/\partial
c)_{k}}{u+c-(\partial\omega/\partial k)_{c}}.$ (21)

This equation was first obtained by El [20] from the small-amplitude limit of
the Whitham equations. Its derivation from Hamilton’s equations demonstrates
its much wider applicability region. For example, it was applied in Ref. [30]
to propagation of localized wave packets along arbitrary large scale smooth
simple waves. In case of motion of the small-amplitude edge of a DSW it should
be solved with the initial condition

$k(\overline{c}_{0})=0$ (22)

which means that the wave breaking occurs at the boundary with the undisturbed
region where $c=\overline{c}_{0}$ and in the Whitham approximation the DSW at
the initial wave breaking moment shrinks to a point without any oscillations
inside.

To solve this equation, it is convenient to introduce the variable [14]

$\alpha(c)=\sqrt{1+\frac{k^{2}(c)}{4c^{2}}},\qquad
k(c)=2c\sqrt{\alpha^{2}(c)-1},$ (23)

so that Eq. (21) transforms to

$\frac{d\alpha}{dc}=-\frac{(1+\alpha)(2/p-1+2\alpha)}{c(1+2\alpha)},$ (24)

which can be easily solved with the initial condition
$\alpha(\overline{c}_{0})=1$ to give [27]

$\frac{c(\alpha)}{\overline{c}_{0}}=\left(\frac{2}{1+\alpha}\right)^{\frac{p}{3p-2}}\left(\frac{2+p}{2-p+2p\alpha}\right)^{\frac{2(p-1)}{3p-2}}.$
(25)

The function $c=c(\alpha)$ can be inverted in three particular cases [27]:

$\begin{split}&p=1,\quad\alpha(c)=\frac{2\overline{c}_{0}}{c}-1,\\\
&p=2,\quad\alpha(c)=\frac{1}{2}\left(\sqrt{1+8\left(\frac{\overline{c}_{0}}{c}\right)^{2}}-1\right),\\\
&p=\frac{1}{2},\quad\alpha(c)=\frac{1}{16}\left(5\sqrt{\frac{\overline{c}_{0}}{c}\left(25\frac{\overline{c}_{0}}{c}-16\right)}+25\frac{\overline{c}_{0}}{c}-24\right).\end{split}$
(26)

Substitution of these formulas into Eq. (23) yields explicit expressions for
the function $k(c)$.

###  II.3 Path of the small amplitude edge

The small amplitude edge propagates with the group velocity

$\frac{dx}{dt}=\frac{d\omega}{dk}=\frac{2}{p}(c-\overline{c}_{0})+c(\alpha(c)-\alpha^{-1}(c)).$
(27)

According to Ref. [23], this equation must be compatible with the solution
(18), where $c$ equals to the local value of the sound velocity at the point
of location of the wave packet. This condition leads to the linear equation

$\frac{c(\alpha-1)(2\alpha+1)}{\alpha}\frac{dt}{dc}-\frac{2+p}{p}t=\overline{x}^{\prime}(c),$
(28)

which should be solved first with the initial condition
$t(\overline{c}_{0})=0$. Easy integration yields

$\begin{split}t=t_{1}(c)&=\frac{1}{\sqrt{c(\alpha(c)-1)}}\left(\frac{\alpha(c)+1}{2-p+2p\alpha(c)}\right)^{\frac{1}{3p-2}}\\\
&\times\int_{\overline{c}_{0}}^{c}dc^{\prime}\frac{\overline{x}_{1}(c^{\prime})}{\sqrt{c^{\prime}(\alpha(c^{\prime})-1)}(2\alpha(c^{\prime})+1)}\\\
&\times\left(\frac{2-p+2p\alpha(c^{\prime})}{\alpha(c^{\prime})+1}\right)^{\frac{1}{3p-2}}\end{split}$
(29)

for the motion of the wave packet along the branch $x_{1}(c)$: at the moment
$t$ the small-amplitude edge is located at the point with the local value of
the sound velocity $c$. The expression (29) is correct up to the moment
$t_{m}=t_{1}(c_{m})$ when the edge reaches the point with the minimal value
$c_{m}$ of the local sound velocity. It is worth noticing that the
dispersionless solution (18) preserves the minimal value $c_{m}$ in the
distribution $c=c(x,t)$, so $c_{m}$ does not depend on $t$. For $t>t_{m}$ the
edge propagates along the solution (18) corresponding to the second branch
$\overline{x}_{2}(c)$. Accordingly, Eq. (28) should be solved with the initial
condition $t(c_{m})=t_{m}$ and the solution reads

$\begin{split}t=t_{2}(c)&=\frac{1}{\sqrt{c(\alpha(c)-1)}}\left(\frac{\alpha(c)+1}{2-p+2p\alpha(c)}\right)^{\frac{1}{3p-2}}\\\
&\times\Bigg{\\{}\int_{\overline{c}_{0}}^{c_{m}}dc^{\prime}\frac{\overline{x}_{1}(c^{\prime})}{\sqrt{c^{\prime}(\alpha(c^{\prime})-1)}(2\alpha(c^{\prime})+1)}\\\
&\times\left(\frac{2-p+2p\alpha(c^{\prime})}{\alpha(c^{\prime})+1}\right)^{\frac{1}{3p-2}}\\\
&+\int_{c_{m}}^{c}dc^{\prime}\frac{\overline{x}_{2}(c^{\prime})}{\sqrt{c^{\prime}(\alpha(c^{\prime})-1)}(2\alpha(c^{\prime})+1)}\\\
&\times\left(\frac{2-p+2p\alpha(c^{\prime})}{\alpha(c^{\prime})+1}\right)^{\frac{1}{3p-2}}\Bigg{\\}}.\end{split}$
(30)

Substitution of the above two formulas into
$x(c)=\overline{x}_{1,2}(c)+\frac{1}{p}[(2+p)c-2\overline{c}_{0}]t_{1,2}(c)$
gives the coordinate of the edge at the moment $t(c)$, but we do not need
these expressions here.

##  III Number of solitons

Integration of formula (10) with the use of the dispersion relation (19) gives
the total number of wave crests containing in the DSW,

$N=\frac{1}{8\pi}\int_{0}^{\infty}\frac{k^{3}}{\sqrt{c^{2}+k^{2}/4}}\,dt,$
(31)

so our task is to evaluate this integral. To this end we pass to integration
over $c$ and substitute (23) for $k=k(c)$,

$N=\frac{1}{\pi}\int\frac{c^{2}}{\alpha(c)}(\alpha^{2}(c)-1)^{3/2}\frac{dt}{dc}dc.$
(32)

In fact, this integral consists of two parts: first with integration from
$\overline{c}_{0}$ to $c_{m}$ with $t=t_{1}(c)$ and second with integration
from $c_{m}$ to $\overline{c}_{0}$ with $t=t_{2}(c)$. The derivative $dt/dc$
can be excluded with the use of Eq. (28) and as a result we obtain

$\begin{split}N=&\frac{1}{\pi}\int_{\overline{c}_{0}}^{c_{m}}\frac{c(\alpha-1)^{1/2}(\alpha+1)^{3/2}}{2\alpha+1}\\\
&\times\left[\frac{2+p}{p}(t_{2}-t_{1})+\overline{x}_{2}-\overline{x}_{1}\right]dc.\end{split}$
(33)

Integration over $c$ can be replaced by integration over $\alpha$ with help of
Eq. (24),

$\begin{split}N=&\frac{p}{\pi}\int_{1}^{\alpha_{m}}\frac{c^{2}(\alpha)(\alpha^{2}-1)^{1/2}}{2-p+2p\alpha}\\\
&\times\left[\frac{2+p}{p}(t_{2}-t_{1})+\Phi_{2}-\Phi_{1}\right]d\alpha,\end{split}$
(34)

where $\alpha_{m}=\alpha(c_{m})$ and
$\Phi(\alpha)\equiv\left.\overline{x}^{\prime}(c)\right|_{c=c(\alpha)}$.

When we substitute here the expressions (29) and (30) for $t_{1}$ and $t_{2}$
as well as Eq. (25) for $c(\alpha)$, we obtain the double integral which can
be simplified by integration by parts:

$\begin{split}&\frac{2+p}{\pi}\int_{1}^{\alpha_{m}}\frac{c^{2}(\alpha)(\alpha^{2}-1)^{1/2}}{2-p+2p\alpha}(t_{2}-t_{1})d\alpha\\\
&=\frac{1}{\pi}\int_{c_{m}}^{\overline{c}_{0}}\frac{c\alpha\sqrt{\alpha^{2}-1}(\overline{x}^{\prime}_{2}-\overline{x}^{\prime}_{1})}{2\alpha+1}dc.\end{split}$
(35)

The remaining terms are transformed to

$\begin{split}&\frac{p}{\pi}\int_{1}^{\alpha_{m}}\frac{c^{2}(\alpha)(\alpha^{2}-1)^{1/2}}{2-p+2p\alpha}(\Phi_{2}-\Phi_{1})d\alpha\\\
&=\frac{1}{\pi}\int_{c_{m}}^{\overline{c}_{0}}\frac{c(\alpha+1)\sqrt{\alpha^{2}-1}(\overline{x}^{\prime}_{2}-\overline{x}^{\prime}_{1})}{2\alpha+1}dc,\end{split}$
(36)

and the sum of the above two expressions yields

$N=\frac{1}{\pi}\int_{c_{m}}^{\overline{c}_{0}}c(\alpha^{2}-1)^{1/2}(\overline{x}^{\prime}_{2}-\overline{x}^{\prime}_{1})dc.$
(37)

At last, we replace integration over the interval $c_{m}\leq
c\leq\overline{c}_{0}$ by integration over $x$ and take into account Eq. (23)
to obtain the final formula

$N=\frac{1}{2\pi}\int_{-\infty}^{\infty}k(\overline{c}(x))dx,$ (38)

where $\overline{c}(x)$ is the initial distribution of the local sound
velocity. In this expression the function $k(c)$ is obtained by substitution
of the function $\alpha(c)$ into Eq. (23). For example, in case $p=1$ of the
standard NLS equation we have for $\alpha(c)$ the first expression in Eq.
(30), hence $k(c)=4\sqrt{\overline{c}_{0}(\overline{c}_{0}-c)}$ and Eq. (38)
reduces to

$N=\frac{2}{\pi}\int_{-\infty}^{\infty}\sqrt{\overline{c}_{0}(\overline{c}_{0}-\overline{c}(x))}dx.$
(39)

This formula coincides with formula (9) for $N_{+}$, since for the simple wave
initial conditions we have
$r_{-}=-\overline{c}_{0},r_{+}=2\overline{c}(x)-\overline{c}_{0}$.

We have obtained Eq. (38) by direct calculation without artificial
introduction of the distribution of wave numbers $k(\overline{c}(x))$ along a
smooth initial state, as it was done in Refs. [14, 15]. Thus, our derivation
makes the assumption about existence of such a distribution quite plausible
(see [21]).

##  IV Comparison with numerical solutions

Figure 2:  The number of solitons produced from the initial pulse (40) for
$l=20$ and different values of $a$. Solid lines correspond to (a) $p=1/2$; (b)
$p=1$; and (c) $p=2$ calculated according to the analytical formula (38).

Here we support our analytical theory by comparison of its predictions with
the results of numerical solutions of the gNLS equation. We take the initial
distribution of the local sound velocity $c=\rho^{p/2}$ in the form

$\overline{c}(x)=\left[1-\frac{a}{\cosh^{2}(x/l)}\right]^{p/2},$ (40)

where $a$ denotes the depth of the initial dip in the density distribution and
$l$ is its half-length. The initial disturbance must be a simple wave in our
approach, so we define the initial distribution of the flow velocity by the
formula

$\overline{u}_{0}(x)=\frac{2}{p}(\overline{c}(x)-\overline{c}_{0}).$ (41)

Equations (4) allow us to find the initial field $\psi(x,0)$.

In our numerical experiments we choose $l=20$ and change $a$ in the interval
$0.1\leq a\leq 0.9$ with the step $\Delta a=0.1$. Such a choice makes the
numerical calculations not too time consuming and clearly demonstrates the
dependence of the number of solitons on $a$. We have done these calculations
for $p=1/2,1$ and 2, when the function $k(\overline{c}(x))$ is given by the
explicit formulas (23) and (26). The results of our calculations are presented
in Fig. 2. As one can see, the agreement is very good.

##  V Conclusion

We have shown that the method of calculation of the number of solitons
produced from an initial pulse of the simple wave type works very well for the
generalized NLS equation having various physical application. The resulting
formula (38) has the structure suggested earlier in Ref. [14, 15] on the basis
of some suppositions about the properties of solutions of Whitham modulation
equations and our derivation makes these suppositions quite plausible. Thus,
the developed method and the general formula (38) become a useful tool for
predictions of the number of solitons in experiments performed with media
whose evolution is described by nonlinear wave equations not belonging to a
specific class of completely integrable equations.

###### Acknowledgements.

The work was supported by the Foundation for the Advancement of Theoretical
Physics and Mathematics “BASIS”.

## References

  * [1] V. E. Zakharov, S. V. Manakov, S. P. Novikov, and L. P. Pitaevskii, The Theory of Solitons: The Inverse Scattering Method, (Nauka, Moscow, 1980) (translation: Consultants Bureau, 1984).
  * [2] A. C. Newell, Solitons in Mathematics and Physics, (SIAM, Philadelphia, 1985).
  * [3] V. I. Karpman, Phys. Lett. A 25, 708 (1967).
  * [4] C. S. Gardner, J. M. Greene, M. D. Kruskal, and R. M. Miura, Phys. Rev. Lett. 19, 1095 (1967).
  * [5] S. Jin, C. D. Levermore, D. W. McLaughlin, Comm. Pure Appl. Math., 52, 613 (1999).
  * [6] V. E. Zakharov, A. B. Shabat, Zh. Eksp. Teor. Fiz., 64, 1627 (1973) [Sov. Phys. JETP, 37, 823 (1973)].
  * [7] E. P. Gross, Nuovo Cimento, 20, 454 (1961).
  * [8] L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 40, 646 (1961) [Sov. Phys. JETP, 13, 451 (1961)].
  * [9] T. Tsuzuki, J. Low Temp. Phys. 4, 441 (1971).
  * [10] A. M. Kamchatnov, R. A. Kraenkel, B. A. Umarov, Phys. Lett. A 287, 223 (2001).
  * [11] M. J. Ablowitz, D. J. Kaup, A. C. Newell, H. Segur, Stud. Appl. Math., 53, 249-315 (1974).
  * [12] A. M. Kamchatnov and R. A. Kraenkel, J. Phys. A: Math. Gen. 35, L13 (2002).
  * [13] A. M. Kamchatnov, R. A. Kraenkel, B. A. Umarov, Phys. Rev. E 66, 036609 (2002).
  * [14] G. A. El, A. Gammal, E. G. Khamis, R. A. Kraenkel, and A. M. Kamchatnov, Phys. Rev. A 76, 053813 (2007).
  * [15] G. A. El, R. H. J. Grimshaw, N. F. Smyth, Physica D 237, 2423 (2008).
  * [16] A. V. Gurevich, L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 65, 590 (1973) [Sov. Phys. JETP, 38, 291 (1974)].
  * [17] G. B. Whitham, Proc. R. Soc. Lond. A 283, 238 (1965).
  * [18] G. A. El, M. A. Hoefer, Physica D 333, 11 (2016).
  * [19] A. M. Kamchatnov, Usp. Fiz. Nauk, 191, 52 (2021) [Physics-Uspekhi, 64, 48 (2021)].
  * [20] G. A. El, Chaos, 15, 037103 (2005).
  * [21] A. M. Kamchatnov, Chaos, 30, 123148 (2020).
  * [22] A. V. Gurevich, L. P. Pitaevskii, Zh. Eksp. Teor. Fiz. 93, 871 (1987) [Sov. Phys. JETP, 66, 490 (1987)].
  * [23] A. M. Kamchatnov, Phys. Rev. E 99, 012203 (2019).
  * [24] A. M. Kamchatnov, Zh. Eksp. Teor. Fiz., 159, 76 (2021) [JETP, 132, 63 (2021)].
  * [25] W. Ketterle, M. W. Zwierlein, Riv. Nuovo Cimento, 31, 247 (2008).
  * [26] E. A. Kuznetsov, M. Yu. Kagan, and A. V. Turlapov, Phys. Rev. A 101, 043612 (2020).
  * [27] M. A. Hoefer, J. Nonlinear Sci. 24, 525 (2014).
  * [28] L. D. Landau and E. M. Lifshitz, Fluid Mechanics, (Pergamon, Oxford, 1959).
  * [29] G. B. Whitham, Linear and Nonlinear Waves, (Wiley Interscience, New York, 1974).
  * [30] A. M. Kamchatnov and D. V. Shaykin, Phys. Fluids, 33, 052120 (2021).