text:
# 1 Introduction

## 1 Introduction

Some localized solutions of a special (and, actually, quite a narrow) class of nonlinear evolution equations are known to be able to move conserving their shape and interacting with each other like particles. The typical examples are the solitons of Nonlinear Schrödinger equation (NLS) and Korteveg-de-Vriese equation (KdV), the kinks of sine-Gordon equation (SG) to name a few, ,. It is worthy of mentioning that these particular solutions conserve all the intergals of motions prescribed by dynamic equations. However, as the evolution equations mentioned above are the result of approximation of more complicated dynamic equations, the solitonic properties of their solutions are the approximation as well. The question up to which extent they survive in ‘reality’ leads to interesting physical problems that urge people to go beyond the theory of exactly solvable models ,,.

Fortunately, the exact results obtained within refine mathematical models can be compared with the real physical experiments. In addition, numerical simulations within more realistic models are performed for many equations supporting solitonic solutions . As a result, one can realize under which conditions the concept of soliton is a good approximation to interpret physical data. For example, the NLS equation is used to describe the phenomenon of Bose-Einstein condensation of dilute trapped gases  and sharp pulses of light in optical fibers , and the KdV equation is used in hydrodynamics to model the moving packets of surface waves in channels . The SG equation models surface grow and reconstruction in material science .

In this article, we discuss the transport properties of excited states (generally of electron transition nature) in a medium (a periodic structure, such as a crystal or semiconductor structure, some periodic biological tissues, etc.). Under certain conditions, a coherent state, which involves both the electronic (or spin) excitations and the medium elastic excitations, can appear in such structures. This state turns out to be localized inside the medium and can travel ballistically through it. In Biological Physics, such a state is known as Davydov soliton , whereas in Condensed Matter Physics, one can mention, for example, Spin-Peierls Systems , Charge Density Waves , and anomalous transport of excitonic packets in semiconducting crystals and heterostructures . In particular, we are motivated by the ballistic transport of excitons in pure 3D crystals, such as Cu20 . However, the ballistic transport of localized coherent structures and a set of nonlinear equations being used to describe it can be found in modeling of different physical phenomena , , . For example, modeling of the coherent excited states of electronic origin ends up with NLS equation, whereas the coherent states of the elastic medium are described by the wave equation with nonlinear anharmonic terms (NLW equation) ,,. In plasma physics, this is the case of Zakharov system of equations ,,, while the Davey-Stewartson system can be mentioned in the context of fluid mechanics . Note that the two equations comprising such a system, NLS and NLW, are quite different from the mathematical point of view, but they are coupled with each other ,, and this is the origin of the wealth of physical effects they can describe.

Not surprisingly for the moving soliton-like states in periodic media, a shape of experimentally registered signals is far from being look like a “true” soliton (e.g., $`1/\mathrm{cosh}\left(L_0^1(xvt)\right)`$ ) or kink. One can assume that nonzero temperature, scattering on impurities and just a noise factor can be responsible for this. In some cases, however, the character of changement hints on some regular dynamic reason for these effects. For example, an asymmetric form of the signal with the pronounced sharp front and a long lasting tail behind the single soliton are observed on experiment ,,. Note that the subsonic ballistic transport is under consideration, i.e., $`v<c_s`$, where $`v`$ is the speed of the packet and $`c_s`$ is the (longitudinal) sound speed in the periodic medium. The similar effects, such as delocalized solitons, were also observed in different physical settings. For example, in a kind of pump and probe experiments with sharp pulses of light in fibers, the solitonic shape of a probe pulse was significantly changed by Stimulated Raman Scattering ; see also . In fluid mechanics, one can mention the problem of wave breaking and appearance of the so-called peakons and compactons in the soliton-like solutions of generalized KdV equations ; see also . In essence, the non-steady solitons are found to be able to move with an acceleration and generate “tails” behind themselves.

In contrast, the transition to diffusion regime is possible for the ballistic exciton-phonon droplets, and, in the case of amplification of the soliton-like state formed inside the exciton-phonon packet, its ballistic velocity is found to be almost unchanged during such a non-stationary process . These facts hint at different physics controlling the moving coherent state formed by excitons and phonons.

In this article, we assume that the coupling between collective excitations of the periodic medium can be responsible for such effects in Condensed Matter Physics, and a proper dynamic model can reproduce them even at $`T0`$. To support this hypothesis, we consider a field model that admits a soliton-like solution and can be easily generalized without a loss of physical clarity.

This is a two field model, in which the excitations of the medium are modeled by a nonideal Bose-gas coupled with long wavelength modes of the displacement field. For example, two (or even several) interacting boson fields can model different phenomena in Condensed Matter Physics . In the case of semiconductors, one can often disregard the influence of free fermions or fermionic complexes on the processes under consideration. It turns out that many processes involving photons, excitons, and phonons can be described by use of the language of interacting Bose-fields . Moreover, if there is a branch of optically inactive excitons in a crystal, one can even exclude the photons from simple models dealing with such excitons. In addition, the lifetime of a moving exciton with $`\mathrm{}k_0m_\mathrm{x}c_s`$ can be large enough in semiconductor structures, so that the transport properties of a packet of moving excitons can be observed experimentally ,.

## 2 Hamiltonian of the model

The Hamiltonian of the medium can be taken in the following form:

$$\widehat{H}=H_{\mathrm{gas}}(\widehat{\psi },\widehat{\psi }^{})+H_{\mathrm{ph}}(\widehat{u},\widehat{\pi })+H_{\mathrm{int}}(\widehat{u},\widehat{\psi }^{}\widehat{\psi }).$$
(1)
Here, $`\widehat{\psi }`$, $`\widehat{\psi }^{}`$ are the Bose-gas field operators; they stands in the non-relativistic Hamiltonian $`H_{\mathrm{gas}}`$, whereas $`\widehat{u}`$ is the displacement field operator and $`\widehat{\pi }`$ is the momentum density operator conjugate to $`\widehat{u}`$. This pair stands in a Phonon Hamiltonian $`H_{\mathrm{ph}}`$. We consider a nonideal gas of phonons in the long-wavelength approximation and take into account the acoustic branch only, e.g., the medium is a 3D crystal of the volume $`V=LS_{}`$.

The Bose-gas Hamiltonian has the following form:

$$H_{\mathrm{gas}}=𝑑𝐱\stackrel{~}{E}_g\widehat{\psi }^{}\widehat{\psi }(𝐱)+\frac{\mathrm{}^2}{2m}\widehat{\psi }^{}\widehat{\psi }(𝐱)+\frac{\nu _0}{2}\left(\widehat{\psi }^{}(𝐱)\right)^2\left(\widehat{\psi }(𝐱)\right)^2+\frac{\nu _1}{3}\left(\widehat{\psi }^{}(𝐱)\right)^3\left(\widehat{\psi }(𝐱)\right)^3,$$

where $`\nu _0>0`$ is the strength of two particle repulsive interaction, and $`\nu _1>0`$ is the strength of three particle one. As the two particle interaction ($`\nu _0(\psi ^{})^2\psi ^2`$ ) can be strongly renormalized because of interaction with other fields, we include the hard core interaction term (modeled by repulsion $`\nu _1(\psi ^{})^3\psi ^3`$ in $`H_{\mathrm{gas}}`$). We assume that the ‘bare’ characteristic energies of the particle-particle interactions satisfy the following inequality:

$$0<\nu _1/a_\mathrm{x}^6<()\nu _0/a_\mathrm{x}^3\mathrm{const}\mathrm{Ry}_\mathrm{x},\mathrm{const}10,$$

see  for discussion. Here, $`a_\mathrm{x}`$ and $`\mathrm{Ry}_\mathrm{x}`$ are the exciton Bohr radius and characteristic Rydberg energy, respectively. We count the energy of a free particle from $`\stackrel{~}{E}_g=E_g\mathrm{Ry}_\mathrm{x}>0`$, so that $`E_𝐤=\stackrel{~}{E}_g+\mathrm{}𝐤^2/2m`$. (For semiconducting materials, $`E_g`$ is the fundamental gap.)

As only the longitudinal phonons $`𝐮_l`$ interact with the Bose-field $`\widehat{\psi }`$, $`\widehat{\psi }^{}`$ in our model, one can exclude the transversal phonons ($`𝐮_t=0`$) from $`H_{\mathrm{ph}}`$. Then, it can be reduced to the following simple form:

$$H_{\mathrm{ph}}=\frac{\widehat{\pi }^2(𝐱)}{2\rho }+\frac{\rho c_l^2}{2}_j\widehat{u}_s_j\widehat{u}_s(𝐱)+\frac{\rho c_l^2}{3}\kappa _3_j\widehat{u}_s_j\widehat{u}_s_j\widehat{u}_s(𝐱)+\frac{\rho c_l^2}{4}\kappa _4\left(_j\widehat{u}_s\right)^4d𝐱,$$

where $`c_l`$ is the longitudinal sound speed of the crystal, and the dimensionless parameters $`\kappa _3`$ and $`\kappa _4`$ account for cubic and quartic nonlinearities. (In this article, we will not take into account a quartic term, so $`\kappa _40`$.)

We take the gas-phonon interaction in the form of Deformation Potential:

$$H_{\mathrm{int}}=\sigma _0_j\widehat{u}_j(𝐱)\widehat{\psi }^{}\widehat{\psi }(𝐱)+\vartheta _0_j\widehat{u}_j(𝐱)\left(\widehat{\psi }^{}\mathrm{\Delta }\widehat{\psi }(𝐱)\right)d𝐱,$$
(2)
where $`\sigma _0`$ and $`\vartheta _0`$ are the coupling constants. Note that this Hamiltonian is equivalent to

$$\stackrel{~}{H}_{\mathrm{int}}=\sigma _0_j\widehat{u}_j(𝐱)\widehat{\psi }^{}\widehat{\psi }(𝐱)+\vartheta _0_j\widehat{u}_j(𝐱)\widehat{\psi }^{}\widehat{\psi }(𝐱)d𝐱.$$
(3)
Here, we choose $`\sigma _0>0`$ and do not fix the sign of $`\vartheta _0`$. Developing a theory, we have a freedom with the sign of $`\vartheta _0`$; its value, however, can be roughly estimated as $`|\vartheta _0|\mathrm{}^2/2m`$. (Recall that $`\sigma _0E_g`$.)

If the exciton Bohr radius $`a_\mathrm{x}`$ is no more than several times larger than the lattice constant $`a_l`$, the second term in (2) can be important too. (In fact, this is an intermediate case between Frenkel exciton and Wannier one with $`a_\mathrm{x}a_l`$.) Then, the following coupling terms appear in Heisenberg equations

$$\stackrel{~}{E}_g\widehat{\psi }\stackrel{~}{E}_g\widehat{\psi }+\sigma _0_j\widehat{u}_j\widehat{\psi },$$
(4)
$$(\mathrm{}^2/2m)\mathrm{\Delta }\widehat{\psi }(\mathrm{}^2/2m)\mathrm{\Delta }\widehat{\psi }\vartheta _0_j\widehat{u}_j\mathrm{\Delta }\widehat{\psi }.$$
(5)
In terms of the lattice analog of the medium Hamiltonian, we made the hopping matrix element $`t`$ ($`t_{ij}`$ in $`\widehat{H}_{\mathrm{gas}}t_{ij}\widehat{\psi }_i^{}\widehat{\psi }_j`$) of the ‘lattice’ boson $`\widehat{\psi }_j`$ to be dependent on the ‘lattice’ deformation field $`\widehat{u}_j`$,

$$t_{ij}t_{ij}(u_i,u_j)t_{ij}+\stackrel{~}{\vartheta }_0(u_iu_j).$$

Meanwhile, the energy on a cite, $`\epsilon _{0,i}=\epsilon _0`$ (in $`\widehat{H}_{\mathrm{gas}}\epsilon _{0,i}\widehat{\psi }_i^{}\widehat{\psi }_i`$) depends on the lattice displacements too, f.ex.,

$$\epsilon _{0,i}\epsilon _0+\stackrel{~}{\sigma }_0(u_{i+1}u_{i1}),$$

see  for discussion.

## 3 Dynamic equations

Our aim is to investigate a special class of solutions of Hamiltonian (1), namely, we will search for the localized moving excitations (or packets). In the case of $`\vartheta _0=0`$, such a solution exists and the packet can move through the crystal with the constant velocity $`v`$ saving its shape in a manner as the solitons and kinks do in nonlinear media ,. To simplify the equations of motion, we choose the quasi-1D approximation. This means $`\widehat{\psi }(𝐱,t)\widehat{\psi }(x,t)`$, $`\widehat{𝐮}(𝐱,t)\widehat{u}_x(x,t)`$. In fact, we assume that the packet is inhomogeneous along the $`Ox`$ axis only, and $`vOx`$. Thus, it occupies all the area of $`S_{}`$. (For a discussion on the validity of 1D approximation, see .) Then, we can write the Heisenberg equations of motion as

$$i\mathrm{}_t\widehat{\psi }(x,t)=$$

$$=(\stackrel{~}{E}_g\frac{\mathrm{}^2}{2m}_x^2\vartheta _0_x\widehat{u}_x(x,t)_x^2+\nu _0\widehat{\psi }^{}\widehat{\psi }(x,t)+\nu _1\widehat{\psi }^{\mathrm{\hspace{0.17em}2}}\widehat{\psi }^2(x,t))\widehat{\psi }(x,t)+\sigma _0_j\widehat{u}_j(x,t)\widehat{\psi }(x,t),$$
(6)
$$\left(_t^2c_l^2_x^2\right)\widehat{u}_x(x,t)$$

