Spectral Properties of Exact Polarobreathers in Semiclassical Systems

: In this paper, we study the spectral properties of polarobreathers, that is, breathers carrying charge in a one-dimensional semiclassical model. We adapt recently developed numerical methods that preserve the charge probability at every step of time integration without using the Born–Oppenheimer approximation, which is the assumption that the electron is not at equilibrium with the atoms or ions. We develop an algorithm to obtain exact polarobreather solutions. The properties of polarobreathers, both stationary and moving ones, are deduced from the lattice and charge variable spectra in the frequency–momentum space. We consider an efﬁcient approach to produce approximate polarobreathers with long lifespans. Their spectrum allows for the determination of the initial conditions and the necessary parameters to obtain numerically exact polarobreathers. The spectra of exact polarobreathers become extremely simple and easy to interpret. We also solve the problem that the charge frequency is not an observable, but the frequency of the charge probability certainly is an observable.


Introduction
The subject of a charge moving when coupled to a lattice vibration has a long history [1,2].This phenomenon is very different from Ohmic conductivity in metals and semiconductors, where electrons in the conduction band or holes in the valence band move freely within a crystal.In insulators, where all the bands are either occupied or empty, but also in metals and semiconductors, an extra charge outside the bands may be attached to some atom or ion and have a low probability of moving to another position.The deformation it produces in the lattice can travel, and is known as a polaron [3][4][5].The charge-electron, or hole-changes its position to the same state in a neighboring atom.Polarons were first conceived as either a static solution or coupled to phonons, which were set in motion by an external electric field, but there are other possibilities as explained below.There is a probability that a charge may change to the neighboring site due to the overlapping of the electron wavefunctions in each atom, described by the transfer integral.It is clear that if the distance between atoms decreases, there will be an increase in the transition probability, and vice versa.Even for a symmetric vibration, the increase in the transition probability is much larger for shorter distances than the other way around, leading to an exponential dependence on the change in the potential barrier [6].
The level of description in which an electron or hole can be ascribed to a single atom is called the tight-binding approximation.It allows for a semiclassical treatment, in which the heavier atoms are described classically and the much lighter electron is described as a quantum particle.A system can be described in terms of the atomic states, and the quantum Hamiltonian is expressed in terms of the energies of these states and the transfer matrices; the latter often depend on the distance exponentially.
Large displacements of the atoms or ions in a crystal imply that the linear description is no longer valid.This nonlinearity may lead to the apparition of localized entities or, when traveling, solitary waves.Notable examples at the macroscopic scale are tsunamis, which travel thousands of kilometers at jet speed in a different manner to usual sea waves.Only slightly less impressive are tidal bores, solitary waves that travel upriver when the incoming tide encounters the outcoming one.Moreover, rogue waves in the sea appearing apparently from nowhere have been proven to exist by systematic satellite observations.These solitary waves can be kinks, solitons [7,8], or breathers, also called intrinsic localized modes [9][10][11].Solitons and kinks have a constant profile in the moving frame, where solitons become zero at ±∞, while kinks only at one of ±∞, and tend to a constant at the other infinity.Breathers also have a vibrating profile.Sometimes the distinction depends on the variables used.Important steps in breather theory were the mathematical proof of their existence [12,13] and Floquet and band theory for breathers [14][15][16], which later expanded to dark breathers [17,18] and multibreathers [19,20].They have been obtained with classical molecular dynamics [21,22] and ab initio molecular dynamics [23,24].Recently, the spectral theory of exact moving breathers was developed, showing that breathers admit a very simple description in their moving frame, where they have usually a single frequency [25,26].On the one hand, this makes possible the interpretation of numerical spectra, and on the other, it facilitates the integration methods.
The large displacements bringing neighboring atoms closer and increasing the probability of charge transfer between atoms in semiclassical systems lead to a strong coupling between charge and vibrations.The combined entities are called polarobreathers [27][28][29], or solectrons in some systems [30,31].The mathematical methods are often within the Born-Oppenheimer approximation [32], that is, supposing that the electrons are always in an equilibrium state with the configuration of atoms around them, due to their dynamics being much faster than the dynamics of the atoms.The authors have overcome this limitation by developing efficient numerical methods for semiclassical systems, based on splitting methods that conserve the charge probability at each step of time integration [33].On the experimental side of polarobreathers, it was recently discovered that by bombarding layered silicates and other materials with alpha particles, it was possible to measure electrical currents in the absence of an electric field; the electrons or holes were carried by nonlinear lattice excitations [34,35].
The purpose of this paper is to apply the theory of exact polarobreathers in their moving frame to a semiclassical system, using numerical integration methods developed in Ref. [33].In this way, we clarify the polarobreather spectra, facilitate their description, develop the mathematical methods for obtaining numerically exact polarobreathers, and compare their spectra with the approximate ones.The results can help to identify bands in the spectra of some materials where nonlinear vibrations couple to nonfree charges.The lattice part of the model without the charge is given by a Frenkel-Kontorova model [36,37] with the Lennard-Jones interaction potential, because it has been proven extremely good for producing long-lived breathers in two-dimensional systems [38,39].
The layout of the article is as follows.Section 2 presents the semiclassical system and its transformation into a Hamiltonian one.Section 3 deals with an important problem of the charge description: that the charge frequencies are not an observable of the system.In Section 4, a review of the theory of exact traveling excitations is presented.The linearization of the system is performed in Section 5. Very useful calculations, known as tail analysis, which advance many properties of nonlinear excitations, are presented in Sections 6 and 7 for stationary and moving polarobreathers, respectively.In Section 8, the developed numerical methods to solve dynamical equations and obtain exact polarobreathers are described.Sections 9 and 10 analyze stationary and moving polarobreathers, respectively, following a similar pattern: generation, description of the spectrum of approximate polarobreathers, obtainment of exact solutions and interpretation of their spectra, stability analysis, and path continuation.The article ends with the conclusions.

