{'fraction_non_alphanumeric': 0.03392946682266421, 'fraction_numerical': 0.014708340297509611, 'mean_word_length': 4.640584492104643, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Propagation of extremely short electromagnetic pulses in a doubly-resonant medium

Propagation of extremely short electromagnetic pulses in a homogeneous doubly-resonant medium is considered in the framework of the total Maxwell-Duffing-Lorentz model, where the Duffing oscillators (anharmonic oscillators with cubic nonlinearities) represent the dielectric response of the medium, and the Lorentz harmonic oscillators represent the magnetic response. The wave propagation is governed by the one-dimensional Maxwell equations. It is shown that the model possesses a one-parameter family of traveling-wave solutions with the structure of single or multiple humps. Solutions are parametrized by the velocity of propagation. The spectrum of possible velocities is shown to be continuous on a small interval at the lower end of the spectrum; elsewhere the velocities form a discrete set. A correlation between the number of humps and the velocity is established. The traveling-wave solutions are found to be stable with respect to weak perturbations. Numerical simulations demonstrate that the traveling-wave pulses collide in a nearly elastic fashion.


Introduction
The recent demonstration of artificial materials (metamaterials) with the left oriented triplet of electric E, magnetic H and wave vector k of electromagnetic field [1,2,3,4] stimulated studies of nonlinear optical phenomena in such materials [5,6,7,8,9,10,11]. Nonlinear dynamics of extremely short optical pulses in left-handed materials was the subject of particular interest in several recent papers [7,12,13,14]. The first experimental realization of the left-handed property based on the resonant response of the artificial material to both electric and magnetic fields was described in [1]. To mention just one of the latest experimental achievements, Valentine et al [2] were able to observe the negative refractive index in the balk material in the optical range. A theoretical description of the electromagnetic wave interaction with such double resonance materials (DRM) was considered in [15,16,17,18,19]. Presence of two frequency intervals with different orientation of ( E, H, k) triplets is a characteristic feature of such materials.
Most of the studies of electromagnetic pulse propagation in DRM has been conducted in the slowly varying envelope approximation. On the other hand, there is a broad area of nonlinear optical phenomena taking place in the limit of extremely short pulses, when the slowly varying envelope approximation is not valid [20]. The case of extremely short electromagnetic pulses offers a new type of nonlinear interaction, when different frequency components of electromagnetic pulses have different orientations of the ( E, H, k) triplets.
The design of currently available DRM is based upon the use of embedded metallic structures whose size is on the same order as the spatial size of an extremely short electromagnetic pulse. Therefore a theoretical and numerical investigation of the currently existing DRM would require 3D computer simulation on Maxwell's equations that takes into account the strong inhomogeneity of composite materials. Recently, there have been introduced some qualitatively different approaches to design of DRM, including the use of multilevel atoms [21,24,25]; the latter gives rise to a spatially homogeneous medium. Possibilities of experimental realizations of such an approach were recently discussed in [22,23]. As a first step in the theoretical investigation of electrodynamics of homogeneous DRM in this paper we study a simple model of a homogeneous doubly-resonant medium. Even under such simplification, dynamics of extremely short pulses turn out to be quite complex.

Basic equations
The system of equations that describe interaction of coherent light with a medium consisting of molecules (considered as harmonic oscillators) is known as the Maxwell-Lorentz model [27]. In this work we use a version of the Maxwell-Lorentz system that is extended to account for simultaneous magnetic and electric resonances, with the magnetic susceptibility being linear, while the electric polarization being nonlinear. Consider the general form Maxwell's equations: For simplicity, we consider transverse electromagnetic plane waves propagating along the z-axis with the electric field E = (E(z, t), 0, 0) and the magnetic field B = (0, B(z, t), 0). Then the Maxwell equations transform to the scalar form: The system (2.4) must be closed by two additional equations describing the interaction of the electric and magnetic fields with the DR medium. As usual, it is convenient to avoid the 4π-factors by changing the units for M and P :M = 4πM,P = 4πP. In the sequel we drop the tildes over M and P.
Assume that the medium polarization is defined by the plasma oscillation of electron density, P tt = ω 2 p E Here ω p is an effective parameter characterizing polarizability of the medium; in the case of metallic nanostructures it would be the effective plasma frequency. To account for the dimensional quantization due to the confinement of the plasma in nanostructures one should include the additional term ω 2 D P , where ω D is the frequency of dimensional quantization. We take into account nonlinearity in the lowest order of P , which is P 3 . A more accurate analysis, based on a quantum mechanical approach [29] and experimental measurements [30] confirms validity of this assumption. Therefore we consider the modeling equation for the medium polarization dynamics in the following form where κ is a constant of anharmonisity. To account for magnetic resonances we use the standard model [15] M tt + ω 2 Here β is a constant characterizing magnetization. We represent equations (2.4), (2.5) and (2.6) in a dimensionless form by introducing τ = t/τ 0 (τ 0 = 1/ω p is the characteristic time), η = z/z 0 (z 0 = cτ 0 is the characteristic distance), q = P/P 0 (P 0 = ω p / √ κ is the maximal achievable medium polarization). It is convenient to normalize remaining variables as follows: m = M/P 0 , e = E/P 0 , h = H/P 0 . The system of dimensionless equations then takes the following form: where γ = κ/ |κ| ω 2 p , ω 1 = ω D /ω p , ω 2 = ω T /ω p . The system possesses the following conserved quantity: which is positive-definite for β < 1. For the traveling-wave solutions the conservation relation (2.8) yields conservation of electromagnetic energy 1 2 e 2 + h 2 dη = const (see [28] for details). A natural question arises is whether the system in (2.7) possesses any solitary-wave solutions. We address this issue in the following section.

Solitary wave solutions
Consider a traveling wave solution of (2.7), i.e., a solution that is a function of the variable ζ = τ − η/V. Then the PDEs in (2.7) become ODEs, and one obtains the following system: Upon the integration of equations (3.9) and (3.10) once, we get the algebraic conservation relations We are interested in a traveling-wave solution on the zero background, hence h = m = q = e = 0 at ±∞; therefore the constants of integration R = S = 0. This yields the following expressions for h and e h = a 1 m + a 2 q (3.13) e = a 2 m + a 1 q (3.14) where We insert expressions (3.13) and (3.14) for h and e into the equations (3.11) and (3.12) for q and m and obtain the following system of second order equations: This system can be diagonalized with respect to the second derivatives by the means of the transformation Instead of the second order system (3.16) we will consider the following equivalent 4 × 4 first order system Obviously [0, 0, 0, 0] (the zero background) is the only equilibrium solution (the critical point) of the system. The pulse solutions are the trajectories of the system (3.18) that start and end at the equilibrium (homoclinic orbits). Thus, the investigation of solitary pulses is mathematically equivalent to studying homoclinic solutions.

Structure of solitary waves
To investigate the structure of homoclinic solutions, we linearize the system in (3.18) near the critical point The characteristic equation of the matrixÃ on the right-hand side is given by Therefore, the values of p 2 coincide with the eigenvalues of the matrix −A. It is easy to see that Thus, the condition makes det A < 0, causing A to have eigenvalues of opposite signs, which is a necessary condition for the existence of homoclinic orbits. Indeed, in the case of A having eigenvalues of the opposite signs, the 4 × 4 matrixÃ has two pure imaginary eigenvalues (square roots of the negative eigenvalue of −A) and one negative, and one positive eigenvalues. Therefore the nonlinear system has one-dimensional stable and unstable manifolds, and a two-dimensional center manifold (corresponding to the imaginary eigenvalues). It was first noticed in [14] that the nonlinear system in (3.18) has Hamiltonian structure. If the kinetic, E, and potential, U, energies and the Hamiltonian, H, are introduce as follows then the system (3.18) takes the form Since the Hamiltonian is a conserved quantity, ∂ ζ H = 0, any trajectory issued from the critical point [0, 0, 0, 0] stays on the zero energy level surface H = 0 for all time. Note that the surface H = 0 is a 3D manifold in R 4 . The intersection of this 3D hypersurface with the hyperplanes Q 1 = 0 and M 1 = 0 is a curve Γ in the QM-plane (the figure-eight shaped curve on the left in Fig. 1 and 2). If for a given V there exits a homoclinic trajectory of (3.18), then on this trajectory E + U = 0 and since E ≥ 0, necessarily U ≤ 0. At the extrema of U its gradient is zero: By eliminating M from the equations above, we obtain the cubic equation whose roots are easily found: Thus, ∇U = 0 at the points thus producing the figure eight level curves. We already encountered this inequality above, see (4.21). After some algebra it can be rewritten as the following constraints on the traveling wave velocity: Thus, a possible velocity of the propagating pulse is bounded below.
5 Numerical study of solitary waves   A priori it is not clear why any homoclinic solution would possess this symmetry, and it is quite likely that there exist non-symmetric homoclinic orbits; we plan to investigate them elsewhere. The characteristic property of a time reversal orbit is that at the symmetry point Q 1 = M 1 = 0, and consequently the kinetic energy E must be zero; i.e., the symmetry point lies on the curve Γ, see (4.25). Moreover at the symmetry point the trajectory is orthogonal to Γ (for an illustration, see the QM diagrams on the left of Fig. 1 and Fig. 2).
Our search algorithm for finding solitary-wave solutions is based on the following minimization idea. If for a given value of the propagation velocity V there exists a homoclinic orbit with the time-reversal symmetry, then at some point both the kinetic and potential energies are zero. The algorithm takes the initial condition u 0 from a domain S on the  where E is the kinetic energy, u(ζ|u 0 ) is the solution of (3.18) with the initial condition u(0|u 0 ) = u 0 , recall that u = [Q, M, Q 1 , M 1 ]. The parameter τ 0 is the expected "width" of the pulse. Since E = 0 at ζ = 0 we take ζ > ζ 0 to obtain a nontrivial solution for the energy minimization problem. Computation of any particular value of min ζ 0 <ζ<ζ 0 +τ 0 E[u(ζ|u 0 )] involves a numerical solution of the nonlinear system of ODEs.
The search of the optimal initial datum is stochastic and is organized via a version of simulated annealing [31]. On each step the initial datum is obtained by sampling a random distribution with the density determined by the results of the previous step (see [28,32] for more detail). If Φ(V ) = 0 then there exists a homoclinic solutions with velocity V. When the kinetic energy possesses several local minima along the trajectory the corresponding solitary wave has the multi-hump structure. Figures 1 and 2 illustrate this phenomenon. The figure-eight shaped curves on the left correspond to the E = 0 cross-section of the potential energy landscapes; the curves inside the domains represent the Newtonian particle trajectories in the QM configuration space. The graphs on the right show the profiles of the corresponding solitary wave solutions. Fig. 1 illustrates a typical one-hump solution. In contrast, the trajectory shown in Fig. 2 has a point of the nearest approach to the boundary where the kinetic energy attains a local minimum. The resulting solution has a two-hump structure. Multi-hump solutions correspond to more complicated trajectories. Each of these trajectories has the return point at which it has the normal incidence with the E = 0 contour. For the fixed set of physical parameter values, the shape of the potential energy landscape is controlled by the pulse velocity V via the coefficients a 1 and a 2 in (3.15). We investigated numerically the set V of values of V which give rise to homoclinic orbits; in some sense one might think of these V s as the "spectrum" of the problem. For numerous applications with soliton-like solutions the velocity value is known to change continually (a continuous spectrum). However, for the Maxwell-Duffing model under consideration our numerical investigation demonstrates that the spectrum V contains both an interval of a continuous spectrum and a discrete subset of parameter values V for which a wave solution exists. One of the principal issues is to understand the correspondence between types of solitary wave solutions and values of V ∈ V.
We first investigated numerically the distribution density of the values of V , which give rise to homoclinic orbits. For all numerical computations of this section we adopted the following values of the nondimensional physical parameters: Our numerical computations show that on a rather small interval V c = [0.73, 0.7642] at the low end of the spectrum, every attempt of computing a homoclinic orbit was successful (20 orbits per bin). These results stay consistent with the refinement of the computational grid size δV.
All the solitary wave solutions in V c are of the one-hump variety; note, however that the single-hump solitons are not exclusively confined to the lower end of V. Elsewhere the spectrum density is very low, and the solitons are mostly of a multi-hump kind. Somewhat arbitrarily, we define a hump as a local maximum of the electric field e, which is at least 50% of the global maximum.
Next we studied the distribution of the different type of solitary wave solutions on the interval of velocities [0.73, 0.95]. The figure (Fig. 3(b)) gives a very clear idea of the placement of solitons according to the number of humps, which ranges from one to ten. Some typical soliton profiles for four-and eight-hump solutions with V = 0.84322 and V = 0.91416 respectively are collected in Fig. 4.
where N is the number of humps. As follows from this figure, in the log-log coordinates the energy increase is very well approximated by a linear function. The least square fit of the data from this figure shows that the energy increases approximately as a polynomial of fifth degree in V.

Formation, stability, and interaction of solitary waves: computer simulations
In this section through direct numerical simulations we study evolution of waves as well as wave interactions. We consider formation of solitary wave solutions from arbitrary initialboundary condition, stability of traveling waves under small perturbations and stability under strong perturbations due to wave collisions. In all numerical simulations of this section we use the same values of physical parameters (5.28) as in Sec. 5.
Numerically we solve the signaling problem for (2.7). In other words we give boundary conditions on either one or both ends of the spatial interval (0, L); as initial conditions we assign zero values for all the variables, which corresponds to propagation in a quiescent medium. For solving the initial-boundary value problem for the system in (2.7) we devised a simple fractional step numerical method.
Because the first two equations in (2.7) are hyperbolic PDEs while the rest are ODEs the choice of the fractional steps is extremely natural: on the first half-step we propagate the PDE part of the governing equations, and on the second half-step we march according to the system of ODEs. The resulting ODE system is solved by using the midpoint rule, while the PDEs are solved by the explicit McCormack method [33]. The midpoint rule and the McCormack method are both second order accurate. To increase the accuracy of the fractional step method we utilize the Strang split approximation [34], which results in the second order convergence of the final numerical scheme.
For many of the solitary-wave solutions discussed in Sec. 5 we ran direct numerical simulations on the model with these solitary waves as input pulses. All the waves tested, even the ones of a rather intricate shape, propagate with constant speed and without any shape distortion. See, for example, Fig. 6 where propagation of an eight-hump soliton is depicted.  Figure 10: Evolution of the 15 exp(−0.1τ 2 ) Gaussian. A "sharp" Gaussian quickly evolves into a near-solitary wave, leaving some disturbance in the wake. The speed of the nearsolitary wave is significantly higher than the speed of propagation of the radiation; thus the wave quickly leaves the disturbance behind This suggests that the solitary waves are (nonlinearly) stable with respect to numerical perturbations. We remark that although because of the scale of Fig. 6, the pulses appear rather singular, they are in fact completely smooth and numerically resolved. The numerical resolution of this computation is 8 mesh intervals per unit length, which provides about 60 mesh points per each hump of the traveling wave. Similarly, fine computational meshes are employed in all the simulations below.
The issue of stability can be addressed analytically by studying the linearization of the system of partial differential equations (2.7) about arbitrary traveling wave solutions and analyzing the corresponding linear evolution operator. Our analysis showed that this operator is skew-Hermitian in L 2 with the appropriate norm. Therefore the spectrum of the evolution operator is pure imaginary and the traveling wave solutions are neutrally linearly stable (see [28] for detail). To further elucidate the issue of stability we consider stability with respect to a finitesize perturbation in the initial wave shape. This situation is illustrated in Fig. 7. To the two-hump numerical soliton we add a rather substantial perturbation and employ the thus obtained functions as boundary data for the system of partial differential equations (2.7). As the result of evolution, the solution relaxes to the solitary wave shape followed by a low-amplitude "continuous radiation".
Stability with respect to strong perturbations due to collision of two traveling wave solutions is illustrated in Fig. 8. We take two solitary waves obtained by numerical solutions of ODEs and use these solutions as the boundary conditions for the PDEs. The left part of the Fig. 8  A soliton nature of solutions of (2.7) is further confirmed by the set of numerical simulations we present next. We consider propagation of solutions with the pulses in Fig. 9(a) given as a series of boundary conditions at the x = 0 boundary. The soliton of velocity V = 0.75 (see Fig. 1) propagates in a stable fashion, while its least-squares approximation by a Gaussian 9.65 exp −0.1t 2 approaches the soliton shape after shedding a small amount of residual continuous radiation. These time evolutions are not included for space saving (the energy dissipation curves for these cases show conservation of electromagnetic energy, see Fig. 9(b)).
In the next three figures (Fig. 10 -12) we present evolutions of the larger Gaussian pulses from Fig. 9(a). Evolution of the sharpest Gaussian (σ = √ 5) is displayed in Fig. 10. Very fast the solution forms a solitary wave that moves with constant velocity with no shape change. It is followed by low magnitude oscillations whose leading edge also moves with constant speed. During the evolution, the oscillatory part disperses more and more. This part of the solutions appears to be of a nonlinear nature; it will be studied separately. We note that although because of the scale of the figure, the pulse appears very sharp, it is in fact completely smooth with "width" about 20 and about 150 computational mesh points within the pulse.
The evolution of a wider Gaussian, σ = √ 10 (see Fig. 11) is similar with a very interesting distinction. Now the leading soliton is trailed by a slower low amplitude soliton. The latter is followed by low amplitude oscillations that again lag behind and disperse. The waves become more separated over time because the solitary waves with higher amplitude have higher velocities. Finally, the widest Gaussian, σ = 5 √ 2, develops into a train of four solitons, see Fig. 12.
To characterize the energy exchange between the propagating pulse and the medium, in Fig. 9(b) we present plots of the total electromagnetic energy as a function of time for all the input profiles from Fig. 9(a). For the soliton solution there is a dynamic equilibrium between the energy stored in the medium and the electromagnetic energy of the pulse. In case of the input impulse being not a soliton, the balance between the medium and the pulse is violated, which leads to the dissipation of the electromagnetic energy into the medium.

Concluding Remarks
In this paper we considered propagation of extremely short pulses in a nonlinear medium, which is characterized by both electric and magnetic resonance responses. Interaction of the electromagnetic field with the medium was described in the framework of the Maxwell-Duffing model. In particular we employed the classical Maxwell-Lorenz model for describing the magnetic resonance, [15]. For describing the interaction of the electric field component with the medium we used a generalized Maxwell-Lorentz model which takes into account cubic anharmonism of the polarization response (i.e., the Maxwell-Duffing system). Our findings demonstrate that the model supports a wide array of traveling-wave solutions. We investigated the structure and properties of these solutions through a combination of analysis and numerical modeling. We determined that the family of traveling-wave solutions is parameterized by one parameter, which is the velocity of a steady wave solution, normalized by the speed of light in vacuum. The spectrum V contains both an interval of a continuous spectrum and a discrete subset of parameter values for which a traveling-wave solution exists. Computer modeling demonstrated a multi-hump structure of these solutions. Their multi-hump nature suggests to characterize solitary wave solutions by a number of humps (types). All types are determined by not overlapping sets of velocities.
Direct numerical simulations showed that solitary-wave solutions are dynamically stable. This dynamical stability is consistent with the analysis of the system linearized about solitary wave solutions [28]. Stability of these solutions with respect to strong perturbations was studied by means of solitary wave collisions. Computer simulations indicated nearly elastic nature of scattering followed by a residual excessive radiation and a phase shift. In addition to traveling-wave solutions, numerical simulations demonstrated presence of another type of nonlinear oscillatory solutions with extended tail.
(Proposition 301), and by the Russian Foundation for Basic Research through grant 06-02-16406. Roytburd's work was partially supported by the National Science Foundation, while working at the Foundation. Part of his work was performed during a sabbatical leave at the Lawrence Berkeley National Laboratory. The authors would like to thank M. Stepanov for the enlightening discussions and for the valuable help in preparation of this manuscript.