$$c_l^22\kappa _3_x^2\widehat{u}_x_x\widehat{u}_x(x,t)=\rho ^1\sigma _0_x\left(\widehat{\psi }^{}\widehat{\psi }(x,t)\right)+\rho ^1\vartheta _0_x\left(_x\widehat{\psi }^{}_x\widehat{\psi }\right).$$
(7)
After some energy (and momentum) was pumped into the medium, an excited state of it appears near the boundary where the external energy (and momentum) was absorbed, see Fig. 1. Due to initial conditions, the excited state is a localized one, and it can be modeled as a droplet consisting of excitons (Bose-particles) and phonons. Moreover, it can acquire an average momentum directed along the $`Ox`$ axis because of unidirectional phonon production during the thermalization stage. Thus, the exciton-phonon droplet starts to move . The important assumption, however, is the appearance of a moving coherent field, or a condensate, from the localized excited state of a medium. For example, if the medium is a 3D crystal or an array of channels embedded into a matrix, it has to be coupled with a thermostat at low temperature. Moreover, the medium has to posses an inner mechanism of fast thermalization of the electronic excited states, for example, by emission of phonons .

It is known  that nonlinear elastic lattices support several types of localized excitations; some of them turn out to have the so-called nonzero direct current, $`u_xu_\mathrm{o}(x\overline{v}t).`$ For example, a moving packet can consist of two parts, such as

$$\widehat{u}_xu_\mathrm{o}(x\overline{v}t)+\underset{k}{}C(x\overline{v}t)\mathrm{e}^{i(kx\omega _kt)}\widehat{b}_k+\mathrm{h}.\mathrm{c}.,$$

and $`\overline{v}`$ is the group velocity of the packet. On the other hand, the Bose-gas has a property of multiple occupancy, so that the coherent mode can be introduced through the substitute $`\widehat{\psi }(x,t)\mathrm{\Psi }_\mathrm{o}(x,t)`$ and $`𝑑𝐱|\mathrm{\Psi }_\mathrm{o}(x,t)|^21`$. Thus, in the limit of $`T0`$ and within the long wavelength approximation, one can explore how Eqs. (6), (7) support the localized coherent solutions. In this article, we use the language and technique of Bose-Einstein condensation  that are suitable to take into account both the effect of many-particle coherency and the effect of many-particle nonlinearity. Note that the two parts of the coherent field, $`\mathrm{\Psi }_\mathrm{o}(x,t)`$ and $`u_\mathrm{o}(x,t)`$, are taken into account selfconsistently and we do not integrate out the phonons.

Substituting the coherent fields ($`c`$-functions) into operator equations (6) and (7), we obtain the following system of dynamic equations ($`\vartheta _0=(\mathrm{}^2/2m)\overline{\vartheta }_0`$ and $`\overline{\vartheta }_0`$ is dimensionless):

$$i\mathrm{}_t\mathrm{\Psi }_\mathrm{o}(x,t)=(\mathrm{}^2/2m)\overline{\vartheta }_0_xu_\mathrm{o}(x,t)_x^2\mathrm{\Psi }_\mathrm{o}(x,t)+$$

$$+(\stackrel{~}{E}_g\frac{\mathrm{}^2}{2m}_x^2+\nu _0|\mathrm{\Psi }_\mathrm{o}|^2(x,t)+\nu _1|\mathrm{\Psi }_\mathrm{o}|^4(x,t))\mathrm{\Psi }_\mathrm{o}(x,t)+\sigma _0_xu_\mathrm{o}(x,t)\mathrm{\Psi }_\mathrm{o}(x,t),$$
(8)
$$_t^2u_\mathrm{o}(x,t)c_l^2_x^2u_\mathrm{o}(x,t)c_l^22\kappa _3_x^2u_\mathrm{o}_xu_\mathrm{o}(x,t)=$$

$$=\rho ^1\sigma _0_x\left(\mathrm{\Psi }_\mathrm{o}^{}\mathrm{\Psi }_\mathrm{o}(x,t)\right)+\rho ^1\vartheta _0_x\left(_x\mathrm{\Psi }_\mathrm{o}^{}_x\mathrm{\Psi }_\mathrm{o}(x,t)\right).$$
(9)
This system can be considered as a generalization of the Zakharov system which appeared in Plasma Physics . The main difference is that, in the case of collective phenomena in Condensed Matter Physics, Eqs. (8) and (9) define the main part of the moving packet. (It is called the coherent Bose-core in this article.) In fact, in the case of $`T0`$ and, generally, in all the nonstationary cases, these two equations have to be coupled with another system of equations on the so-called out-of-condensate excitations . The dynamics of these excitation states can strongly influence the dynamics of the “parent” condensate.

In this article, we rely on a kind of adiabatic hypothesis to simplify the solution of the problem. If the occupancy of the coherent mode, $`N_\mathrm{o}(t)=𝑑𝐱|\mathrm{\Psi }_\mathrm{o}(x,t)|^2`$, and the occupancy of the out-of-condensate cloud $`\delta N(t)=𝑑𝐱\delta \psi ^{}\delta \psi (x,t)`$, change in time slowly and $`N_\mathrm{o}(t)\delta N(t)`$, one can try to find some quasi-stationary solution for the coherent core $`\mathrm{\Psi }_\mathrm{o}(x,t)`$ and balance its energy and momentum by a kind of “leakage” from the Bose-core into the incoherent out-of-condensate cloud and tail, $`N_\mathrm{o}N_\mathrm{o}\delta N(t)`$. This approach can be called the propagation of an exciton-phonon “comet” in a periodic medium, and it is valid within a finite time interval which is estimated below.

Recall that the standard ballistic ansatz with $`v=\mathrm{const}`$,

$$u_\mathrm{o}(x,t)=u_\mathrm{o}\left(xtv\right),$$

$$\mathrm{\Psi }_\mathrm{o}(x,t)=\mathrm{exp}\left(i(\phi _\mathrm{c}+k_0x)\right)\mathrm{exp}(i\mathrm{\Omega }(k_0)t)\psi _\mathrm{o}\left(xtv\right),$$

where

$$\mathrm{}\mathrm{\Omega }(k_0)=\stackrel{~}{E}_g+(\mathrm{}^2k_0^2/2m)+\mu ,\mathrm{}k_0=mv,$$

does not work properly in the case of $`\vartheta _00`$. Moreover, with such a ballistic ansatz, we always have $`v=v_\mathrm{s}`$, where $`v_{\mathrm{s},x}_x\phi _\mathrm{c}(x,t)=\mathrm{const}`$ is the superfluid velocity of the Bose-core. As the anomalous ballistic transport of exciton-phonon packets in semiconductor and quasi-1D structures resembles, to some extent, the well-known effect of superfluidity , we apply the methods of Many-Particle Physics to Eqs. (8) and (9). We assume that the following substitutes for the ballistic velocity and the chemical potential,

$$vv(x,t)v\left(1\pm \zeta (v)\psi _\mathrm{o}^2(x,t)+\mathrm{}\right),$$
(10)
$$\mu \mu (x,t)\mu \left(1\pm \stackrel{~}{\zeta }(v)\psi _\mathrm{o}^2(x,t)+\mathrm{}\right),$$
(11)
are more appropriate to model the nonlinear dynamics of the ballistic boson-phonon packets in the case of $`\vartheta _00`$. However, within the model with Hamiltonian (3), one can start, for example, from

$$vv(x,t)v_1\left(_xu_\mathrm{o}(x,t)\right)$$

and derive expansions (10) and (11).

To simplify the dynamic equations, we choose the following ansatz:

$$u_\mathrm{o}(x,t)=u_\mathrm{o}\left(xtv_1\left(_xu_\mathrm{o}(x,t)\right)\right),$$
(12)
$$\mathrm{\Psi }_\mathrm{o}(x,t)=\mathrm{exp}\left(i\phi _\mathrm{c}+ik_0_2\left(_xu_\mathrm{o}\right)xi\mathrm{\Omega }(k_0,_xu_\mathrm{o})t\right)\psi _\mathrm{o}\left(xtv_1\left(_xu_\mathrm{o}(x,t)\right)\right),$$
(13)
where

$$\mathrm{}\mathrm{\Omega }(k_0,_xu_\mathrm{o})\stackrel{~}{E}_g+(\mathrm{}^2k_0^2/2m)_3\left(_xu_\mathrm{o}(x,t)\right)+\mu _4\left(_xu_\mathrm{o}(x,t)\right).$$

For $`_j`$, we use the following expansion:

$$_j1+c_j\overline{\vartheta }_0_xu_\mathrm{o}(x,t)+\overline{c}_j\left(\overline{\vartheta }_0_xu_\mathrm{o}(x,t)\right)^2,j=1,\mathrm{},4,$$

and $`\mathrm{}k_0=mv`$. Note that, in some sense, ansatz (13) is the well-known representation

$$\mathrm{\Psi }_\mathrm{o}(x,t)=\sqrt{n_\mathrm{o}(x,t)}\mathrm{exp}\left(i\phi _\mathrm{c}(x,t)\right),$$

with the inhomogeneous fields $`n_\mathrm{o}(x,t)`$, $`v_{\mathrm{s},x}(x,t)_x\phi _\mathrm{c}(x,t)`$, and $`\mu (x,t)_t\phi _\mathrm{c}(x,t)`$. In this article, we will not write out and solve a set of nonlinear hydrodynamic equations  on these variables. Instead, we are interested in a microscopic approach aimed to clarify the role of the coupling parameters, $`\vartheta _0`$ and $`\sigma _0`$. Obviously, for a packet moving nonuniformly, one might expect some deviation from the simple low

$$\phi _\mathrm{c}(𝐱,t)=\phi _\mathrm{c}+𝐩_s𝐱/\mathrm{}(\stackrel{~}{E}_g+𝐩_\mathrm{s}^2/2m+\mu )t/\mathrm{}.$$

If a good initial approximation for such a deviation can be found, one can finish with a kind of relatively simple NLS equation on the envelope functions of the Bose-core, $`\mathrm{exp}(i\phi _\mathrm{c})\psi _\mathrm{o}(x)`$ and $`_xu_\mathrm{o}(x)`$. Such a mathematical structure do appears from Eqs. (8) and (9) within the validity of the adiabatic hypothesis if one starts from quasistationary ansatz (12) and (13).

Recall that the following relation between $`\psi _\mathrm{o}(x)`$ and $`_xu_\mathrm{o}(x)`$ is a good approximation for the selfconsistent boson-phonon coherent state at $`\vartheta _0=0`$,

$$_xu_\mathrm{o}(x,t)=f(\psi _\mathrm{o}^2(x,t))A_2\psi _\mathrm{o}^2(x,t)+A_4\psi _\mathrm{o}^4(x,t).$$

In fact, this approximation is valid in a more general case (although it is a part of a more complicated relation, such as $`_xu_\mathrm{o}=f(\psi _\mathrm{o},_x\psi _\mathrm{o})`$ ). Then, for the coherent phase $`\phi _\mathrm{c}(x,t)`$, we can write out the following symmetric representation:

$$\phi _\mathrm{c}(x,t)\phi _\mathrm{c}+k_0\left\{1+\zeta _2(v)\psi _\mathrm{o}^2(x,t)+\mathrm{}\right\}x$$

$$\left(\stackrel{~}{E}_g+\mathrm{}k_0^2/2m\{1+\zeta _3(v)\psi _\mathrm{o}^2(x,t)+\mathrm{}\}\mu \{1+\zeta _4(v)\psi _\mathrm{o}^2(x,t)+\mathrm{}\}\right)t/\mathrm{}.$$
(14)
The ballistic velocity of the condensate, $`xxv(x,t)t`$, is modified as follows:

$$xxv\left\{1+\zeta _1(v)\psi _\mathrm{o}^2(x,t)+\mathrm{}\right\}t.$$
(15)
Here, we introduce the parameters $`\zeta _j(v)c_j\zeta (v)\vartheta _0`$ and expect $`\zeta _1(v)\zeta _2(v)\zeta _3(v)\zeta _4(v)`$. If $`c_21`$, we get immediately $`c_1=1+c_22`$ and $`c_3=1+2c_23`$ within the quasistationary approximation for the dynamic equations. The next field terms ( $`\overline{c}_j\vartheta _0^2\psi _\mathrm{o}^4(x,t)`$ ) in expansions (14), (15) can be defined through $`c_2`$ and $`\overline{c}_2`$ as well. However, they are small in a weakly nonlinear case and we will disregard them. The constant $`c_4`$ will be defined to simplify the quasistationary equation on $`\psi _\mathrm{o}^2(x)`$.

To clarify the dynamic properties of the generalized ballistic ansatz and define the unknown $`c_j`$, we will calculate the energy and momentum of a moving soliton-like solution with the factors $`_j(x,t)`$. Meanwhile, the envelope functions, $`\mathrm{exp}(i\phi _\mathrm{c})\psi _\mathrm{o}(x)`$ and $`_xu_\mathrm{o}(x)`$, are found from the quasistationary equations one has to derive. These equations read

$$\mu \left\{1+c_4\overline{\vartheta }_0_xu_\mathrm{o}(x,t)+\mathrm{}\right\}\psi _\mathrm{o}(x,t)(\mathrm{}^2/2m)\overline{\vartheta }_0_xu_\mathrm{o}(x,t)_x^2\psi _\mathrm{o}(x,t)+$$

$$+(\frac{\mathrm{}^2}{2m}_x^2+\nu _0\psi _\mathrm{o}^2(x,t)+\nu _1\psi _\mathrm{o}^4(x,t))\psi _\mathrm{o}(x,t)+\sigma _0_xu_\mathrm{o}(x,t)\psi _\mathrm{o}(x,t),$$
(16)
$$v^2_x^2u_\mathrm{o}(x,t)\left(1+c_1\overline{\vartheta }_0_xu_\mathrm{o}+\mathrm{}\right)^2c_l^2_x^2u_\mathrm{o}(x,t)$$

$$c_l^2\left(2\kappa _3_x^2u_\mathrm{o}_xu_\mathrm{o}(x,t)+3\kappa _4_x^2u_x\left(_xu_\mathrm{o}(x,t)\right)^2\right)$$

$$\rho ^1\sigma _0_x\left(\psi _\mathrm{o}^2(x,t)\right)+\rho ^1\vartheta _0_x\left(_x\mathrm{\Psi }_\mathrm{o}^{}_x\mathrm{\Psi }_\mathrm{o}(x,t)\right).$$
(17)
Note that we did not include the terms depending on time explicitly into these equations. The consequent divergencies, however, can be cured by taking into account the incoherent parts of the packet. In this article, we propose the effect of slow dynamic redistribution of the occupation numbers between the condensate and out-of-condensate cloud, while both of these states can be easily calculated in the quasistationary approximation. Recall that dynamic mean field interaction between the condensate and noncondensed part of the packet leads to the damping of both of them . On the other hand, the time interval, during which the exciton-phonon packet with the coherent Bose-core moves in the periodic medium, is always finite, especially for the “clean” samples without lattice imperfections. Therefore, the adiabatic hypothesis can be used in our case.