The Model
We consider a model of N particles, atoms or ions, with a classical and a quantum Hamiltonian imposing periodic boundary conditions.The classical Hamiltonian is given by In this equation, the variable u n is the separation of the particle n from the equilibrium position using the lattice unit as length unit; p n represents the momentum of the particle n, p n = mv n = m un equal to p n = un because the mass of the particle is chosen as the mass unit.Therefore, is the sum of the on-site energies of the classical system with respect to their equilibrium positions.It represents the interaction of the system with other systems, as other layers in silicates, or simply other atom chains in a crystal.Due to the periodicity of the crystal, the simplest form for the on-site energy for a particle is given by the first-order Fourier series, with the condition that u n = 0 is a stable equilibrium with zero energy: , known as the Frenkel-Kontorova model.The depth of the potential well is taken as unit of energy, i.e., U 0 = 1, and the period T 0 of small oscillations of a particle in the potential well U(u n ) is taken as time unit, i.e., ω 0 = 2π.The interaction energy ) is the sum of interaction energies between the particles of the system, that is, the energy of the interaction of the system with itself, described as pairs of Lennard-Jones potential between nearest neighbors.The Lennard-Jones potential has the generic physical property of growing to infinity when the distance between particles becomes zero and tending to a constant value, i.e., with zero derivative or zero force, when the particles separate.The Lennard-Jones potential is given by: The interaction potential has a minimum at r = σ with depth ε, being the ratio between the interaction energy and the on-site energy.We use ε = 0.05, a value that brings about the extraordinary mobility of breathers, both in the system without charge [38,39] and with it [33].We consider σ = 1, that is, also equal to the lattice constant, which is justified by the fact that there are no forces on the lattice particles at the equilibrium position.We also add the well depth ε to the interaction potential in order to have zero energy at the equilibrium distance.
The quantum Hamiltonian [40] for an electron or hole added to the system is given by: with the expectation value: n is the probability that the charge is located at site n with energy E n , and the sum of probabilities adds to one, i.e., The complex variable c n is the n-th component of |ψ = ∑ n c n |n , where |n represents the state for which the charge is completely localized at site n.The possibility of expanding the wave function in a basis of localized states |n , that is, states for which the electron is tight-bound to a single atom of the system, consists of the tight-binding approximation.
The evolution of the wave function |ψ is described by the Schrödinger equation ih∂ t |ψ = Ĥc |ψ .The Hamiltonian Ĥc in (2), acting on |ψ = ∑ m c m |m , leads to: where n|c m |m = c n δ n,m .The left hand side of the Schrödinger equation becomes ih∂ t |ψ = ∑ n ih ċn |n .Identifying the components on the basis of vectors {|n } in both sides of the evolution equation, we obtain the following evolution equation for the probability amplitudes c n 's: where τ is the reduced Planck constant h = h/2π in scaled units, with h being the Planck constant.In Refs.[25,41], which study a particular model for the movement of potassium ions in a silicate layer, the scaled units for length, energy, mass, and time are u L = 5.19 Å, the interatomic distance; u E = k e e 2 /u L 2.77 eV, the electric potential energy between two ions with unit charge e; u M = 39.1 amu, the mass of an ion; and u T = (u M u 2 L /u E ) 1/2 0.2 ps, the derived unit from the previous ones.In those units, τ = h/u E u T 0.00119.In this article, we do not refer to a specific model, but we use τ = 0.001 to have a correct order of magnitude.Note that the ratio of masses between a proton and an electron, approximately 2000, is coherent with the fastest movement of the electron with respect to the atoms in the lattice by a factor of 1/τ.Equation (3) is the expected value of the Hamiltonian Ĥc (2) when the system is in the state |ψ = ∑ n c n |n , i.e., H c = ψ| Ĥc |ψ .Ĥc is represented on the basis of |n by a matrix with elements n| Ĥc |m .E n = n| Ĥc |n are, therefore, the diagonal elements.The off-diagonal matrix elements, also called the transfer matrix elements, are J m,n = m| Ĥc |n , and are related with the probability of transition between the sites n and m.We suppose that they are zero, except for the nearest-neighbor particles, and that they depend on the distance between particles in the form: Therefore, the transfer elements and the probability of transition between sites increase rapidly when the particles become closer.The parameter values are α = 15, a value corresponding to the n = 3 orbitals; and I 0 = J 0 exp(−α) = 5 × 10 −4 , a good value for an insulator at equilibrium.Both values are also very convenient for producing traveling polarobreathers [33].E n represents the charge energy at the n-th site, and in general, it may depend on the variables of the lattice.We keep it in the general notation, but in the present paper, we consider it uniform and constant, which makes possible to set it to zero.However, the nonzero value is useful, as explained below in Section 3.
The expected Hamiltonian for the lattice and the charge is given by H = H L + H c , with Hamiltonian equations for the lattice variables: un = p n and ṗn = − ∂H ∂u n .
Therefore, the system of governing mathematical equations of the charge-lattice interactions reads as: where we have omitted the explicit dependence on t from the variables.Considering the time-dependent variables a n and b n , the real and imaginary parts of √ 2τc n , where the normalization is required, the governing Equations ( 8) and ( 9) can be written in the canonical Hamiltonian form with the Hamiltonian: which is the sum of the lattice and charge Hamiltonians (1) and (3), respectively, in variables u n , a n , p n , and b n .For the components u n , a n , p n , and b n of the canonical variables, the canonical Hamiltonian equations derived from (10) are: for all n = 1, . . ., N. Note that the total probability (4) is conserved along the solutions of Equations ( 11)-( 14) as the sum: is conserved.We solve the canonical Hamiltonian Equations ( 11)-( 14) with the exact charge probability (4) conserving, symplectic numerical method of Ref. [33] described in Section 8.1, to obtain the solution for the charge amplitude c n = (a n + ib n )/ √ 2τ and its probability |c n | 2 = (a 2 n + b 2 n )/2τ.In addition, in Section 8.2, we describe the numerical algorithm based on nonlinear least squares for the computation of numerically exact time-periodic stationary and moving polarobreathers.

Frequency Shift of the Charge Amplitudes {c n }
Let us define new probability amplitudes { cn } = {c n exp(−iµt)} with real µ that does not depend on n.This change conserves the density matrix [40], that is, the products cn c * m = c n c * m do not change.Consider a state vector |φ = ∑ n c n |n and a quantum observable, corresponding to a Hermitian operator Â.The expected value of the observable is given by A = φ| Â|φ = ∑ n,m A n,m c n c * m .It is obvious that c n and cn bring about the same expected values as well as derivatives of the expected values, as in (8).There is no physical form to detect the difference between the two sets of probability amplitudes; the kets |φ and | φ represent the same state.They are not, however, solutions of the same Equation ( 9), as can be seen by supposing that c n is a solution, and substituting c n = cn exp(iµt) in (9) we obtain: Multiplying (15) by exp(−iµt), we arrive at the equation: where E 0 = µτ.Therefore, multiplying the solution c n of (9) by exp(−iµt) brings about a new solution to Equation ( 9), but with the energy level shifted up by E 0 = µτ.This is a valuable property; exact solutions c n of (9) do not need to be periodic, but there may exist a solution cn = c n exp(−iµt), which is periodic and a solution of the same evolution Equation ( 9) with the energy shift E 0 = µτ to E n , and therefore, with a frequency shift µ = E 0 /τ for all the frequencies.As we have seen above, the magnitude E 0 does not appear in any physical observable.It is equivalent to fixing the gravitational potential energy at some specific height.
The actual solution could be obtained as c n = cn exp(iµt) from ( 16), but this is actually irrelevant, since we simply can obtain the actual energies and frequencies of the original system (9) by subtracting E 0 and E 0 /τ from the energy and frequency values, respectively.