To find a soliton-like solution of Eqs. (16) and (17), we reduce Eq. (17) to the following form:

$$(c_l^2v^2)_x^2u_\mathrm{o}(x,t)c_l^22\overline{\kappa }_3_x^2u_\mathrm{o}_xu_\mathrm{o}(x,t)=$$

$$=\rho ^1\overline{\sigma }_0_x\left(\psi _\mathrm{o}^2(x,t)\right).$$
(18)
Note that the interaction vertices $`\kappa _3`$ and $`\sigma _0`$ are slightly renormalized in this equation, $`\overline{\kappa }_3=\overline{\kappa }_3(\vartheta _0)`$, $`\overline{\sigma }_0=\overline{\sigma }_0(\vartheta _0)`$. For example, we can write $`\overline{\kappa }_3`$ as follows ($`c_21`$, $`c_12`$)

$$\overline{\kappa }_3=|\kappa _3|+2(v^2/c_l^2)|\overline{\vartheta }_0|<\mathrm{\hspace{0.17em}0},$$
(19)
and $`|\overline{\kappa }_3|<|\kappa _3|`$ for $`\overline{\vartheta }_0<0`$. Here, we choose $`\kappa _3<0`$ ($`\kappa _4>0`$) and assume the inequality

$$|\kappa _3|>|\overline{\vartheta }_0|$$

to be valid. Then, the sign and the order of value of $`\kappa _3`$ remains the same after the renormalization. For example, we take $`|\kappa _3|510`$, whereas $`|\overline{\vartheta }_0|12`$, and, for $`v<c_l`$, we can estimate $`v^2/c_l^20.50.9`$.

One can represent the solution of Eq. (18) in the following form:

$$_xu_\mathrm{o}(x,t)A_2\psi _\mathrm{o}^2(x,t)+A_4\psi _\mathrm{o}^4(x,t),$$

where

$$A_2\overline{\sigma }_0/\rho (c_l^2v^2)=\gamma (v)\left(\overline{\sigma }_0/Mc_l^2\right)a_l^3>0,$$

$$A_4\gamma (v)|\overline{\kappa }_3|A_2^2>0,$$

and $`\rho M/a_l^3`$, $`\gamma (v)=c_l^2/(c_l^2v^2)>1`$. Immediately, we can rewrite the factor $``$ for the ballistic velocity ($`\overline{\vartheta }_0<0`$) as follows:

$$_1(_xu_\mathrm{o}(x,t))_1\left(|\mathrm{\Psi }_\mathrm{o}(x,t)|^2\right)$$

$$\mathrm{\hspace{0.17em}1}+|\zeta _1(v)|\psi _\mathrm{o}^2(x,t)\overline{\zeta }_1(v)\psi _\mathrm{o}^4(x,t).$$
(20)
For $`\overline{\vartheta }_0>0`$, we estimate

$$_1(_xu_\mathrm{o}(x,t))_1\left(|\mathrm{\Psi }_\mathrm{o}(x,t)|^2\right)$$

$$\mathrm{\hspace{0.17em}1}|\zeta _1(v)|\psi _\mathrm{o}^2(x,t)+\overline{\zeta }_1(v)\psi _\mathrm{o}^4(x,t).$$
(21)
We use the following representation for $`|\mathrm{\Psi }_\mathrm{o}(x,t)|^2`$:

$$|\mathrm{\Psi }_\mathrm{o}(x,t)|^2=\psi _\mathrm{o}^2(xv(x,t)t)\mathrm{\Phi }_\mathrm{o}^2f^2(\overline{x}/L_0),$$

where $`f(\overline{x}/L_0)`$ is the dimensionless ‘shape’ function, $`|f(x/L_0)|1`$, $`\mathrm{\Phi }_\mathrm{o}`$ is the amplitude of the coherent state, and $`xv(x,t)t\overline{x}`$. The characteristic width of the condensate is estimated as $`(36)L_0`$. Although the terms $`\psi _\mathrm{o}^4`$ in Eqs. (20) and (21) could be important in strong nonlinear regimes (as well as other omitted contributions, e.g., $`(_x\psi _\mathrm{o})^2`$ ), we assume $`\zeta _1(v)\mathrm{\Phi }_\mathrm{o}^2<1`$ and $`\zeta _1(v)\mathrm{\Phi }_\mathrm{o}^2\overline{\zeta }_1(v)\mathrm{\Phi }_\mathrm{o}^4`$ in this article. Then, we can roughly estimate the dynamic factors $`_j`$ as

$$_j\left(|\mathrm{\Psi }_\mathrm{o}(x,t)|^2\right)1\pm c_j\zeta (v)\mathrm{\Phi }_\mathrm{o}^2f^2(x,t),c_j1,$$
(22)
and the question is how strong the factor $``$ can deviate from 1. If $`|\overline{\vartheta }_0|1`$ and $`\gamma (v)510`$ ($`v<c_l`$), we can estimate the dimensionless factor $`c_j\zeta _j(v)/a_l^3`$ that enters Eq. (22) as follows

$$c_j\zeta _j(v)/a_l^3c_j|\overline{\vartheta }_0|\gamma (v)\left(\overline{\sigma }_0/Mc_l^2\right)1.$$

Then, the most important parameter in Eq. (22) remains to be $`a_l^3\mathrm{\Phi }_\mathrm{o}^2<\mathrm{\hspace{0.17em}1}`$.

Now, we can qualitatively describe the dynamics of bosonic core of the packet, see Fig. 2. The top of the soliton $`|\mathrm{\Psi }_\mathrm{o}(xv(x,t)t)|^2\mathrm{\Phi }_\mathrm{o}^2f^2(\overline{x}/L_0)`$ moves with

$$vv\pm (\zeta _1(v)\mathrm{\Phi }_\mathrm{o}^2)v,f_{\mathrm{top}}(x/L_0)1,$$

whereas the slopes of $`|\mathrm{\Psi }_\mathrm{o}(x,t)|^2`$ move with

$$vv\pm 10^1(\zeta _1(v)\mathrm{\Phi }_\mathrm{o}^2)v,f_{\mathrm{slope}}(x/L_0)10^1.$$

The ballistic velocity of the tails of $`|\mathrm{\Psi }_\mathrm{o}(x,t)|^2`$ remains to be almost unchanged, i.e. $`=v`$. One can introduce the important parameter $`\delta v_{\mathrm{top}}`$ as

$$\delta v_{\mathrm{top}}\pm c_1|\zeta (v)|\mathrm{\Phi }_\mathrm{o}^2v$$

and the dimensionless ratio

$$\frac{\delta v_{\mathrm{top}}t}{L_\mathrm{o}}\pm c_1|\zeta (v)|\mathrm{\Phi }_\mathrm{o}^2vt/L_\mathrm{o}.$$
(23)
Then, the time scale $`\mathrm{\Delta }t`$, during which our quasistationary ansatz can be used to describe the dynamics of the core of the total solitonic packets, can be roughly estimated from $`|\delta v_{\mathrm{top}}|\mathrm{\Delta }t/L_\mathrm{o}\mathrm{\hspace{0.17em}1}`$. In fact, parameter (23) depends on time because the dynamic leakage from the Bose-core into the tail and coma of the “comet” leads to $`\mathrm{\Phi }_\mathrm{o}(t)`$ and $`L_\mathrm{o}(t)`$. To estimate the values of $`|\mu |`$, $`a_l^3\mathrm{\Phi }_\mathrm{o}^2`$, and $`L_0`$, we have to solve the stationary equation on $`\psi _\mathrm{o}(x)`$.

## 4 Solution of stationary equation

As a result of all the simplifications we made, we obtain this equation in the following form (we choose $`c_4=1`$):

$$\mu \left\{1+\zeta (v)\psi _\mathrm{o}^2(x)\overline{\zeta }(v)\psi _\mathrm{o}^4(x)\right\}\psi _\mathrm{o}(x)=$$

$$=(\frac{\mathrm{}^2}{2m}\{\mathrm{\hspace{0.17em}1}+\zeta (v)\psi _\mathrm{o}^2(x)\overline{\zeta }(v)\psi _\mathrm{o}^4(x)\}_x^2|\stackrel{~}{\nu _0}(v)|\psi _\mathrm{o}^2(x)+\stackrel{~}{\nu _1}(v)\psi _\mathrm{o}^4(x))\psi _\mathrm{o}(x).$$
(24)
This equation resembles combined Ablowitz-Ladik and Nonlinear Schrödinger Equations taken in the continuous limit ,. Note that the boson interaction vertices are strongly renormalized because of interaction with the coherent phonon field $`_xu_\mathrm{o}`$. We assume the following conditions to be valid

$$\nu _0>0\stackrel{~}{\nu _0}(v/c_l,\overline{\sigma }_0)<0,\nu _1>0\stackrel{~}{\nu _1}(v/c_l,\overline{\sigma }_0,\overline{\kappa }_3)>0.$$
(25)
If the dimensionless parameter $`\zeta (v)\mathrm{\Phi }_\mathrm{o}^2a_l^3\mathrm{\Phi }_\mathrm{o}^2`$ is small enough, one can reduce Eq. (24) to Nonlinear Schrödinger Equation. We use the representation

$$\psi _\mathrm{o}(x)=\mathrm{\Phi }_\mathrm{o}f(\beta (\mathrm{\Phi }_\mathrm{o})x,\eta (\mathrm{\Phi }_\mathrm{o})),\mu =|\mu |<0,$$

with unknown functions $`\beta (\mathrm{\Phi }_\mathrm{o})L_0^1`$ and $`\eta (\mathrm{\Phi }_\mathrm{o})`$. Then, the following simple equation appears as an approximation of Eq. (24):

$$|\mu |f(x)|\overline{\nu }_0(v)|\mathrm{\Phi }_\mathrm{o}^2f^3(x)+\overline{\nu }_1(v)\mathrm{\Phi }_\mathrm{o}^4f^5(x)\frac{\mathrm{}^2}{2m}\beta ^2(\mathrm{\Phi }_\mathrm{o})_x^2f(x).$$
(26)
This equation is the so-called “subcritical” NLS equation. In Condensed Matter Physics, it is used in the theory of superfluidity and Bose-Einstein condensation . It is also applied in nonlinear optics for light pulses in the medium with a cubic-quintic nonlinearity , and it is known as the Lienard equation in the theory of exactly solvable nonlinear equations .

Here, we apply it to describe the coherent state of excitons and phonons with macroscopic occupancy. The following approximation is used for the renormalized vortices: $`\overline{\nu }_0\stackrel{~}{\nu _0}`$ and $`\overline{\nu }_1(v)\stackrel{~}{\nu _1}(v)+|\stackrel{~}{\nu _0}(v)|\zeta (v)`$. To estimate their strength, we use the following formulas:

$$|\stackrel{~}{\nu _0}(v)|=\gamma (v)\left(\overline{\sigma }_0/Mc_l^2\right)(\sigma _0a_l^3)\nu _0>0,$$
(27)
which as valid at $`\gamma (v)>\gamma _\mathrm{o}`$, or, equivalently, $`v_\mathrm{o}<v<c_l`$ (for discussion, see ,),

$$\stackrel{~}{\nu _1}(v)=(\gamma (v)|\overline{\kappa }_3|)\left\{\gamma (v)\left(\frac{\overline{\sigma }_0}{Mc_l^2}\right)\right\}^2(\sigma _0a_l^6)+\nu _1,$$
(28)
and

$$\overline{\nu }_1\mathrm{\Phi }_\mathrm{o}^4\stackrel{~}{\nu _1}\mathrm{\Phi }_\mathrm{o}^4+|\stackrel{~}{\nu _0}|\mathrm{\Phi }_\mathrm{o}^2\left(\zeta (v)\mathrm{\Phi }_\mathrm{o}^2\right).$$
(29)
Note that in the case of $`\vartheta _0<0`$ we obtain the enhanced three particle vertex, and the nonlinear corrections are important to estimate its value, $`\nu _1\stackrel{~}{\nu _1}\overline{\nu }_1>0`$.

The localized solution of Eq. (26) exists if the following inequalities are valid :

$$|\mu |<\mu ^{}=\frac{3}{16}\frac{|\stackrel{~}{\nu _0}|^2}{\overline{\nu }_1},\mathrm{\Phi }_\mathrm{o}^2<\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}=\frac{3}{4}\frac{|\stackrel{~}{\nu _0}|}{\overline{\nu }_1}.$$
(30)
Then, $`|\mu |=|\mu |(\mathrm{\Phi }_\mathrm{o})`$ and we have

$$|\mu |/\mu ^{}=2\mathrm{\Phi }_\mathrm{o}^2/\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}\left(\mathrm{\Phi }_\mathrm{o}^2/\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}\right)^2<1.$$

As a rough estimate, we can use the inverse formula $`\mathrm{\Phi }_\mathrm{o}^2(\mu )0.5(|\mu |/\mu ^{})\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}`$ valid at $`|\mu |/\mu ^{}1`$.

The representation $`|\mu |=(\mathrm{}^2/2m)L_0^2`$ leads to an easy estimate of the characteristic width of the soliton. Indeed, we can introduce the length $`L_{}`$ through $`\mu ^{}=(\mathrm{}^2/2m)L_{}^2`$ and obtain the following representation

$$\beta (\mathrm{\Phi }_\mathrm{o})\beta (\mu )=\sqrt{\frac{2m}{\mathrm{}^2}\mu ^{}(|\mu |/\mu ^{})},\beta (|\mu |)=L_{}^1\sqrt{|\mu |/\mu ^{}}.$$
(31)
Then, $`L_0=L_{}\sqrt{\mu ^{}/|\mu |}`$ and, always, $`L_0>L_{}`$. To estimate the important microscopic parameters $`\mu ^{}`$ and $`L_{}`$, we write the renormalized vertices in the following form:

$$|\stackrel{~}{\nu _0}(v)|=\stackrel{~}{\alpha }_0(v)\mathrm{Ry}_\mathrm{x}a_\mathrm{x}^3\mathrm{and}\overline{\nu }_1(v)=\overline{\alpha }_1(v)\mathrm{Ry}_\mathrm{x}a_\mathrm{x}^6,$$
(32)
and estimate $`\stackrel{~}{\alpha }_0\overline{\alpha }_110^1`$. (In theory, $`\stackrel{~}{\alpha }_0`$ can vary from $`10^2`$ to $`1`$ with changement of $`v`$ within $`v_\mathrm{o}<v<c_l`$ ,.) Then, we have

$$\mu ^{}(v)=\frac{3}{16}\frac{\stackrel{~}{\alpha }_0^2}{\stackrel{~}{\alpha }_1}\mathrm{Ry}_\mathrm{x}10^2\mathrm{Ry}_\mathrm{x},L_{}^2(v)\frac{\stackrel{~}{\alpha }_1}{\stackrel{~}{\alpha }_0^2}a_\mathrm{x}^210a_\mathrm{x}^2.$$
(33)
As the following inequality is valid

$$\zeta _j(v)\mathrm{\Phi }_\mathrm{o}^2<\zeta _j(v)\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}},$$

we can define the upper limit of $`\delta v_{\mathrm{top}}/v\zeta _1(v)\mathrm{\Phi }_\mathrm{o}^2`$ in our model. We have

$$\zeta _1(v)\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}(v)|\overline{\vartheta }_0|c_1\left\{\gamma (v)\frac{\overline{\sigma }_0}{Mc_l^2}\right\}\left(a_l^3\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}\right)\mathrm{\hspace{0.17em}0.1}1,$$
(34)
where

$$a_l^3\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}(v)(\stackrel{~}{\alpha }_0(v)/\overline{\alpha }_1(v))(a_l^3/a_\mathrm{x}^3)0.1.$$
(35)
As a result, one can define the meaning of the weakly nonlinear case in our model. For example, if the parameter $`\mathrm{\Phi }_\mathrm{o}^2/\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}<0.10.5`$, the assumption

$$\overline{\zeta }(v)\mathrm{\Phi }_\mathrm{o}^4\zeta (v)\mathrm{\Phi }_\mathrm{o}^2$$

we made to simplify the factors $`_j(_xu_\mathrm{o}(x,t))`$ to $`_j(\psi _\mathrm{o}^2(x,t))`$ is correct.

We write the exact solution of Eq. (26) in the following form:

$$\psi _\mathrm{o}^2(x)=\frac{\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}\left(|\mu |/2\mu ^{}\right)}{\sqrt{1|\mu |/\mu ^{}}\mathrm{cosh}^2\left(\beta (|\mu |)x\right)+(1/2)(1\sqrt{1|\mu |/\mu ^{}})}$$
(36)
Note that

$$\psi _\mathrm{o}^2(x=0)=\mathrm{\Phi }_\mathrm{o}^2(|\mu |)=\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}(1\sqrt{1|\mu |/\mu ^{}}).$$

If $`|\mu |\mu ^{}`$, we can use the following asymptotics ( $`\eta (|\mu |/\mu ^{})0`$):

$$\psi _\mathrm{o}(x)2\mathrm{\Phi }_\mathrm{o}\mathrm{exp}(\beta (\mathrm{\Phi }_\mathrm{o})|x|)\mathrm{at}|x|>2\beta (\mathrm{\Phi }_\mathrm{o})^1.$$
(37)
In addition, the coherent phonon part $`u_\mathrm{o}(x,t)`$ can be rewritten as follows

$$_xu_\mathrm{o}(x)\gamma (v)\left(\overline{\sigma }_0/Mc_l^2\right)(a_l^3\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}})\left(\mathrm{\Phi }_\mathrm{o}^2/\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}\right)f^2\left(\beta (|\mu |)x\right)+$$

$$+\gamma (v)|\overline{\kappa }_3|\left\{\gamma (v)\left(\overline{\sigma }_0/Mc_l^2\right)(a_l^3\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}})\right\}^2\left(\mathrm{\Phi }_\mathrm{o}^2/\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}\right)^2f^4\left(\beta (|\mu |)x\right).$$
(38)
We apply the 3D normalization condition on the macroscopic wave function $`\mathrm{\Psi }_0(x,t)`$ as follows

$$|\mathrm{\Psi }_0|^2(x,t)𝑑𝐱=N_\mathrm{o}(t)1.$$

Substituting the quasistationary ansatz into the integrand, we have

$$S\psi _\mathrm{o}^2\left(xv(1+\zeta _1(v)\psi _\mathrm{o}^2(x,t)+\mathrm{})t\right)𝑑x=N_\mathrm{o}(t)$$
(39)
$$S\psi _\mathrm{o}^2(x^{})𝑑x^{}=N_\mathrm{o}.$$
(40)
Here, $`N_\mathrm{o}1`$ has a meaning of the number of particles inside the quasistationary condensate, (it is a macroscopic number), and $`SL_{}^2S_{}`$, where $`S_{}`$ is the cross section area of the crystal. Note that, due to the symmetry of $`f(x/L_0)`$ against $`xx`$, the value of $`N_\mathrm{o}`$ is conserved, but as a first approximation.

To understand how the macroscopic wave function is formed from the microscopic parameters of the theory (they are $`\mu ^{}`$, $`\mathrm{\Phi }_\mathrm{o}^{}`$, $`L_{}`$, and the renormalized vertices $`\stackrel{~}{\nu }(v)`$), it is useful to introduce the parameter $`N^{}=SL_{}\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}`$ as well. We estimate $`N^{}`$ as follows:

$$N^{}(v)\left(1/\sqrt{2\overline{\alpha }_1(v)}\right)(S/a_\mathrm{x}^2)1.$$
(41)
Then, $`N^{}(v)(S/a_\mathrm{x}^2)`$ is always the macroscopic number in 3D case ($`\overline{\alpha }_1(v)10^110^2`$). In fact, Eq. (40) leads to an algebraic equation that relates the macroscopic parameters with the microscopic ones. As a result, $`\sqrt{|\mu |/\mu ^{}}`$ can be defined as a function of $`\mathrm{exp}(2N_\mathrm{o}/N^{})`$. To avoid cumbersome formulas, we use simple estimates:

$$\sqrt{|\mu |/\mu ^{}}=f(N_\mathrm{o}/N^{})N_\mathrm{o}/N^{}\mathrm{and}L_0/L_{}(N_\mathrm{o}/N^{})^1.$$
(42)
They are valid up to $`N_\mathrm{o}/N^{}0.10.5`$.

Thus, we can estimate the important dimensionless parameter that characterizes the deviation of the ballistic velocity $`v(x,t)`$ from $`v=\mathrm{const}`$ and the chemical potential $`\mu (x,t)`$ from $`|\mu |=\mathrm{const}`$. We can write

$$\mathrm{\Phi }_\mathrm{o}^2/\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}=\mathrm{f}(|\mu |/\mu ^{})(1/2)(|\mu |/\mu ^{})(1/2)(N_\mathrm{o}/N^{})^2.$$
(43)
This means

$$\mathrm{\Phi }_\mathrm{o}^2/\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}10^310^2\mathrm{while}N_\mathrm{o}/N^{}(v)10^2\mathrm{\hspace{0.17em}10}^1,$$

and

$$\mathrm{\Phi }_\mathrm{o}^2/\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}0.2\mathrm{\hspace{0.17em}0.4}\mathrm{while}N_\mathrm{o}/N^{}(v)0.60.7\mathrm{\hspace{0.17em}1}.$$

As a result, the reasonable values for the macroscopic parameter $`\zeta (v)\mathrm{\Phi }_\mathrm{o}^2\delta v_{\mathrm{top}}/v`$ can be estimated as

$$\left\{\zeta (v)\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}\right\}\left(\mathrm{\Phi }_\mathrm{o}^2/\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}\right)10^310^2,$$
(44)
(the microscopic $`\zeta (v)\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}0.11`$). Now, we can rewrite parameter (23) to emphasize the interplay between the microscopic and macroscopic parameters of the theory:

$$\frac{\delta v_{\mathrm{top}}t}{L_\mathrm{o}}\pm \mathrm{\hspace{0.17em}0.5}c_1\left\{|\zeta (v)|\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}}\right\}\frac{vt}{L_{}}(N_\mathrm{o}/N_{})^3,\zeta (v)\vartheta _0\sigma _0.$$
(45)
On Fig. 2, we show how the coherent part of the packet moves in the medium if one switch on the interaction $`\vartheta _00`$. The initial state is taken to be symmetric (the exact solution for $`\sigma _00`$, $`\vartheta _0=0`$ is used), and its evolution is obtained by use of quasistationary ansatz (12), (13) which conserves the amplitude and the characteristic width of the moving soliton. Although it is a rough estimate of the real dynamics of the Bose-core, it gives some understanding of how the microscopic exciton-phonon interaction controls the dynamics. In addition, it can be easily adjusted to a more realistic case of the moving core and out-of-condensate cloud.

## 5 Energy of the moving soliton

Recall that the standard ballistic boson-phonon soliton, $`\mathrm{\Psi }_\mathrm{o}(x,t)=\mathrm{exp}\left(i\phi _\mathrm{c}(x,t)\right)\varphi _\mathrm{o}(xvt)`$ and $`_xu_\mathrm{o}(xvt)`$ with $`\phi _\mathrm{c}(x,t)=\phi _\mathrm{o}+k_0x\omega _0t`$, is a stationary state ($`\vartheta _0=0`$), that is

$$E_\mathrm{o}=𝑑𝐱𝒯_0^0(x,t)=\mathrm{const}\mathrm{and}𝐏_\mathrm{o}P_{\mathrm{o},x}=𝑑𝐱𝒫_x(x,t)=\mathrm{const}.$$

At $`\vartheta _00`$, we obtained an analog of such a solution, starting from Eqs. (12), (13) with $`vv(x,t)`$, $`k_0k_0(x,t)`$, and $`\mu \mu (x,t)`$. In this article, the quasistationary approximation is the most important one we used to simplify the dynamic equations for the core of the packet. Therefore, the exact calculation of $`E_\mathrm{o}`$ leads to $`_tE_\mathrm{o}0`$ for $`\mathrm{\Psi }_\mathrm{o}(x,t)=\mathrm{exp}\left(i\phi _\mathrm{c}(x,t,\psi _\mathrm{o}(x,t))\right)\psi _\mathrm{o}(xv(x,t)t)`$ with $`N_\mathrm{o}=\mathrm{const}`$ and $`L_0=\mathrm{const}`$. Within the validity of the adiabatic hypothesis (slow exchange between $`N_\mathrm{o}(t)`$ and $`\delta N(t)`$ while $`N_\mathrm{o}(t)\delta N(t)`$ ), one has, first, to calculate $`_tE_\mathrm{o}`$ of the quasistationary Bose-core and, second, to cure the obtained divergency by $`_tE_{\mathrm{tot}}=_tE_\mathrm{o}+_t\delta E=0`$ ($`T0`$).

To calculate the energy of the moving condensate, we have to integrate the zero component of the energy-momentum tensor $`𝒯_0^0`$ over the spatial coordinates. We have the following formula (written in the laboratory frame):

$$𝒯_0^0(𝐱,t)=\stackrel{~}{E}_g\mathrm{\Psi }_\mathrm{o}^{}\mathrm{\Psi }_\mathrm{o}+\frac{\mathrm{}^2}{2m}\mathrm{\Psi }_\mathrm{o}^{}\mathrm{\Psi }_\mathrm{o}+\frac{\nu _0}{2}\left(\mathrm{\Psi }_\mathrm{o}^{}\right)^2\mathrm{\Psi }_\mathrm{o}^2+\frac{\nu _1}{3}\left(\mathrm{\Psi }_\mathrm{o}^{}\right)^3\mathrm{\Psi }_\mathrm{o}^3+$$

$$+\frac{\rho }{2}(_tu_\mathrm{o})^2+\frac{\rho c_l^2}{2}(_xu_\mathrm{o})^2+\frac{\rho c_l^2}{3}\kappa _3(_xu_\mathrm{o})^3+\sigma _0_xu_\mathrm{o}\mathrm{\Psi }_\mathrm{o}^{}\mathrm{\Psi }_\mathrm{o}+\overline{\vartheta }_0\frac{\mathrm{}^2}{2m}_xu_\mathrm{o}\mathrm{\Psi }_\mathrm{o}^{}\mathrm{\Psi }_\mathrm{o}.$$

After integration of $`𝒯_0^0`$, we conclude that there are no terms $`t^1`$ in $`E_\mathrm{o}(t)`$,

$$E=E_\mathrm{o}+(\delta E)(vt/L_0)^2+\mathrm{}e_0N_\mathrm{o}+\delta e_0N_\mathrm{o}\left(\delta v_{\mathrm{top}}t/L_0\right)^2.$$
(46)
In fact, the dependence on $`N_\mathrm{o}`$ is a nonlinear one,

$$Ee_0(N_\mathrm{o}/N_{})N_\mathrm{o}+\delta e_0(N_\mathrm{o}/N_{})N_\mathrm{o}\left\{\delta v_{\mathrm{top}}(N_\mathrm{o})t/L_0(N_\mathrm{o})\right\}^2.$$

We found $`\delta e_0>0`$ for both $`\vartheta _0>0`$ ($`\delta v_{\mathrm{top}}<0`$) and $`\vartheta _0<0`$ ($`\delta v_{\mathrm{top}}>0`$). Indeed, we can write out the expansion

$$\delta e_0\delta e_0^{(0)}+\delta e_0^{(1)}(\zeta (v)\mathrm{\Phi }_\mathrm{o}^2)+\mathrm{}\delta e_0^{(0)}>\mathrm{\hspace{0.17em}0},$$

and $`\delta e_0^{(1)}\delta e_0^{(0)}`$. Recall the structure of the stationary part, $`E_\mathrm{o}=e_0N_\mathrm{o}=E_\mathrm{x}+E_{\mathrm{int}}+E_{\mathrm{ph}}`$. We have