Review of the Properties of Exact Moving Breathers and Solitons
Before analyzing the properties of polarobreathers, it is useful to review the properties of exact moving excitations-breathers or solitons.See Refs.[25,26] for details.We focus on the lattice variables u n , but similar analysis applies to |c n | 2 and c n .Figure 1 illustrates the properties below.
Exact traveling wave: An exact traveling wave with velocity V b is characterized by a function or sum of functions of the form: with f a 2π periodic on its second one, and with the following condition below.

Fundamental time and step:
There exist a fundamental time T F and an integer number called the step s, so the profile f repeats exactly after a time T F .The velocity is then Localized solutions: If f is localized in the first variable, it may be a breather, soliton, or kink.We will refer to it as a breather for simplicity.Then, we write the second variable as Ω b t.If it is delocalized, it is an exact extended traveling wave.
Fundamental frequency: ω F = 2π/T F .The relevant frequencies are integer multiples of ω F , among which is the frequency associated to the velocity, i.e., the frequency at which a traveler in the moving frame encounters particles 2πV b = sω F .
Frequency in the moving frame: Ω is the frequency in the moving frame, because if an observer moves with the breather, then n − V b t ≈ 0, and there remains a single frequency Ω, Ω b for the breather.We do not take into account the variation due to the discreteness.

Exactness:
The exactness condition implies that Ω = mω F for the integer m, or m b for the breather.

Harmonic modes:
The breather, as any function of (n, t), can be expressed as a sum of harmonic waves, also called modes u q = A q exp i[qn − ω q t] , with ω q being the laboratory frequency of each mode.Resonant modes: Resonant modes with the breather are all the modes that advance the step s in the time T F , that is, they are also exact with the same step and fundamental frequency.They can be written as Resonant lines: Therefore, the laboratory frequencies of the resonant modes are given by ω q = mω F + qV b .They form parallel straight lines called resonant lines, with slope V b , and cross the vertical axis (q = 0) at Ω m = mω F , for the m integer.All modes in a resonant line travel with speed V b and have the same frequency in the moving frame Breather line: All the modes in the spectrum of the localized exact solution u n are within one resonant line, with m = m b , that is, the breather line.Note that the breather line ends at q = π and reappears at q = π − 2π = −π as a different resonant line.

Breather frequency in the moving frame:
The breather line crosses the vertical axis (q = 0) at Ω b = m b ω F , the frequency of the breather in the moving frame.
Transformation to the moving frame: Changing the frame of reference to the one moving with the breather, it is equivalent to the transformation ω → ω MF = ω − qV b .The resonant lines and the breather line become horizontal lines with the frequencies in the moving frame.

Soliton or kink:
If Ω b = 0, the excitation is a kink or soliton, that is, a static profile in the moving frame.
Wing: If a resonant line crosses the dispersion band, there is a resonant phonon, i.e., a solution to the linearized equation, that may bring about a wing.It is an extended wave that travels with the localized solution, becoming a pterobreather, pterosoliton, or pterokink.The wing may be an integral part of the solution to the nonlinear dynamical equations, that is, it will not exist without the wing [42].

Commensurability condition:
The breather line frequency change for ∆q = 2π is ∆Ω b = V b ∆q = 2πV b , the velocity frequency, which is given by sω F .Then, is a rational number.This important property allows for the determination of m b , s, ω F , Ω b , and V b = ∆Ω b /2π, by simple inspection of the breather spectrum.∆Ω b is also the difference in frequency between the breather line and its continuation when π → −π.The moving frame frequency and the velocity frequency are commensurate.

Linear Approximations
In this section we obtain, the linearized equations corresponding to the dynamical Equations ( 8) and ( 9).We will use them to obtain the dispersion relations (DRs), also called the phonon bands, both for the variables u n and c n , while there is no dispersion relation for ρ n = |c n | 2 .They provide the vibrational modes at small amplitude, the phonons, and they are necessary to understand the spectrum of larger-amplitude polarobreathers formed with modes that separate from the linear ones, as observed in Sections 9 and 10.In general, the larger the amplitude, the larger the separation from the phonon band.
The dispersion relations are also useful for the tail analysis used in Sections 6 and 7, which supposes that a localized nonlinear solution has a core of large amplitude that decreases at sites further away from the core, forming a tail of decreasing amplitude.A few sites away from the core, when the amplitude is small enough, the lattice variables would abide by the linear dispersion relations, with a decreasing exponential solution valid only to one side of the core and at some distance from it.This simple method is extremely useful for predicting the properties of nonlinear excitations, such as the increase or decrease in frequency with respect to phonons with the same wavenumber.In this way, the properties of the polarobreathers in Sections 9 and 10 can be easily understood.
If we expand the terms in the dynamical Equations ( 8) and ( 9) to the first order in their variables, using and also seeing (6), we obtain: Keeping only the linear terms in ( 18) and ( 19), the lattice and the charge decouple, and we obtain a fully linearized system of equations: The solutions of linear equations have an exponential form with imaginary exponents if they are bounded, that is, they have the form ).The substitution of these expressions into the equations above leads to the dispersion relations: We have included these dispersion relations in Figure 1 and in all XTFT plots for reference.Note that the variables u n and c n become decoupled at the linear limit.Note also that wavenumbers q = 0 and q = ±π correspond to stationary solutions.For q = 0, the particles vibrate in phase with the same frequency, and for q = ±π, they vibrate with an alternate pattern.
There could be some doubt about which variables the linearization should be performed in (18), because only the density matrix elements ρ n,m = c n c * m appear.However, the linearization with respect to ρ n,m and u n results in the linear terms on ρ n,m becoming zero, and (8) does not change.There is no such problem for (19), as it is already linear in c n terms.
In Sections 6 and 7, using the system of linearized Equations ( 20) and ( 21), we perform the tail analysis of stationary and moving localized excitations, respectively.We will compare the results in these sections with the present one.