$$e_0(\mathrm{}^2k_0^2/2m)(\mathrm{\hspace{0.17em}1}+\mathrm{const}\zeta (v)\mathrm{\Phi }_\mathrm{o}^2+\mathrm{})\mathrm{const}^{}|\mu |+e_{\mathrm{ph}},$$
(47)
where $`\mathrm{const},\mathrm{const}^{}1`$ and $`e_{\mathrm{ph}}e_{\mathrm{ph}}^{(2)}e_{\mathrm{ph}}^{(3)}`$. Here, the harmonic contribution into $`E_{\mathrm{ph}}`$ is $`e_{\mathrm{ph}}^{(2)}N_\mathrm{o}>0`$, and

$$e_{\mathrm{ph}}^{(2)}\frac{M(v^2+c_l^2)}{2}\left(\gamma (v)\frac{\sigma _0}{Mc_l^2}\right)^2(a_l^3\mathrm{\Phi }_\mathrm{o}^{\mathrm{\hspace{0.17em}2}})\left\{0.5|\mu |/\mu ^{}\right\}.$$

We can compare formula (47) with our estimate of $`\delta e_0^{(0)}`$:

$$\delta e_0(\mathrm{}^2\stackrel{~}{k}_0^2/2m)(c_3^2/4c_1^2)\mathrm{const}_1+|\mu |\mathrm{\hspace{0.17em}3}\mathrm{const}_2+e_{\mathrm{ph}}^{(2)}\mathrm{const}_3>\mathrm{\hspace{0.17em}0},$$
(48)
where all the constants are small (as a rough estimate, $`\mathrm{const}_j10^1`$). Then, $`\delta e_0e_0`$.

The result (46) means that changement of the shape of the Bose-core, which was included into the generalized ballistic ansatz to satisfy the dynamic equations with $`\vartheta _00`$, costs energy for both $`\delta v_{\mathrm{top}}>0`$ and $`\delta v_{\mathrm{top}}<0`$ during the observation time. Note that the width of the soliton, $`(36)L_0(N_\mathrm{o})`$, $`L_0^2|\mu |(N_\mathrm{o})`$, and the number of Bose-particles forming the soliton, $`N_\mathrm{o}`$, are conserved as it was presented on Fig. 2, but this is only the first approximation within the adiabatic hypothesis. It is possible, however, to balance the energy of the total boson-phonon packet. (Here, the image of an exciton-phonon “comet” helps: the Bose-core with $`N_\mathrm{o}1`$ is only a part of the moving delocalized packet.)

First, we can calculate the value of $`\{_{N_\mathrm{o}}E_\mathrm{o}\}\delta N`$. We take it as

$$\left\{_{N_0}E_\mathrm{o}\right\}\delta N\stackrel{~}{e_0}\delta N,$$

where we have the following formula for $`\stackrel{~}{e_0}`$, ($`\stackrel{~}{e_0}e_0`$ in Eq.(47) ):

$$\stackrel{~}{e_0}(\mathrm{}^2k_0^2/2m)(\mathrm{\hspace{0.17em}1}+\mathrm{const}_{\mathrm{new}}\zeta (v)\mathrm{\Phi }_\mathrm{o}^2+\mathrm{})\mathrm{const}_{\mathrm{new}}^{}|\mu |+\stackrel{~}{e}_{\mathrm{ph}},$$
(49)
and $`\stackrel{~}{e}_{\mathrm{ph}}3e_{\mathrm{ph}}^{(2)}5e_{\mathrm{ph}}^{(3)}`$. Second, we assume that there is a “leakage” from the moving localized state, and such a leakage occurs dynamically,

$$\delta N(t)\delta N_\mathrm{o}(t)\mathrm{const}N_\mathrm{o}\left(\delta v_{\mathrm{top}}t/L_0\right)^2<\mathrm{\hspace{0.17em}0}.$$
(50)
Then, the tail behind the condensate grows as

$$\delta \widehat{N}_{\mathrm{x},\mathrm{tail}}\delta N_\mathrm{o}(t)=\mathrm{const}^{}N_\mathrm{o}\left(\delta v_{\mathrm{top}}t/L_0\right)^2>\mathrm{\hspace{0.17em}0}.$$
(51)
To simplify the model, we disregard the creation of inside excitations, or the out-of-condensate cloud around the Bose-core, in this article. Then, $`\mathrm{const}\mathrm{const}^{}`$. This means the following terms,

$$\stackrel{~}{e_0}\mathrm{const}N_\mathrm{o}\left(\delta v_{\mathrm{top}}t/L_0\right)^2+\delta e_0N_\mathrm{o}\left(\delta v_{\mathrm{top}}t/L_0\right)^2$$
(52)
and

$$E_{\mathrm{tail}}(t)\mathrm{}\omega _\mathrm{x}\delta N_{\mathrm{x},\mathrm{tail}}(t)+\mathrm{}\omega _{\mathrm{ph}}\delta N_{\mathrm{ph},\mathrm{tail}}(t),$$
(53)
have to be included into the (total) energy of the moving packet. The time interval $`\mathrm{\Delta }t`$, during which our quasistationary solution for the core $`+`$ tail can be used to describe the transport of the total packet, can be roughly estimated from the condition $`|\delta v_{\mathrm{top}}|\mathrm{\Delta }t/L_\mathrm{o}\mathrm{\hspace{0.17em}1}`$, (see Eq. (45) and Fig. 3). To estimate the excitation energies $`\mathrm{}\omega _\mathrm{x}`$ and $`\mathrm{}\omega _{\mathrm{ph}}`$ in Eq. (53), we have to discuss the excitation spectrum of the moving condensate.

## 6 Excitations and the tail of soliton

Within the quasistationary approximation, the following decomposition of the field operators is used ,:

$$\widehat{\psi }(𝐱,t)=\mathrm{\Psi }_\mathrm{o}(x,t)+\delta \widehat{\psi }(𝐱,t),$$

$$\widehat{u}(𝐱,t)=u_\mathrm{o}(x,t)+\delta \widehat{u}(𝐱,t),$$

where the fields $`\delta \widehat{\psi }`$ and $`\delta \widehat{u}(𝐱,t)`$ describe the out-of-condensate particles. To consider the case in which the moving condensate emits excitations during the observation time, one has to introduce the fluctuation of the condensate, $`\delta \mathrm{\Psi }_\mathrm{o}(x,t)`$ and $`\delta u_\mathrm{o}(x,t)`$, and the fluctuation of the quasistationary out-of-condensate cloud described by $`\delta \widehat{\psi }^{}`$ and $`\delta \widehat{u}_j^{}`$. If the number of the “lost” particles, $`\delta N`$, is small during the observation time (but $`\delta N(t)`$ is continuously growing), one can consider the moving coherent packet as a quasistable one. Therefore, we can write the field decomposition as follows:

$$\widehat{\psi }(𝐱,t)=\mathrm{exp}\left(i\phi _c(x,t,\psi _\mathrm{o}(x,t))\right)\{\psi _\mathrm{o}(xv(x,t)t)+\delta \widehat{\psi }_0(xv(x,t)t,𝐱_{},t)\}+$$

$$+\delta \mathrm{\Psi }_\mathrm{o}(x,t)+\delta \widehat{\psi }^{}(𝐱,t),$$
(54)
$$\widehat{u}_j(𝐱,t)=u_\mathrm{o}(xv(x,t)t)\delta _{1j}+\delta \widehat{u}_{0,j}(xv(x,t)t,𝐱_{},t)+$$

$$+\delta u_\mathrm{o}(x,t)\delta _{1j}+\delta \widehat{u}_j^{}(𝐱,t).$$
(55)
In this article, we do not take into account the fluctuational parts in this decomposition and, strictly speaking, we can only estimate the the kinetic effects, such as the proposed $`_tN_\mathrm{o}(t)0`$ and $`_t\delta n(t)=_t\delta \widehat{\psi }^{}\delta \widehat{\psi }(t)0`$. However, the linear quasistationary equations on $`\delta \widehat{\psi }_0`$ and $`\delta \widehat{u}_{0,j}`$ can be used to find the excitation spectrum of the outside excitations that is required for Eq. (53). Then, assuming the quasistability of the Bose-core during the observation time $`\mathrm{\Delta }t`$, one can switch on the mechanism of occupancy redistribution between the macroscopically occupied mode $`\mathrm{\Psi }_\mathrm{o}(x,t,N_\mathrm{o})`$ and the out-of-condensate excitations, e.g., $`\delta n_k=\delta \widehat{\psi }_k^{}\delta \widehat{\psi }_k`$.

Thus, we “dress” the quasistationary parts of the out-of-condensate fields by Eqs. (54) and (55) and obtain the following set of equations ($`xv(x,t)tx`$):

$$i\mathrm{}_t\delta \widehat{\psi }_0(x,𝐱_{},t)=|\mu |\delta \widehat{\psi }_0(𝐱,t)\frac{\mathrm{}^2}{2m}\mathrm{\Delta }\delta \widehat{\psi }_0(𝐱,t)\frac{\mathrm{}^2}{2m}\zeta (v)\psi _\mathrm{o}^2(x)\mathrm{\Delta }\delta \widehat{\psi }_0(𝐱,t)+$$

$$+\left\{(\nu _0+\stackrel{~}{\nu _0}+|\mu |\zeta (v))\psi _\mathrm{o}^2(x)+(2\nu _1+\stackrel{~}{\nu _1})\psi _\mathrm{o}^4(x)\right\}\delta \widehat{\psi }_0(x,𝐱_{},t)+$$

$$+\left\{\nu _0\psi _\mathrm{o}^2(x)+2\nu _1\psi _\mathrm{o}^4(x)\right\}\delta \widehat{\psi }_0^{}(𝐱,t)+\sigma _{\mathrm{eff}}\psi _\mathrm{o}(x)\delta \widehat{𝐮}_0(x,𝐱_{},t),$$
(56)
$$\left(\left\{_tv(x,t)_x\right\}^2c_l^2\mathrm{\Delta }\right)\delta \widehat{u}_{0,x}(x,𝐱_{},t)$$

$$c_l^22\kappa _3_xu_\mathrm{o}(x)_x^2\delta \widehat{u}_{0,x}(x,𝐱_{},t)c_l^22\kappa _3_x^2u_\mathrm{o}(x)_x\delta \widehat{u}_{0,x}(x,𝐱_{},t)=$$

$$=\rho ^1\stackrel{~}{\sigma }_{\mathrm{eff}}_x(\psi _\mathrm{o}(x)\delta \widehat{\psi }_0(x,𝐱_{},t)+\mathrm{h}.\mathrm{c}.).$$
(57)
Here, we simplified the boson-phonon coupling terms and, roughly, $`\sigma _{\mathrm{eff}}\stackrel{~}{\sigma }_{\mathrm{eff}}\sigma _0`$. Note that the interaction between $`\delta \widehat{\psi }_0`$ and $`\delta \widehat{u}_{0,x}`$ is mediated through the condensate $`\psi _\mathrm{o}(x,t)_xu_\mathrm{o}(x,t)`$ in the linear theory at $`T0`$.

For the outside excitations, the first approximation of Eqs. (56) and (57) consists of taking uncoupled excitons and phonons. However, they live in a half of the medium ($`x<0`$) and form a tail behind the condensate. Then, we have the excitonic excitations with

$$\mathrm{}\omega _\mathrm{x}(k_x)=|\mu |+|\mu |(k_xL_0)^2,|k_x|L_0<1,$$

$$\delta \widehat{\psi }(x,t)\left(\sqrt{V/2}\right)^1\mathrm{exp}(i\phi _k+ik_xxi\omega _\mathrm{x}t)\widehat{a}_k,x<0.$$

This is exactly a free exciton with the energy $`\stackrel{~}{E}_g+\mathrm{}^2(k_0|k_x|)^2/2m`$ in the laboratory frame. It moves behind the condensate with $`p_x=\mathrm{}(k_0|k_x|)>0`$. For the acoustic phonons, we have ($`k_x^{}>0`$):

$$\mathrm{}\omega _{\mathrm{ph}}(k_x^{})=\mathrm{}\omega _{\mathrm{ac}}(k_x^{})v\mathrm{}k_x^{},$$

$$\delta \widehat{u}_x(x,𝐱_{},t)\mathrm{const}(k^{})\mathrm{exp}\left(i\phi _k^{}+ik_x^{}xi\{\omega _{\mathrm{ac}}(k_x^{})k_x^{}v\}t\right)\widehat{b}_k^{}+\mathrm{h}.\mathrm{c}.,$$

where $`\mathrm{const}(k^{})\left(\sqrt{\rho (V/2)\omega _{\mathrm{ac}}(k^{})}\right)^1`$. (After $`xxvt`$, we have free phonons in the tail, and $`k_x^{}`$ is the same in the laboratory frame.) The ‘resonance’ equation for the excitations written in the comoving frame has the following form:

$$|\mu |+|\mu |(k_xL_0)^2=\mathrm{}(c_lv)k_x^{},$$
(58)
and it can be solved easily. Here, we use the following estimate for the ‘relevant’ phonons:

$$\left(k_x^{}L_0\right)0.5\frac{\sqrt{|\mu |/mc_l^2/2}}{(1v/c_l)}.$$

For $`k_x^{}>0`$, we have $`\left(k_x^{}L_0\right)(110)\sqrt{|\mu |/mc_l^2/2}`$. In this article, we do not go into detail of how the coupling between the outside excitons and phonons is formed within Eqs. (56), (57). Our aim is to show that it is possible to balance the energy of the total packet by taking into account the growing out-of-condensate tail (51) behind the coherent state (12), (13) with leakage (50).

Returning to Eq. (53), we estimate the characteristic energies of excitations as follows:

$$\mathrm{}\omega _\mathrm{x}\mathrm{}^2\stackrel{~}{k}_0^2/2m\mathrm{and}\mathrm{}\omega _{\mathrm{ph}}\mathrm{const}_k|\mu |,\mathrm{const}_k28.$$
(59)
(The value of $`\mathrm{}\omega _{\mathrm{ph}}`$ comes from the resonance condition at $`k_{\mathrm{ph},x}>0`$.) We assume that the tail consists of excitons and phonons, and the selfconsistency condition between $`\psi _\mathrm{o}^2(x/L_0)\mathrm{\Phi }_\mathrm{o}^2`$ and, roughly, $`_xu_\mathrm{o}(x/L_0)\mathrm{\Phi }_\mathrm{o}^2`$ leads to the coherent emission of the outside excitons and phonons, and $`\mathrm{\Phi }_\mathrm{o}^2\mathrm{\Phi }_\mathrm{o}^2\delta \mathrm{\Phi }_\mathrm{o}^2(t)`$. Then, we have

$$\delta \widehat{N}_{\mathrm{x},\mathrm{tail}}\delta \widehat{N}_{\mathrm{ph},\mathrm{tail}}.$$

As a result, the energy of the total packet

$$E(t)e_\mathrm{o}N_\mathrm{o}(t)+\delta e_0N_\mathrm{o}(t)\left(\delta v_{\mathrm{top}}(t)t/L_0(t)\right)^2\stackrel{~}{e_0}\mathrm{const}N_\mathrm{o}(t)\left(\delta v_{\mathrm{top}}(t)t/L_0(t)\right)^2+$$

$$+\mathrm{}\omega _\mathrm{x}\mathrm{const}^{}N_\mathrm{o}(t)\left(\delta v_{\mathrm{top}}(t)t/L_0(t)\right)^2+\mathrm{}\omega _{\mathrm{ph}}\mathrm{const}^{}N_\mathrm{o}(t)\left(\delta v_{\mathrm{top}}(t)t/L_0(t)\right)^2$$
(60)
is conserved if one takes $`\stackrel{~}{e_0}\mathrm{const}\delta e_0`$ and $`\mathrm{const}\mathrm{const}^{}0.1`$ in Eqs. (52), (53), and (60) (the last value is probably overestimated).

To address the problem of the dynamic stability, one has to estimate how many particles and phonons are emitted out from the moving coherent packet. We start from the same initial conditions as in Fig. 2. However, the occupations of the tail states start growing immediately as $`t>0`$ while the occupancy of the coherent state starts decreasing. Using the conditions $`E(t)=E_\mathrm{o}(t=0)`$ and $`N_\mathrm{o}(t)+\delta N_{\mathrm{tail}}(t)=N_\mathrm{o}(t=0)`$, we obtain the dynamics presented on Fig. 3 by combination of numerics and analytic estimates. For example, after the transport has been observed during the time interval of $`\mathrm{\Delta }t`$ corresponding to $`\delta v_{\mathrm{top}}\mathrm{\Delta }t/L_0(1/4\mathrm{\hspace{0.17em}1}/2)`$, the soliton develops the leakage of $`\delta N(13)\mathrm{\hspace{0.17em}10}^2N_\mathrm{o}`$ into the tail excitations, see Fig. 3. However, the value of $`\delta v_{\mathrm{top}}(N_\mathrm{o})\mathrm{\Delta }t/L_0(N_\mathrm{o})`$ starts to decrease because of this leakage too. For instance, when the adjusted value of $`\delta v_{\mathrm{top}}\left(N_\mathrm{o}(t)\right)\mathrm{\Delta }t^{}/L_0\left(N_\mathrm{o}(t)\right)1`$, the leakage is esimated as $`\delta N(12)\mathrm{\hspace{0.17em}10}^1N_\mathrm{o}`$, see Fig. 3. Therefore, the condensate can be considered as a quasistable one during the time interval of $`\mathrm{\Delta }t^{}`$. The following formulas can be used for estimates

$$\mathrm{\Phi }_\mathrm{o}^2(t)S_{}\mathrm{\hspace{0.17em}2}L_0(t)N_\mathrm{o}(t),$$
(61)
$$\delta n_{\mathrm{tail}}(t)S_{}\stackrel{~}{L}q(t)N_\mathrm{o},$$
(62)
where $`\delta n_{\mathrm{tail}}(t)q^{}\mathrm{\Phi }_\mathrm{o}^2(t)`$, $`\stackrel{~}{L}=\stackrel{~}{q}L_0`$. If $`q(t)1`$ (e.g., $`q10^1`$ within the time interval $`\mathrm{\Delta }t^{}`$ on Fig. 3) and $`\stackrel{~}{q}1`$ (e.g., $`\stackrel{~}{q}1050`$), one can roughly estimate the dimensionless $`q^{}`$ as $`10^1q`$. The dynamics presented on Fig. 3 is in a qualitative agreement with experimental data . (It is interesting to mention that the similar problems of dynamic stability appear in the theory of the so-called Embedded Solitons in nonlinear optics .)

## 7 Conclusion

When a quasistationary exciton-phonon soliton moves in a periodic medium, an asymmetric form of the ballistic signal can be developed as a dynamic effect due to exciton-phonon interaction. This effect has the tendency to be accumulated with time, and it can be clearly seen after the packet has traveled the distance of $`nL_0`$, $`n1`$. Here, the width of the packet is $`(36)L_0`$, and the value of $`n`$ depends on the strength of relevant interaction parameters, such as $`\vartheta _0`$ and $`\sigma _0`$. It can be estimated by the following formula (see Eq. (44) ):

$$n\left(|\delta v_{\mathrm{top}}|/v\right)^110^310^2.$$

To a first approximation, we introduced the coherent part of the packet and described its dynamics by use of the selfconsistent generalization of $`v=\mathrm{const}v(x,t)`$ and $`\phi _\mathrm{c}(x,t)=k_0x\omega _0t\stackrel{~}{\phi }_\mathrm{c}(x,t)`$ in such a way:

$$v\stackrel{\vartheta _00}{}v(_xu_\mathrm{o}(x,t))v(|\mathrm{\Psi }_\mathrm{o}(x,t)|^2)$$

$$\phi _\mathrm{c}(x,t)\stackrel{\vartheta _00}{}\stackrel{~}{\phi }_\mathrm{c}(x,t,_xu_\mathrm{o}(x,t))\phi _\mathrm{c}(x,t,|\mathrm{\Psi }_\mathrm{o}(x,t)|^2).$$

This substitute was applied within the quasistationary approximation to simplify and solve the corresponding dynamic equations on the coherent Bose-core of the total packet. The incoherent part of it turned out to be equally important in dynamics of the total packet. Indeed, the quasistationary solution for the Bose-core cannot move in a periodic medium as a single soliton, but with the leakage into the out-of-condensate excitation states. Thus, the question on whether a moving conserving solution (i.e., a soliton conserving at least $`E_\mathrm{o}`$, $`P_{\mathrm{o},x}`$, and $`N_\mathrm{o}`$) exists in the model with $`\vartheta _00`$ and $`\sigma _00`$ remains open. In this article, we showed that an exciton-phonon “comet” is a better image for the moving ballistic packet, and the core of it (a kind of the “nucleus” of such a comet) can be modeled by a coherent state, or a condensate.

The ansatz we used for the soliton-like solution can be applied within the validity of the adiabatic hypothesis (i.e., it is not the exact one of the dynamic equations). As a result, we could obtain an explicit dependence on time ( $`\left(\delta v_{\mathrm{top}}t/L_0\right)^2`$ ) for the exact energy and momentum of the moving soliton if it were the only one component of the moving packet. However, the technique of Bose-Einstein condensation allows to proceed with the approximate solution obtained in this article. It is known that the interaction between Bose-condensate and non-condensed particles leads to an effective dissipation for the condensate wave function at $`T0`$ . We assume that an effective dissipation appears naturally in description of the transport properties of the coherent part, $`\mathrm{\Psi }_\mathrm{o}(x,t)`$ and $`_xu_\mathrm{o}(x,t)`$, at $`T0`$. If the initial state of the packet was taken as a pure coherent state with a nonzero momentum, the moving condensate emits excitations, i.e., $`N_\mathrm{o}N_\mathrm{o}\delta N(t)`$ during the observation time. Then, the out-of-condensate cloud of collective excitations and, in particular, the tail consisting of the out-of-condensate excitons and phonons grow around and behind the moving coherent state, e.g., $`\delta \widehat{n}(x,t)0`$, $`x<3L_0`$.

We argue that the emission effect can make the dynamics of the total moving packet conservative. Therefore, the coherent part of such a packet has to be considered as a quasistable core (a “nucleus” of the “comet”) of the moving exciton-phonon droplet provided $`\delta N_{\mathrm{tail}}(t)N_\mathrm{o}(t)`$ during the observation time. In addition, the coherent phase $`\phi _\mathrm{c}(x,t)`$ prescribed to the macroscopic wave function $`\mathrm{\Psi }_\mathrm{o}(x,t)`$ cannot be taken as a regular field under such conditions because of dephasing effects . Its fluctuations has to be taken into account to clarify the coherent properties of the Bose-core. On the other hand, the rigorous approach to the dynamics of solitons with emission lies beyond the quasistationary approximation we used in this articles. A set of kinetic equations has to be applied to treat this problem in detail.

## 8 Acknowledgements

One of the authors (D.R.) thanks I. Loutsenko for critical reading of the manuscript.


format: text

text:
# Efficient schemes for reducing imperfect collective decoherences

## Abstract

We propose schemes that are efficient when each pair of qubits undergoes some imperfect collective decoherence with different baths. In the proposed scheme, each pair of qubits is first encoded in a decoherence-free subspace composed of two qubits. Leakage out of the encoding space generated by the imperfection is reduced by the quantum Zeno effect. Phase errors in the encoded bits generated by the imperfection are reduced by concatenation of the decoherence-free subspace with either a three-qubit quantum error correcting code that corrects only phase errors or a two-qubit quantum error detecting code that detects only phase errors, connected with the quantum Zeno effect again.

Information processing with quantum bits (qubits), e.g., quantum computing, quantum cryptography, and quantum gambling is, a novel technique that solves some classically intractable problems -. However, in order to make quantum information processing involving many qubits practical, some methods for reducing decoherence (MRDs) are indispensable. Among these, there are quantum error correcting codes (QECCs) -, decoherence-free subspaces (DFSs) -, the quantum Zeno effect (QZE) , ***The QZE was discovered by Misra and Sudersan . The use of the QZE for combating decoherence was first suggested by Zurek , and it is a part of a scheme considered by Barenco et al. . and dynamical suppression of decoherence .

If the provisos for DFSs are fulfilled, DFSs are more efficient than QECCs or the QZE in respect to the amount of other necessary resources as well as the number of qubits. The robustness of DFSs against perturbation of replica symmetry is shown in Refs. . Indeed, qubits in collective decoherence with imperfect replica symmetry can be preserved with concatenation of DFSs with QECCs . However, the efficiencies of various MRD’s depend on the decoherence model. So devising an optimal scheme that appropriately combines existing MRDs for a given decoherence model will be important in the design of quantum information processors. In this paper, we propose a scheme that is efficient when each pair of qubits undergoes imperfect collective decoherence with different baths (cluster decoherence ). We start with a subspace composed of two qubits which is decoherence-free against a certain interaction that generates only phase errors. Other interactions assumed small but non-negligible make the encoded states leak out of the DFS, and the interactions generate phase errors in the encoded qubits. The leakage is reduced by the QZE. The phase errors in the encoded qubits are corrected by concatenating the DFS with a three-qubit QECC that corrects only phase errors  or by concatenating the DFS with a two-qubit quantum error detecting code that detects only phase errors and by QZE again.

The dynamics of the qubits and bath is governed by

$`𝐇_T=𝐇_S+𝐇_B+𝐇_I,`$ (1)
where $`𝐇_T`$, $`𝐇_S`$, and $`𝐇_B`$ denote the total, the system, and the bath (or environment) Hamiltonian, respectively, and $`𝐇_I`$ is the interaction Hamiltonian. First, we consider the following simple model.

$`𝐇_S`$ $`=`$ $`ϵ(\sigma _1^z+\sigma _2^z),`$ (2)
$`𝐇_I`$ $`=`$ $`\lambda (\sigma _1^z+\sigma _2^z)V_z,`$ (3)
where $`\sigma _i^z`$ ($`i=1,2`$) are Pauli spin operators, $`V_z`$ is the bath operator coupled to the degree of freedom, and $`𝐇_B`$ is arbitrary. The type of Hamiltonian in Eq. (3), which corresponds to a special case of the spin-boson problem, has been used by many authors to model decoherence despite its simplicity : This model describes a decohering mechanism with only phase errors. Amplitude errors would involve a time scale much longer than that of phase errors in some real physical systems . (Later, we will treat more general models.) We can easily see that a subspace spanned by $`|01`$ and $`|10`$ satisfies the DFS condition in the case of the interaction given by Eq. (3): $`(\sigma _1^z+\sigma _2^z)|01=0|01`$ and $`(\sigma _1^z+\sigma _2^z)|10=0|10`$. It is noted that another condition should be additionally satisfied in order that some subspaces become decoherence-free: the system Hamiltonian $`𝐇_S`$ does not make qubits leak out of the subspace. Otherwise, $`𝐇_S`$ needs to be eliminated to satisfy this condition by the method proposed in Ref. . Therefore,

$`𝐇_T[(\alpha |01+\beta |10)`$ $`|\mathrm{\Psi }_b(0)]`$ (4)
$`=(\alpha |01+\beta |10)`$ $`[0(ϵ+\lambda V_z)+𝐇_B]|\mathrm{\Psi }_b(0),`$ (5)
and as a result

$`\mathrm{exp}[i𝐇_TT_0](\alpha |01+\beta |10)|\mathrm{\Psi }_b(0)`$ (6)
$`=(\alpha |01+\beta |10)\mathrm{exp}[i𝐇_BT_0]|\mathrm{\Psi }_b(0).`$ (7)
Here we can see that the qubits indeed do not decohere. So Span$`[|01,|10]`$ (the subspace that $`|01`$ and $`|10`$ span) can be used to encode one qubit. That is, we can encode a qubit $`\alpha |0+\beta |1`$ into, for example, the state