Tail Analysis of Stationary Localized Excitations
For the tail of a stationary localized excitation, we propose the ansatz: (24) which is valid for n > 0 (for n < 0, we change the sign in front of n in the first exponentials).Note that both ξ L,c = 0 implies no localization and extended waves.With the ansatz (24) and assumption of the constant charge energy E n = E 0 , for all n, from the Equations ( 20) and ( 21) we obtain the following equations for the lattice and charge frequencies: respectively, where q is the wavenumber.The momentum is given by p c = τq, with τ = τ/2π or h in physical units.However, we will use both terms for q when there is no confusion, as they represent the same physical observable in different units.
If the frequencies ω L,c are not real, that will imply decaying solutions; therefore, either both ξ L,c = 0, where we obtain extended solutions and the linear dispersion relations for the lattice and charge (22) and (23), or sin(q) = 0, i.e., q = 0 or ±π, and for q = 0 from (25) and (26), we obtain: From (27), it can be seen that the localization drives the frequencies below the linear spectrum (above for negative ω L ), and further away the larger the localization is.
For q = ±π, there is the opposite effect, i.e., localization increases the frequencies above the linear spectrum (below for negative ω L ).From ( 25) and ( 26): , Note that the localization parameters ξ L,c do not need to be the same in (25) and ( 26) and ( 27) and (28) for the lattice and the charge.A good estimate would be 2ξ c = ξ L , as in the lattice Equation ( 8), the charge amplitude c n products c n c * n±1 appear.Furthermore, has the same pattern as As the ±π modes are not actually traveling, the two modes will appear mixed, i.e., where A, B ∈ C. In this way, we can obtain the time dependence on the charge probability If we write AB * = |A||B| exp(iδ) for δ = arg(A) − arg(B), we obtain: In this way, we will observe in the charge probability spectrum two frequencies: ω ρ = 0, corresponding to a stationary deformation, and ω ρ = 4 I 0 τ cosh(ξ c ), where the charge probability spectrum is independent of the E 0 value.Note that as |c n | 2 is an observable, its frequencies can be measured; from them, we can deduce the difference in frequencies from the q = π and q = 0 modes of c n , but not their actual value, because there are no physical means of knowing which is the shift E 0 /τ of the c n frequencies.We choose for our convenience the value of E 0 that makes c n periodic with a commensurate period with the lattice one.However, that selection has no physical consequences, and we will plot the c n spectrum corresponding to E 0 = 0 for simplicity.

Tail Analysis of Moving Localized Excitations
We can repeat the tail analysis above for the moving solutions.In this case, the trial solutions in (24) are changed to: These ansätze are valid for n > V b t or for n < V b t changing ξ L,c ∈ R + to −ξ L,c .Considering first the ansatz (29) for u n , in the moving frame n = V b t, the frequency is that is, u n+1 is vibrating with the same frequency Ω L but with a smaller amplitude and a difference of phase q.The general solution for a traveling breather would be a sum of terms as u n in (29), with the same frequency Ω L in the moving frame.The parameter ξ L would be also dependent on each mode.The larger ξ L is, the more localized the specific localized mode would be, and likewise for c n with the ansatz (30).We also assume that the general solution is exact, that is, it repeats after some fundamental time T F displaced an integer number of sites s = V b T F , called the step.Then, all the modes (29) have the same properties, that is, V b , T F , and step s.The fundamental frequency is defined as ω F = 2π/T F .The exactness condition implies that Ω L = m L ω F for the integers m L and V b = s(ω F /2π).See Ref. [25] for more details.The fundamental period T F could, in principle, be different for the lattice and the charge amplitude, but for the common system to be exact, there should be integer multiples of both periods to obtain a common fundamental period T F and step s.
Substituting the ansatz (29) for u n into the linear Equation (20), we obtain: Compare (31) with Equation ( 25) of Section 6. Expanding the left hand side of (31), we arrive at two equations for the real and imaginary part: For compactness, we use the laboratory frequency: of the moving mode with wavenumber q.
For delocalized waves, with ξ L = 0, from (32), we recover the lattice dispersion relation (22) and the group velocity of the linear extended waves: for ξ c = 0, The latter equation indicates that there is only one mode with a given V b : the one where the slope of the dispersion relation is precisely V b .This is why it is not possible to have a coherent wave at the linear limit, as every mode has a different velocity.
As cosh(ξ L ) ≥ 1 and increases with the localization parameter ξ L , in addition, if π/2 < |q| < π, cos(q) < 0, and if 0 < |q| < π/2, cos(q) > 0, then from the first equation of (32), we find that Therefore, it can be seen that ω 2 c is above the linear frequency for the same |q| ∈ [π/2, π].However, when |q| becomes smaller than π/2, the two terms of the nonlinear correction have different signs, but they are also clearly above the linear frequency in a close proximity to |q| = π/2, at which point the negative term is zero.The main conclusion is that a solution such as the one proposed is possible below the linear spectrum, closer to q = 0, and above the linear spectrum, closer to |q| = π.The latter would be the case for our system.
Let us analyze the ansatz for c n , given by (30).The substitution into Equation ( 21) with the constant E n = E 0 leads to the following equation: where the charge amplitude laboratory frequency is given by: For the comparison with the stationary case, compare Equation ( 33) with (26) of the previous section.
Equation (33) leads to two equations for the real and imaginary part: The first equation gives the energy of the system, and the second gives the velocity.Note that a change in frequency ωc − E o /τ, corresponding to cn = c n exp(iE 0 /τ), is the frequency with the corresponding energy Ẽc = τ ωc − E 0 for the system with Ẽ0 = 0.
For ξ c = 0, we recover the linear dispersion relation (23) together with the group velocity: Again, there is only one linear mode for a given velocity V b , and coherent wave packages of c n are not possible at the linear limit or close to it.

Numerical Methods
In this section, we describe numerical methods used to solve the canonical Hamiltonian system ( 11)-( 14) and obtain numerically exact polarobreather solutions.

Numerical Integration of the Canonical Hamiltonian Equations
To solve the canonical Hamiltonian Equations ( 11)-( 14) numerically, we consider the exact charge probability (4) conserving, second-order semi-implicit splitting method PQDWDQP from the recently proposed splitting method class specifically developed for the semiclassical Hamiltonian dynamics of charge transfer in nonlinear lattices [33].The symmetric and symplecticitypreserving method PQDWDQP is constructed by splitting the total Hamiltonian (10) in the sum of the following Hamiltonians: H = H Q + H P + H D + H W , where Importantly, corresponding Hamiltonian systems associated to each Hamiltonian ( 34)-(36) can be solved exactly, and solutions are identified with the analytic symplectic flows φ Q t , φ P t , and φ D t , respectively.Meanwhile, we solve the Hamiltonian system of (37) with the symplectic implicit midpoint rule, whose solution we identify with the flow map ψ W h , where h > 0 is the time step.Thus, the numerical method PQDWDQP with the flow map ψ h is a symmetric composition of flows φ Q t , φ P t , and φ D t as well as the flow map ψ W h , i.e., where the symmetry of the method follows from the construction, symplecticity follows from the composition of symplectic maps, and the charge probability conservation (4) follows from the implicit midpoint step ψ W h , which preserves quadratic invariants [43].Advancing from state (u n , a n , p n , b n ) to a new state (U n , A n , P n , B n ) after one time step h for all n = 1, . . ., N, the solution given by the flow φ Q h of the Hamiltonian system associated with the split dynamics (34) reads: Similarly, the solution given by the flow φ P h of the Hamiltonian system associated with the split dynamics (35) reads: Notice that during the applications of the flow maps φ Q h and φ P h , only the lattice variables u n and p n get updated, while the charge variables a n and b n remain unchanged.
To find the exact expression of the flow φ D h of the Hamiltonian system associated with the Hamiltonian (36), we state the associated system of differential equations: where n = 1, . . ., N. From ( 38)- (41), it is easy to see that we obtain decoupled harmonic oscillator Equations ( 39) and ( 41) for the variables a n and b n , which can be solved exactly for all n and any initial condition (a 0 n , b 0 n ), i.e., where ω n = τ −1 E n .Thus, the exact (explicit) solution given by the flow φ D h of the Hamiltonian system (38)-( 41) is in the following form: If all E n = 0, then the flow map φ D h is just an identity map.We conclude this section by listing the explicit representation in the component form of the implicit midpoint map φ W h : taking into account the periodic boundary conditions.Then, where x(T F ) is the numerical solution of ( 11)-( 14) at time T F obtained with exact charge probability (4) conserving the splitting method PQDWDQP and time step h of Section 8.1.
We recall that T F is the fundamental period and the step s = 0 for the stationary polarobreathers, while s > 0 if the traveling polarobreather moves to the right and s < 0 if the traveling polarobreather moves to the left.Note that the objective function in (46) is minimized for x, Lagrange multipliers λ, and the constant charge shift energy value E 0 (parameter value in the system) with the fixed value of T F and s, which are estimated from the approximate solution spectra; see Sections 9 and 10.
The vector function g(x) ∈ R M contains imposed constraints such as the charge probability (4) conservation, without the loss of generality b N = 0, to eliminate the charge rotational invariance by an arbitrary angle θ, i.e., canonical Hamiltonian Equations ( 11)-( 14) are invariant under the transformation cn = c n exp(iθ).Thus, For the computation of exact stationary polarobreathers, we also impose, in addition to (48), that all p n = 0, which leads to time-periodic stationary polarobreather solutions p n (T F ) = p n (0) = 0.The quadratic constraint in (48) demonstrates the necessity for the exact charge probability (4), conserving numerical methods, e.g., the use of PQDWDQP.
We consider the damped Gauss-Newton algorithm to solve the optimization problem (46), i.e., ( 46) is reduced to the regularized linear least squares: where µ ≥ 0 is the regularization parameter adjusted with each iteration, and k is the iteration index, where k = 0, 1, 2, . . ., x k , and E k 0 are the solution and energy E 0 value at the k-th iteration, respectively.J(x k , E k 0 ) ∈ R 4N×4N+1 is the Jacobian matrix of (47), while G(x k ) ∈ R M×4N is the Jacobian matrix of the constraint vector function g(x).The Jacobian matrix G(x k ) can easily be evaluated analytically.However, in theory, the Jacobian matrix J(x k , E k 0 ) can also be evaluated analytically but with high difficulty.Thus, we evaluate it numerically.Accordingly, increments (descent directions) ∆x ∈ R 4N and ∆E 0 ∈ R.
The minimum of the linear least squares (49) is the solution to the following linear system of equations: where I and 0 are identity and zero matrices of appropriate dimensions, and We set the stopping criteria on the maximal absolute errors of functions f (x; E 0 ) and g(x) with tolerance 10 −14 .Thus, the converged solution satisfies (pointwise) time periodicity T F and the step s conditions (47), as well as the system constraints g(x), to high numerical accuracy, which are then preserved by the numerical method PQDWDQP to perform spectral analysis.
For all µ > 0, x k , and E k 0 , the matrix ) has a full row rank and the linear system above has a unique solution such that ∑ N n=1 (a k+1 n ) 2 + (b k+1 n ) 2 > 0 and b k+1 N = 0, which implies that the Jacobian matrix G(x k+1 ) has a full row rank and the linear system has a unique solution with x k+1 and E k+1 0 .To verify this, assume the opposite: = 0 for all n, i.e., increments ∆a n = −a k n and ∆b n = −b k n .From the unique solvability of the linear equations with x k and E k 0 , from the linear equations, we obtain: which yields a contradiction, since τ > 0. Thus, the linear system has a unique solution for all x k and E k 0 if x 0 satisfies all constraints.To update the regularization parameter µ after each iteration step in (49), we compute the gain ratio value: The gain ratio value υ allows us to adjust the value of µ such that it decreases as we approach the minimum of the nonlinear least squares (46).We were able to obtain satisfactory convergence results; the numerical results demonstrated superlinear convergence if good starting x 0 and E 0 0 values are provided, with the update strategy for µ, that is, double the value of µ, if υ < 0.25, remaining the same value for the next iteration, unless υ > 0.75, then a three-times smaller value of µ is considered.As an initial value, we chose µ = 10 −6 max diag J(x 0 , E 0 0 ) T J(x 0 , E 0 0 ) .