$$|\mathrm{\Psi }_{enc}=\alpha |01+\beta |10.$$
(8)
It is clear that this subspace is sufficient for preventing decoherence provided that the system and bath are perfectly governed by Eq. (3). However, in real systems there are some small perturbative interactions that are not included in Eq. (3). When the perturbative interaction is non-negligible, its effect must be reduced by some method that we will describe. Now, let us consider a more general decoherence model:

$`𝐇_S`$ $`=`$ $`ϵ_1\sigma _1^z+ϵ_2\sigma _2^z,`$ (9)
$`𝐇_I`$ $`=`$ $`[(\lambda _1^z\sigma _1^z+\lambda _2^z\sigma _2^z)V_z+(\lambda _1^+\sigma _1^++\lambda _2^+\sigma _2^+)V_++`$ (11)
$`(\lambda _1^{}\sigma _1^{}+\lambda _2^{}\sigma _2^{})V_{}].`$
Here, $`\sigma _i^j`$ ($`j=z,+,`$) are Pauli spin operators and $`V_j`$ are the bath operators coupled to these degrees of freedom. We assume that $`\mathrm{\Delta }ϵϵ_2ϵ_1ϵ_1`$ and $`\mathrm{\Delta }\lambda ^z\lambda _2^z\lambda _1^z\lambda _1^z`$. We also assume that phase damping is dominant $`\lambda _i^z\lambda _i^+`$ and $`\lambda _i^z\lambda _i^{}`$. In the limit when $`\mathrm{\Delta }ϵ`$, $`\mathrm{\Delta }\lambda ^z`$, $`\lambda _i^+`$, and $`\lambda _i^{}`$ vanish, Eq. (11) reduces to Eq. (3). Then let us consider the following. In a short period of time $`T_0/N`$, under the Hamiltonian $`𝐇_T`$, the encoded state evolves into

$`|\mathrm{\Psi }(T_0/N)`$ (12)
$``$ $`[1i𝐇(T_0/N)](\alpha |01+\beta |10)|\mathrm{\Psi }_b(0)`$ (13)
$`=`$ $`(\alpha |01+\beta |10)`$ (18)
$`[10i(T_0/N)(ϵ_1+\lambda _1^zV_z)i(T_0/N)𝐇_B]|\mathrm{\Psi }_b(0)`$
$`+(T_0/N)(\alpha |01+\beta |10)(\mathrm{\Delta }ϵ+\mathrm{\Delta }\lambda ^zV_z)|\mathrm{\Psi }_b(0)`$
$`+(T_0/N)|00(\lambda _1^+\beta +\lambda _2^+\alpha )V_+|\mathrm{\Psi }_b(0)`$
$`+(T_0/N)|11(\lambda _1^{}\alpha +\lambda _2^{}\beta )V_{}|\mathrm{\Psi }_b(0)],`$
where $`|\mathrm{\Psi }_b(0)`$ denotes the bath state. Then we perform a measurement that discriminates between the encoding space Span$`\{|01,|10\}`$ and Span $`\{|00,|11\}`$. This measurement can be implemented by XORing each qubit to an ancilla qubit consecutively . By frequently (i.e., $`N`$ is made large) repeating time evolution by the Hamiltonian in Eq. (18) and the consecutive measurements, we can make effects of the terms involving $`|00`$ and $`|11`$ in Eq. (18) negligible (QZE). Then, after some simple calculations, we obtain

$`|\mathrm{\Psi }(T_0)`$ $``$ $`(\alpha |01+\beta |10)|\mathrm{\Psi }_b`$ (20)
$`+O(T_0)(\alpha |01+\beta |10)|\mathrm{\Psi }_b^{},`$
where $`|\mathrm{\Psi }_b`$ and $`|\mathrm{\Psi }_b^{}`$ are some arbitrary bath states which are not necessarily orthogonal to each other. We can see that overall time evolution generates only phase errors in the encoded bit. In other words, the QZE prevents leakage of the states out of the encoding space while does not prevent time evolution within the encoding space. However, the QZE can be practicable for only fairly stable quantum states. Typical systems are well described by the model assumed here and thus the encoded qubit is fairly stable. Therefore, in this case the QZE can be a suitable choice for protecting the encoded qubit. When $`\mathrm{\Delta }ϵ`$ and $`\mathrm{\Delta }\lambda ^zV_z`$ are negligible, the second term of the right hand side of Eq. (20) is negligible and thus we need no more MRD. The two-qubit DFS in Eq. (8) plus the QZE is sufficient for preservation of one qubit. When they are not, we should reduce the effect of the term. This can be done in two ways, as noted in the Introduction. First, we concatenate the DFS in Eq. (8) with a three-qubit QECC that corrects only the phase errors (Eq. (15) of ). In this case, six qubits are needed to encode one qubit in the proposed scheme. Secondly, we concatenate the DFS in Eq.(8) with a two-qubit quantum code that detects only phase errors. That is,

$`|0_{enc}`$ $`=`$ $`(|0+|1)(|0+|1),`$ (21)
$`|1_{enc}`$ $`=`$ $`(|0|1)(|0|1),`$ (22)
where $`|0`$ and $`|1`$ denote encoded qubits using the two-qubit DFS in Eq. (8) and normalization factors are omitted. Then we preserve the states using the QZE again: We frequently perform measurements that tell us whether the error has occurred or not . In this case, four qubits are needed to encode one qubit in the proposed scheme. The first proposed scheme (two-qubit DFS $`+`$ QZE) is efficient when replica asymmetry is negligible ($`\mathrm{\Delta }\lambda _z0`$ and $`\mathrm{\Delta }ϵ0`$) and other terms ($`\lambda _i^+`$ and $`\lambda _i^{}`$) are small but non-negligible. The second (\[two-qubit DFS $`+`$ QZE\]$`\times `$ three-qubit QECC) and third (\[two-qubit DFS $`+`$ QZE\]$`\times `$ \[two-qubit quantum error detecting code $`+`$ QZE\]) proposed schemes are efficient when replica asymmetry is also non-negligible.

In Duan and Guo’s scheme , the subspace that is orthogonal to the space to which the subspaces leak through the interaction Hamiltonian is adopted as the encoding space. Then the QZE is used for preventing leakage of qubits out of the encoding space. In contrast, in our scheme the qubit is first stabilized using DFS and then leakage is prevented by the QZE and time evolution within the encoding space is corrected by other MRDs (QECC or quantum error detecting code plus QZE). So the encoding space of Duan and Guo’s scheme differs from that of our scheme for a given Hamiltonian. For example, in the case of the model of Eq. (3), the encoding space of Duan and Guo’s scheme is Span$`[|\overline{0}\overline{1}`$, $`|\overline{1}\overline{0}]`$ where $`|\overline{0}=(1/\sqrt{2})(|0+|1)`$ and $`|\overline{1}=(1/\sqrt{2})(|0|1)`$. This differs from the encoding space Span$`[|01,|10]`$ of our scheme. Duan and Guo’s scheme  is more powerful than ours in that theirs is effective for wide classes of decoherence, i.e., for independent and even cooperative decoherence. In contrast, our scheme is a specialized one that is efficient in the case where phase errors are dominant but other errors are still non-negligible.

Here we proposed three schemes that are efficient when each pair of qubits undergoes some imperfect collective decoherence with different baths. In the first scheme, each pair of qubits is encoded in a DFS composed of two qubits. Leakage out of the encoding space generated by the imperfection is reduced by the quantum Zeno effect. In the second scheme, phase errors in the encoded bits also generated by imperfection of replica symmetry are reduced by concatenation of the DFS with a three-qubit QECC that corrects only phase errors. In the third scheme, the same thing is done by concatenation of the DFS with two-qubit quantum error detecting code that detects only phase errors plus the QZE again.

###### Acknowledgements.

This work was supported by the Korean Ministry of Science and Technology through the Creative Research Initiatives Program under Contract No. 98-CR-01-01-A-20.


text:
# 1 Introduction

## 1 Introduction

We consider the following problem.

#### Problem E.

Let $`H_n=\underset{j=1}{\overset{n}{}}\frac{1}{j}`$. Show that, for each $`n1`$,

$$\underset{d|n}{}dH_n+\mathrm{exp}(H_n)\mathrm{log}(H_n),$$
(1.1)
with equality only for $`n=1`$.

The function $`\sigma (n)=_{d|n}d`$ is the sum of divisors function, so for example $`\sigma (6)=12`$. The number $`H_n`$ is called the n-th harmonic number by Knuth, Graham and Patashnik \[9, Sect. 6.3\], who detail various properties of harmonic numbers.

The ‘$`E`$’ in Problem $`E`$ might stand for either ‘easy’ or ‘elementary’. Perhaps ‘$`H`$’ for ‘hard’ would be a better letter to use, since our object is to show the following equivalence.

###### Theorem 1.1

Problem $`E`$ is equivalent to the Riemann hypothesis.

The Riemann hypothesis, stated by Riemann  in 1859, concerns the complex zeros of the Riemann zeta function. The Riemann zeta function $`\zeta (s)`$ is defined by the Dirichlet series

$$\zeta (s)=\underset{n=1}{\overset{\mathrm{}}{}}n^s,$$

which converges for $`\mathrm{}(s)>1`$, and it has an analytic continuation to the complex plane with one singularity, a simple pole with residue $`1`$ at $`s=1`$. The Riemann hypothesis states that the nonreal zeros of the Riemann zeta function $`\zeta (s)`$ all lie on the line $`\mathrm{}(s)=\frac{1}{2}.`$ One reason for the great interest in the Riemann hypothesis, regarded by many as the most important unsolved problem in pure mathematics, is its connection with the distribution of prime numbers, described below. More significantly, the Riemann hypothesis is a special case of questions concerning generalizations of the zeta function ($`L`$-functions) and their connections with problems in number theory, algebraic geometry, topology, representation theory and perhaps even physics, see Berry and Keating , Katz and Sarnak  and Murty.

The connection of the Riemann hypothesis with prime numbers was the original question studied by Riemann . Let $`\pi (x)`$ count the number of primes $`p`$ with $`1<px`$. C. F. Gauss noted empirically that $`\pi (x)`$ is well approximated by the logarithmic integral

$$Li(x)=_2^x\frac{dt}{\mathrm{log}t},$$

which itself satisfies

$$Li(x)=\frac{x}{\mathrm{log}x}+O(\frac{x}{(\mathrm{log}x)^2}).$$

The Riemann hypothesis is equivalent to the assertion that for each $`ϵ>0`$ there is a positive constant $`C_ϵ`$ such that

$$|\pi (x)Li(x)|C_ϵx^{1/2+ϵ}$$

for all $`x2`$, see Edwards \[5, p. 90\]. The force of the Riemann hypothesis lies in the small size of the error term. The strongest form known of the Prime Number Theorem with error term asserts that

$$|\pi (x)Li(x)|C_1x\mathrm{exp}(C_2(\mathrm{log}x)^{3/5ϵ}),$$

for any positive $`ϵ`$, for certain positive constants $`C_1`$ and $`C_2`$ depending on $`ϵ`$; this result is due to Vinogradov and Korobov in 1958.

Problem $`E`$ encodes a modification of a criterion of Guy Robin  for the Riemann hypothesis. Robin’s criterion states that the Riemann hypothesis is true if and only if

$$\sigma (n)<e^\gamma n\mathrm{log}\mathrm{log}n\text{for all}n5041,$$
(1.2)
where $`\gamma 0.57721`$ is Euler’s constant. This criterion is related to the density of primes, as explained in §2. Our aim was to obtain a problem statement as elementary as possible, containing no undefined constants. However the hard work underlying the equivalence resides in the results of Robin stated in §3, where we give a proof of Theorem 1.1.

Before coming to the proof, in the next section we describe how the Riemann hypothesis is related to the sum of divisors function. The connection traces back to Ramanujan’s work on highly composite numbers, and involves several results of Erdős with coauthors.

## 2 Colossally Abundant Numbers

The Riemann hypothesis is encoded in the criterion of Theorem 1.1 in terms of the very thin set of values of $`\sigma (n)`$ that are “large.” The sum of divisors function is given by

$$\sigma (n)=\underset{p^a||n}{}(1+p+p^2+\mathrm{}+p^a)=n\underset{p^a||n}{}(1+\frac{1}{p}+\mathrm{}+\frac{1}{p^a}),$$
(2.3)
where the product is taken over all primes $`p`$ dividing $`n`$ and the notation $`p^a||n`$, means $`p^a`$ divides $`n`$ but $`p^{a+1}`$ does not divide $`n`$. Most values of $`\sigma (n)`$ are on the order of $`Cn`$. The average size of $`\sigma (n)`$ was essentially found by Dirichlet in 1849, and is given by the following result of Bachmann, cf. Hardy and Wright \[7, Theorem 324\].

###### Theorem 2.1

(Bachmann) The average order of $`\sigma (n)`$ is $`\frac{\pi ^2}{6}n.`$ More precisely,

$$\underset{j=1}{\overset{n}{}}\sigma (j)=\frac{\pi ^2}{12}n^2+O(n\mathrm{log}n)$$

as $`n\mathrm{}.`$

The maximal order of $`\sigma (n)`$ is somewhat larger, and was determined by Gronwall in 1913, see Hardy and Wright \[7, Theorem 323, Sect. 18.3 and 22.9\].

###### Theorem 2.2

(Gronwall) The asymptotic maximal size of $`\sigma (n)`$ satisfies

$$\underset{n\mathrm{}}{lim\; sup}\frac{\sigma (n)}{n\mathrm{log}\mathrm{log}n}=e^\gamma ,$$

where $`\gamma `$ is Euler’s constant.

This result can be deduced from Mertens’ theorem, which asserts that

$$\underset{px}{}(1\frac{1}{p})\frac{e^\gamma }{\mathrm{log}x}.$$

as $`x\mathrm{}.`$ A much more refined version of the asymptotic upper bound, due to Robin \[18, Theorem 2\], asserts (unconditionally) that for all $`n3`$,

$$\sigma (n)<e^\gamma n\mathrm{log}\mathrm{log}n+0.6482\frac{n}{\mathrm{log}\mathrm{log}n}.$$
(2.4)
In §3 we will show that

$$H_n+\mathrm{exp}(H_n)\mathrm{log}(H_n)e^\gamma n\mathrm{log}\mathrm{log}n+\frac{7n}{\mathrm{log}n},$$