Stationary Polarobreathers 9.1. Generation of Approximate Stationary Polarobreathers
The methods used to produce polarobreathers in this model (solving ( 11)-( 14)) are extremely efficient; they have been used in previous works in one and two dimensions [33,38,39].In this section, we present the method for the stationary case; the moving polarobreather case will be presented in the following Section 10.
We introduce nonzero initial conditions only for the velocities or momenta: The parameter γ is related with the kinetic energy delivered to the system as K E = 5γ 2 .The reference index n * can be arbitrary because the system is periodic, but we usually take n * = N/2, for visual plotting purposes.The charge wave function is located initially with probability one, with the pattern: That is, the charge is completely localized at n * with probability one.This combined pattern (50) and (51) proves to be very efficient in obtaining quite good stationary solutions with long life.Other methods and patterns have also been investigated, for example, producing a local compression, which often brings about a stationary breather and two traveling ones in opposite directions.Moreover, different patterns for the positions and momenta have been considered.Many of them work quite well, but this is the preferred one for simplicity and for obtaining results.We describe the results for γ = 0.4, corresponding to the kinetic energy K E = 0.8, in scaled units.Other values of γ bring about qualitatively similar results.
We perform the 2D discrete Fourier transform in positions and time (XTFT) of timeseries data on the variables u n , |c n | 2 , and c n , which are represented in Figure 2-top.For reference, the dispersion relations ( 22) and ( 23) for u n and c n are also plotted with gray solid and dashed lines, respectively.Numerical results were obtained, setting all E n = 0 in ( 11)- (14).
As the first two quantities are real, for the XTFT components, it holds that F(−q, −ω) = F * (q, ω) and |F(−q, −ω)| = |F(q, ω)|, which is a symmetry that can be observed in the corresponding two upper plots.Usually, in this case, the negative frequencies are not represented, but here, they are included for comparison with the XTFT of c n , where the symmetry does not hold, as c n is complex.Some main features from the top plots of Figure 2 can be discerned: 1.
For the XTFT of u n , there appear two horizontal lines at ω = 0, centered around q = 0, and some frequency ω L /2π 1.33 above the dispersion relation and centered around q = ±π.This means that the u n breather is composed of a soliton, that is, a static deformation with a displacement largely in phase and a staggered vibration above the u n phonon spectrum, i.e., a nonlinear vibration, as demonstrated in Section 6.The static solution appears due to the asymmetry of the Lennard-Jones potential well, which makes compression harder than expansion, and therefore, oscillations with respect to the equilibrium distance are larger for expansion than for compression.

2.
For ρ n = |c n | 2 , two horizontal lines appear, one at zero frequency, a stationary soliton close to q = 0, that is, with nearest neighbors in phase; and also at frequency ω ρ /2π 0.75, close to the modes q = ±π, that is, with a staggered profile.The soliton here is necessary, as ρ n is a positive quantity, so the vibration has an alternating pattern around a stationary one.This means that there is a small change in probability between neighboring particles with the frequency ω ρ .This was also explained in Section 6.
Depending on the nonlinearity and the system, the interchange of probability will be larger or smaller.

3.
For c n , we find three main frequencies, two of them ±ω c ±0.375, close to the c n phonon spectrum, one above and the other below.The upper one is around q = ±π, i.e., with a staggered profile; and the lower one is around q = 0, that is, with a bell profile.These two frequencies are explained in Section 6.The other two are located at ±(ω L − ω c ).These appear because the quantum Hamiltonian has the time periodicity of u n , which appears in the transfer matrix elements.

4.
As deduced in Section 6 the ρ n = |c n | 2 frequency is equal to the difference between the phonon frequencies ±ω c , in this case being ω ρ = 2ω c .For E n = E 0 = 0, the two frequencies would be E 0 /τ ± ω c .See Section 3 for details.

5.
There are some other lines of weaker intensity in XTFT of |c n | 2 , but especially for c n .We can observe some phonons for u n and even more for c n , where they occupy the whole c n -phonon band, but not for |c n | 2 , as there is no dispersion relation for |c n | 2 , because its evolution depends on the other terms of the density matrix Results on the density matrix for polarobreathers will be discussed and published elsewhere.

Exact Stationary Polarobreathers
The approximate solution described above is good enough to be used as a seed in the numerical algorithm of Section 8.2 for obtaining numerically exact polarobreathers.The two main frequencies are ω L /2π 1.33 and ω ρ /2π 0.75, with corresponding periods: T L 1/1.33 and T ρ 1/0.75.The frequencies of c n are analyzed in a different way, which is explained below.We can observe that the periods and frequencies are approximately commensurate: 9T ρ 12 16T L .Therefore, T F = 12 is a good estimate for the common period for both u n and |c n | 2 , and it is the fundamental time or period taken for the whole system [25].
We suppose that E n = E 0 in Equations ( 12) and ( 14), with initial value , where ω ± c are the upper and lower frequencies above and below the c n -phonon band; see the top right plot of Figure 2. We recall that the value of E 0 also gets adjusted and found in the damped Gauss-Newton method of Section 8.2.This is an essential degree of freedom to obtain exact periodic solutions with the desired numerical accuracy ∼ 10 −14 .
The profile of the exact polarobreathers and their (lack of) change with time for u n and ρ n = |c n | 2 can be seen in Figure 3, where the solution is visualized in time after five fundamental periods.The small asymmetry of ρ n in Figure 3b can be observed.The XTFT of the exact solution is plotted in Figure 2-bottom; the main features are as follows: 1.
All the phonons have disappeared from the dispersion bands.

2.
Extra bands have also disappeared, except for the ones described above, which have become much more defined.

3.
In the XTFT of u n , the zero-frequency component corresponding to the stationary soliton and the frequency ω L slightly above the positive phonon band, centered at ±π (also the symmetric band at −ω L ), have remained.These features are in accordance with the theory described in Section 6.

4.
In the XTFT of c n , already shifted to E 0 = 0 by the numerically found E 0 value, only two bands-±ω c slightly above the c n phonon band at q = ±π and below at q = 0, which correspond to a hard potential case-have remained; see Section 6. 5.
Furthermore, for the XTFT of c n , there is a weak negative frequency −(ω L − ω c ) corresponding to the forcing by the matrix transfer elements J(u n+1 − u n ), which change with ±ω L , the u n frequency.The corresponding band is symmetric in q, but does not include q = 0.This means that it is a stationary wave: the sum of waves traveling in opposite directions with wavenumbers around q ±π/3 or wavelength λ 6.With greater initial kinetic energy, the positive frequencies ω L also appear above the positive phonon band and are centered at q = ±π.The periodicity and asymmetry of the exact solution can be appreciated.
We conclude that although ω ± c are not observable, their difference ω ρ is indeed observable, and appears in the spectrum of the charge probability.Thus, ω ρ may appear in the spectrum of physical systems, providing a valuable insight into the states of extra electrons or holes of the system.

Stability of Exact Stationary Solutions: the Switching Mode
We can numerically obtain the Floquet matrix, that is, the ∂S s x(T F )/∂x(0) at the exact solution.As the system is symplectic, if there is an eigenvalue λ, then 1/λ is also an eigenvalue.As the system is also real, if λ is an eigenvalue, then so are λ * and 1/λ * as the conjugate of 1/λ.Therefore, the eigenvalues come in quadruplets if they are complex with |λ| = 1, and in pairs λ, λ * if |λ| = 1, or λ, 1/λ if λ is real.The perturbation of the system with an eigenvector results in the perturbation growing as λ r with r periods.Therefore, the (linearized) system is only stable if |λ| = 1 for all eigenvalues.When changing a parameter as the frequency, the complex eigenvalues have to leave the unit circle as a quadruplet; therefore, for nonreal λ, two pairs of complex eigenvalues in the unit circle have to first collide for an instability to appear.This will be a Hopf bifurcation to a set of solutions with different periods (arg(λ) ± 1)T F /2π.Most of these eigenvalues correspond to phonons outside the core of the polarobreather.However, two eigenvalues colliding at (−1, 0) can get out of the unit circle in a period-doubling bifurcation.Two eigenvalues can collide at (1, 0) and get out of the circle as two real eigenvalues: one larger than one and another smaller, corresponding to two eigenvectors: one growing and other contracting to conserve the area in the phase space.
The system structurally has four eigenvalues at (1, 0), corresponding to two growth modes, i.e., small changes in amplitude in u n or c n ; and two phase modes, corresponding to two small changes in phase or time origin.Their meaning is that solutions also exist with almost the same period and slightly different amplitudes or phases.
There is a peculiarity in our Floquet matrix calculation: it is performed at T F = 16T L = 9T ρ = 12; therefore, an eigenvalue λ for the u n variables with period T L appears as λ 16 , and for c n , it appears as λ 9 compared to the period for |c n | 2 .This means that instability eigenvalues appear much more grown or contracted than usual, and the numerical imprecision for the eigenvalues at λ = 1 is amplified.
We can observe the eigenvalues of the exact polarobreathers in Figure 4a.The reference circle appears populated with the phonon eigenvalues, and the four structural eigenvalues exist at +1.However, there is also a pair of real instability eigenvalues, which are very unstable, as just explained.For the breather period, they would be reduced to around 1.1 and 0.9.The corresponding two eigenvectors appear in Figure 4b for ρ = |c n | 2 , together with the solution at t = 0. We observe a small asymmetry in the profile of ρ n , and the eigenvectors tend to increase the probability at the neighboring particles away from the center of localization.Long-time simulations confirm this interpretation.After some time, the polarobreather switches a position; after some time, it switches back, and so on.Therefore, in spite of the instability, which gets stabilized by the nonlinearity of the lattice, the breather and charge probability are stationary.They do not disperse or travel, but experience quasiperiodic switches between two neighboring sites.Figure 5 shows the switching behavior of a polarobreather in a long-time simulation.

Moving Polarobreathers
In this section, we follow the scheme of the previous section on stationary polarobreathers.In spite of the additional complications by the mobility, this will be easy to explain based on that section.

Generation of Approximate Moving Polarobreathers
The method used to produce moving polarobreathers in this model is also extremely efficient.We introduce as initial conditions the following initial pattern for the momenta or velocities: while the particle displacements are set at zero.The parameter γ is a measure of the modulus of the initial kinetic energy K E = 3γ 2 .The charge amplitude function with a probability of one is initially located with the pattern (51).Figure 6-top shows the XTFT plots obtained with γ = 0.6, K E = 1.08.They look very similar to Figure 2-top.As already seen in Section 4, the transformation ω L → ω M = ω L − qV b corresponds to the description of the system in the comoving frame, where all the phonons and excitations traveling at the same velocity V b appear as horizontal lines, corresponding to a stationary solution with a single frequency.Quite a few bands appear.For u n and |c n | 2 , real numbers, we only comment on the positive part of the spectrum due to its symmetry.One necessary objective to obtain moving exact polarobreathers is to identify the common fundamental time T F and step s.

1.
We observe the XTFT of u n .There are three localized waves traveling at the same speed, and some phonons.The three localized waves are a soliton; a breather; and with weaker intensity, what we could call a quasilinear breather, very close to the phonon band.We know that the soliton and main breather are characteristic components of this system, with strong asymmetry in the coupling potential.Therefore, we discard the weaker breather as an of the breather not being exact.In Section 4 it is deduced that the breather frequency in the moving frame Ω L is the frequency of the breather line for q = 0.The commensurability relation (17) was also obtained: with m L and the step s L being the integers we have to find.We measure Ω L /2π 1.2257 and V b 0.3065; then Ω L /2πV b = 3.999 4. Therefore, the simplest values are 0.3065, T F,L 3.2626.

2.
For the XTFT of |c n | 2 , also for positive frequencies, we observe four resonant lines with intensity.One is the soliton, characteristic of the XTFT of a positive quantity.The many lines are a sign of high nonlinearity, which corresponds to a very localized |c n | 2 , and is therefore close to one for some n.The linear approximation is not at all valid, and there are many harmonic frequencies.

3.
For the XTFT of c n , we have to observe the frequency differences, as the solution is invariant to a global shift in frequency, as explained in Section 3. We observe many phonons in the dispersion band, two positive bands centered around q = ±π, and two negative bands centered at q = 0. Two bands are very close to the phonon band, indicating a strong interaction with the linear modes.Observing the two truly nonlinear bands, their difference in frequency is Ω c /2π 0.715, corresponding to Ω c = m c ω F,c .Using the commensurability relation (17) as above and with the same V b 0.3064, we obtain: We conclude that 0.1021, T F,c 9.7911.
Thus, we can set the common step s = 3 and estimate the fundamental time T F = 9.79, i.e., T F = 9.79 3T F,L T F,c , as well as finding an approximate value of E 0 0.0028 for the computation of the exact moving polarobreathers of the following section.

Exact Moving Polarobreathers
Using the numerical methods explained in Section 8 and the found step s and fundamental period T F , together with the estimated E 0 value, we can obtain numerically exact moving polarobreathers.Their XTFT is shown in Figure 6-bottom.There is a large difference between the approximate moving polarobreather and the exact one.

1.
For the XTFT of u n , phonons have been eliminated, and only a well localized soliton breather remains.