the right side of which is only slightly smaller than that of (2.4). Comparing this bound with (2.4) shows that the inequality (1.1), if ever false, cannot be false by very much.

The study of extremal values of functions of the divisors of $`n`$ is a branch of number theory with a long history. Let $`d(n)`$ count the number of divisors of $`n`$ (including $`1`$ and $`n`$ itself). Highly composite numbers are those positive integers $`n`$ such that

$$d(n)>d(k)\text{for}1kn1.$$

Superior highly composite numbers are those positive integers for which there is a positive exponent $`ϵ`$ such that

$$\frac{d(n)}{n^ϵ}\frac{d(k)}{k^ϵ}\text{for all}k>1,$$

so that they maximize $`\frac{d(n)}{n^ϵ}`$ over all $`n`$; these form a subset of the highly composite numbers. The study of these numbers was initiated by Ramanujan. One can formulate similar extrema for the sum of divisors function. Superabundant numbers are those positive integers $`n`$ such that

$$\frac{\sigma (n)}{n}>\frac{\sigma (k)}{k}\text{for}1kn1.$$

Colossally abundant numbers are those numbers $`n`$ for which there is a positive exponent $`ϵ`$ such that

$$\frac{\sigma (n)}{n^{1+ϵ}}\frac{\sigma (k)}{k^{1+ϵ}}\text{for all}k>1.$$

so that $`n`$ attains the maximum value of $`\frac{\sigma (k)}{k^{1+ϵ}}`$ over all $`k`$. The set of colossally abundant numbers is infinite. They form a subset of the superabundant numbers, a fact which can be deduced from the definition. Table 1 gives the colossally abundant numbers up to $`10^{18}`$, as given in . Robin showed that if the Riemann hypothesis is false, there will necessarily exist a counterexample to the inequality (1.2) which is a colossally abundant number, cf. \[18, Proposition 1 of Section 3\]; the same property can be established for counterexamples to (1.1). (There could potentially exist other counterexamples as well.)

Superabundant and colossally abundant numbers were studied111Alaoglu and Erdős  use a slightly stronger definition of colossally abundant number; they impose the additional requirement that $`\frac{\sigma (n)}{n^{1+ϵ}}>\frac{\sigma (k)}{k^{1+ϵ}}`$ must hold for $`1k<n`$. With their definition the colossally abundant numbers are exactly those given by (2.5) below. in detail by Alaoglu and Erdős  in 1944. As evidenced in the table, colossally abundant numbers consist of a product of all the small primes up to some bound, with exponents which are nonincreasing as the prime increases. The values of these exponents have a characteristic smooth shape which can be almost completely described (see (2.5) below). In fact these classes of numbers had been studied earlier, by Ramanujan, in his 1915 work on highly composite numbers . The notes in Ramanjuan’s Collected Papers report: “The London Mathematical Society was in some financial difficulty at the time, and Ramanujan suppressed part of what he had written in order to save expense.” \[14, p. 339\]. Only the first 52 of 75 sections were printed. The manuscript of the unpublished part was eventually rediscovered among the papers of G. N. Watson, in “Ramanujan’s Lost Notebook” (see Andrews  and Ramanujan \[15, pp. 280–308\]) and finally published in 1997, in . In it superabundant and colossally abundant numbers are considered in Section 59, as special cases of the concepts of generalised highly composite numbers and superior generalised highly composite numbers, for the parameter value $`s=1`$. Ramanujan derived upper and lower bounds for the maximal order of generalised highly composite numbers in Section 71, assuming the Riemann hypothesis. (The Riemann hypothesis was assumed from Section 40 onward in his paper.) His bounds imply, assuming the Riemann hypothesis, that (1.2) holds for all sufficiently large $`n`$.

The results of Alaoglu and Erdős in their 1944 paper are unconditional, and mainly concerned the exponents of primes occurring in highly composite and superabundant numbers. In considering the exponents of primes appearing in colossally abundant numbers, they raised the following question, which is still unsolved.

#### Conjecture.

(Alaoglu and Erdős) If $`p`$ and $`q`$ are both primes, is it true that $`p^x`$ and $`q^x`$ are both rational only if $`x`$ is an integer?

This conjecture would follow as a consequence222For non-rational $`x`$ consider $`a_1=\mathrm{log}p`$, $`a_2=\mathrm{log}q`$, $`b_1=1`$ and $`b_2=x`$ in the four exponentials conjecture. For non-integer rational $`x`$ a direct argument is used. of the four exponentials conjecture in transcendental number theory, which asserts that if $`a_1,a_2`$ form a pair of complex numbers, linearly independent over the rationals $``$, and $`b_1,b_2`$ is another such pair, also linearly independent over the rationals, then at least one of the four exponentials $`\{e^{a_ib_j}:i,j=1\text{or}2\}`$ is transcendental, see Lang \[10, pp. 8-11\]. Alaoglu and Erdős \[1, Theorem 10\] showed that for “generic” values of $`ϵ`$ there is exactly one maximizing integer $`n`$ , and the exponent $`a_p(ϵ)`$ of each prime $`p`$ in it is

$$a_p(ϵ)=\frac{\mathrm{log}(p^{1+ϵ}1)\mathrm{log}(p^ϵ1)}{\mathrm{log}p}1;$$
(2.5)
furthermore for all $`ϵ>0`$ this value of $`n`$ is maximizing. However for a discrete set of $`ϵ`$ there will be more than one maximizing integer. Erdős and Nicolas \[6, p. 70\] later showed that for a given value of $`ϵ`$ there will be exactly one, two or four integers $`n`$ that attain the maximum value of $`\frac{\sigma (k)}{k^{1+ϵ}}`$. If the conjecture of Alaoglu and Erdős above has a positive answer, then no values of $`ϵ`$ will have four extremal integers. This would imply that the ratio of two consecutive colossally abundant numbers, the larger divided by the smaller, will always be a prime, and that every colossally abundant number has a factorization (2.5) for some nonempty open interval of “generic” values of $`ϵ`$.

Now we can explain the relevance of the Riemann hypothesis to the extremal size of $`\sigma (n).`$ Colossally abundant numbers are products of all the small primes raised to powers that are a smoothly decreasing function of their size. Fluctuations in the distribution of primes will be reflected in fluctuations in the growth rate of $`\frac{\sigma (n)}{n}`$ taken over the set of colossally abundant numbers. Recall that the Riemann hypothesis is equivalent to the assertion that for each $`ϵ>0`$,

$$|\pi (x)Li(x)|<x^{1/2+ϵ}$$

holds for all sufficiently large $`x`$. It is also known that if the Riemann hypothesis is false, then there will exist a specific positive constant $`\delta `$ such that the one-sided inequality

$$\pi (x)>Li(x)+x^{1/2+\delta }$$

is true for an infinite set of values with $`x\mathrm{}.`$ Heuristically, if we choose a value of $`x`$ where such an excess of primes over $`Li(x)`$ occurs, and then take a product of all primes up to this bound, choosing appropriate exponents, one may hope to construct numbers $`n`$ with $`\sigma (n)`$ exceeding $`e^\gamma n\mathrm{log}\mathrm{log}n`$ by a small amount. If the Riemann hypothesis holds, there is a smaller upper bound for excess of the number of primes above $`Li(x)`$, and from this one can deduce a slighly better upper bound for $`\sigma (n)`$. The analysis of Robin  gives a quantitative version of this, as formulated in Propositions 3.1 and 3.2 below.

One can prove unconditionally that the inequality (1.1) holds for nearly all integers. Even if the Riemann hypothesis is false, the set of exceptions to the inequality (1.1) will form a very sparse set. Furthermore, if there exists any counterexample to (1.1), the value of $`n`$ will be very large.

## 3 Proofs

The proof of Theorem 1.1 is based on the following two results of Robin .

###### Proposition 3.1

(Robin) If the Riemann hypothesis is true, then for each $`n5041`$,

$$\underset{d|n}{}de^\gamma n\mathrm{log}\mathrm{log}n,$$
(3.6)
where $`\gamma `$ is Euler’s constant.

#### Proof.

This is Theorem 1 of Robin .

###### Proposition 3.2

(Robin) If the Riemann hypothesis is false, then there exist constants $`0<\beta <\frac{1}{2}`$ and $`C>0`$ such that

$$\underset{d|n}{}de^\gamma n\mathrm{log}\mathrm{log}n+\frac{Cn\mathrm{log}\mathrm{log}n}{(\mathrm{log}n)^\beta }$$
(3.7)
holds for infinitely many $`n`$.

#### Proof.

This appears in Proposition 1 of Section 4 of Robin . The constant $`\beta `$ can be chosen to take any value $`1b<\beta <\frac{1}{2},`$ where $`b=\mathrm{}(\rho )`$ for some zero $`\rho `$ of $`\zeta (s)`$ with $`\mathrm{}(\rho )>\frac{1}{2}`$, and $`C>0`$ must be chosen sufficiently small, depending on $`\rho `$. The proof uses ideas from a result of Nicolas , \[13, Proposition 3\], which itself uses a method of Landau.

We prove two preliminary lemmas.

###### Lemma 3.1

For $`n3`$,

$$\mathrm{exp}(H_n)\mathrm{log}(H_n)e^\gamma n\mathrm{log}\mathrm{log}n.$$
(3.8)
#### Proof.

Letting $`t`$ denote the integer part of $`t`$ and $`\{t\}`$ the fractional part of $`t`$, we have

$$_1^{\mathrm{}}\frac{t}{t^2}𝑑t=\underset{1rtn}{}_1^{\mathrm{}}\frac{dt}{t^2}=\underset{r=1}{\overset{n}{}}(\frac{1}{r}\frac{1}{n})=H_n1.$$

Thus

$$H_n=1+_1^n\frac{t\{t\}}{t^2}𝑑t=\mathrm{log}n+1_1^n\frac{\{t\}}{t^2}𝑑t.$$
(3.9)
From this we obtain

$$H_n=\mathrm{log}n+\gamma +_n^{\mathrm{}}\frac{\{t\}}{t^2}𝑑t,$$
(3.10)
where we have set

$$\gamma :=1_1^{\mathrm{}}\frac{\{t\}}{t^2}𝑑t.$$

This is Euler’s constant $`\gamma =0.57721\mathrm{},`$ since letting $`n\mathrm{}`$ yields

$$\gamma =\underset{n\mathrm{}}{lim}(H_n\mathrm{log}n),$$

which is its usual definition. Now (3.10) gives

$$H_n>\mathrm{log}n+\gamma ,$$

which on exponentiating yields

$$\mathrm{exp}(H_n)e^\gamma n.$$
(3.11)
Finally $`H_n\mathrm{log}n`$ so $`\mathrm{log}(H_n)\mathrm{log}\mathrm{log}n>0`$ for $`n3`$. Combining this with (3.11) yields (3.8).

###### Lemma 3.2

For $`n20`$,

$$H_n+\mathrm{exp}(H_n)\mathrm{log}(H_n)e^\gamma n\mathrm{log}\mathrm{log}n+\frac{7n}{\mathrm{log}n}.$$
(3.12)
#### Proof.

The formula (3.9) implies, for $`n3`$,

$$\mathrm{log}H_n\mathrm{log}(\mathrm{log}n+1)\mathrm{log}\mathrm{log}n+\frac{1}{\mathrm{log}(n+1)}.$$
(3.13)
Next, (3.10) yields

$`\mathrm{exp}(H_n)`$ $`=`$ $`\mathrm{exp}\left(\mathrm{log}n+\gamma +{\displaystyle _n^{\mathrm{}}}{\displaystyle \frac{\{t\}}{t^2}}𝑑t\right)`$ (3.14)
$``$ $`e^\gamma n\mathrm{exp}\left({\displaystyle _n^{\mathrm{}}}{\displaystyle \frac{dt}{t^2}}\right)`$
$`=`$ $`e^\gamma n\mathrm{exp}\left({\displaystyle \frac{1}{n}}\right).`$
Since

$$\mathrm{exp}(x)1+2x\text{for}0x1,$$

we obtain from (3.14) that

$$\mathrm{exp}(H_n)e^\gamma n\left(1+\frac{2}{n}\right).$$
(3.15)
Combining this bound with (3.13) yields

$`\mathrm{exp}(H_n)\mathrm{log}(H_n)`$ $``$ $`e^\gamma n\mathrm{log}\mathrm{log}n+{\displaystyle \frac{e^\gamma n}{\mathrm{log}(n+1)}}+2e^\gamma (\mathrm{log}\mathrm{log}n+1)`$
$``$ $`e^rn\mathrm{log}\mathrm{log}n+{\displaystyle \frac{6n}{\mathrm{log}n}}`$
for $`n10`$. Using $`H_n\mathrm{log}n+1`$ and $`\frac{n}{\mathrm{log}n}\mathrm{log}n+1`$ for $`n20`$ yields (3.12).

### Proof of Theorem 1.1

$``$ Suppose the Riemann hypothesis is true. Then Proposition 3.1 and Lemma 3.1 together give, for $`n5041`$,

$$\underset{d|n}{}de^\gamma n\mathrm{log}\mathrm{log}n<H_n+\mathrm{exp}(H_n)\mathrm{log}H_n.$$

For $`1n5040`$ one verifies (1.1) directly by computer, the only case of equality being $`n=1`$.

$``$ Suppose (1.1) holds for all $`n`$. We argue by contradiction, and suppose that the Riemann hypothesis is false. Then Proposition 3.2 applies. However its lower bound, valid for infinitely many $`n`$, contradicts the upper bound of Lemma 3.2 for sufficiently large $`n`$. We conclude that the Riemann hypothesis must be true.

#### Acknowledgements.

I thank Michel Balazard for bringing Robin’s criterion to my attention, and J.-L. Nicolas for references, particularly concerning Ramanujan, Eric Rains for checking (1.1) for $`1n5040`$ on the computer, and Jim Reeds for corrections. I thank the reviewers for simplifications of the proofs of Lemmas 3.1 and 3.2.

AT&T Labs–Research, Florham Park, NJ 07932-0971, USA
email address: jcl@research.att.com