2.
For the XTFT of |c n | 2 , only a soliton remains, meaning that in the moving frame, it is reduced to a static deformation, corresponding to a charge probability traveling without vibration.The other part of the solution is a uniform probability spread through the lattice, indicating a small probability that the charge could appear at any site n.

3.
For the XTFT of c n , two intensity bands remain: one around q = 0, and a stronger one centered at q = ±π.The difference in frequencies is exactly the breather frequency, indicating that there is no vibration of c n , except the one driven by the change in the transfer elements J n,n+1 , which depend on u n .
To summarize, there is a common fundamental period T F , during which the polarobreather repeats itself with a shift of s = 3 sites.In T F , the breather at u n repeats three times, and in the moving frame performs m L = 4 oscillations for each repetition, totaling 12 repetitions.In the same duration of time T F , the c n variable oscillates m c = 7 times, while advancing s c = 3 sites.There is no vibration in |c n | 2 ; there is just the translation with velocity V b of the static profile with about half the probability located at a single site, and the rest is uniformly shared by the rest of lattice particles.The periodicity of the exact polarobreather solution can be appreciated in Figure 7, where the particle displacements u n (Figure 7a) and the charge probability |c n | 2 (Figure 7b) are illustrated in time separated by five fundamental periods T F = 9.79.

Stability of Moving Polarobreathers
Most of the properties for the eigenvalues of the Floquet matrix for the moving exact polarobreather are the same as in the stationary case (see Section 9.3), as we are observing the system in the moving frame where it is stationary.Eigenvalues of the Floquet matrix are illustrated in Figure 8.There is also a difference; now, we should expect six structural eigenvalues at the eigenvalue (1, 0).They correspond to the growth and phase mode for u n and c n , and two new eigenvalues associated with a small change of the velocity.The latter imply that there exists another solution with a small change in the velocity V b .This solution will not be an exact solution that will appear with the change in T F for constant s.The stability of the exact moving polarobreather can be observed in a very long-time simulation, shown in Figure 9.

Path Continuation
By changing the parameters T F and keeping the step and the ratio between the u n and c n periods, it is possible to obtain exact solutions with similar characteristics with a different fundamental frequency ω F , both for increasing and decreasing frequencies.For the decreasing frequencies, the convergence stops as the u n frequency approaches the phonon band.For increasing frequencies, the path continuation also eventually stops to converge.These results will be published elsewhere.

Conclusions
We considered a classical model without charge that shows long-lived traveling breathers, and used it to construct a semiclassical model adding the coupling to a charged quantum particle-electron or hole.We analyzed the system and advanced its properties through linear and tail analysis.We adapted recent numerical methods developed by the authors to produce efficient symplectic algorithms that preserve symplecticity and charge conservation at every time step.They allow for efficient integration of the dynamical equations without the Born-Oppenheimer approximation, i.e., both the lattice and the charge are out of equilibrium.We constructed initial conditions for the lattice variables and the charge amplitude that bring about long-lived polarobreathers.The analysis of their spectra allows important parameters to be obtained not only for exact stationary solutions, but in particular for traveling ones, such as the step, fundamental frequency, and frequency in the moving frame, both for the lattice and charge variables.We developed a method for obtaining numerically exact polarobreathers, using as an initial seed the approximate solutions and the parameters found in their spectra.An important aspect of this method is dealing with the system invariance under a change in the frequency of the charge amplitude.The spectrum of the exact polarobreathers is much more simplified, leaving only a breather and a soliton for the lattice variables, a soliton for the charge probability, and two breatherlike exact solutions for the charge amplitude, with a frequency difference that matches the breather one and corresponds to a soliton and a breather with a frequency shift.However, for approximate polarobreathers that might appear in real systems, the charge probability shows small oscillations with frequencies that can be related to the frequencies of the charge amplitude.The obtained results may allow for the identification of charge coupling and its properties in the spectrum of real systems.

Figure 1 .
Figure 1.Plot of the XTFT of exact breathers together with the dispersion relations for u n (continuous line) and c n (dashed line), the latter only for reference.(Left) Stationary soliton-breather.(Right) Moving soliton-breather, resonant lines and the breather line, m b = 4.The breather frequency in the moving frame Ω b = m b ω F and the fundamental frequency ω F are also marked.In this case the step s = 1.See Section 4 for details.

Figure 2 .
Figure 2. XTFT of the lattice variables u n , charge probability |c n | 2 , and probability amplitude c n .(Top) Inexact stationary breather generated with γ = 0.4 and kinetic energy K E = 0.8.(Bottom) Exact stationary breather obtained using the approximate polarobreather as an initial seed.See text.

6 .
In the XTFT plot of charge probability ρ n = |c n | 2 , only the bands with frequencies at zero and ω ρ = ω + c − ω − c have remained, due to the election of E 0 = 0, ω ρ = 2ω + c (as well as −ω ρ , due to the symmetry).

Figure 3 .
Figure 3. Exact stationary polarobreather solution illustrated in time separated by 5 fundamental periods T F = 12.(a) Profile of the lattice variables u n .(b) Profile of the charge probability ρ n = |c n | 2 .The periodicity and asymmetry of the exact solution can be appreciated.

Figure 4 .
Figure 4. (a) Floquet eigenvalues of an exact stationary polarobreather.The unstable eigenvalues are so large and small because the common period T F is 16 times the breather period and 9 times the probability period.(b) Profile of the charge probability and the stable and unstable eigenvectors of the charge probability, where the latter corresponds to the switching mode.See text.

Figure 5 .
Figure 5. Long-time simulation of an exact stationary polarobreather.The alternating switching of position and charge can be observed.

Figure 6 .
Figure 6.(Top) XTFT of u n , |c n | 2 and c n for an approximate polarobreather generated with γ = 0.6 and kinetic energy K E = 1.08.(Bottom) XTFT of u n , |c n | 2 and c n for the exact moving polarobreather obtained using the approximate polarobreather as an initial seed.See text.

Figure 7 .
Figure 7. Profiles of an exact moving polarobreather in time separated by 5 fundamental periods T F = 9.79.(a) u n , with high localization and staggered profile.(b) |c n | 2 , corresponding to a charge probability localized at five sites and peaking at one site.

Figure 8 .
Figure 8. Floquet eigenvalues for the moving exact polarobreather.The six structural eigenvalues at (1, 0) can be seen, with some imprecision due to the numerical error and the long time used for the common T F for all the variables.See text.

Figure 9 .
Figure 9. Long-time simulation of an exact moving polarobreather showing its stability.