Asymptotic dynamics of three-dimensional bipolar ultrashort electromagnetic pulses in an array of semiconductor carbon nanotubes

We study the propagation of three-dimensional bipolar ultrashort electromagnetic pulses in an array of semiconductor carbon nanotubes at times much longer than the pulse duration, yet still shorter than the relaxation time in the system. The interaction of the electromagnetic field with the electronic subsystem of the medium is described by means of Maxwell's equations, taking into account the field inhomogeneity along the nanotube axis beyond the approximation of slowly varying amplitudes and phases. A model is proposed for the analysis of the dynamics of an electromagnetic pulse in the form of an effective equation for the vector potential of the field. Our numerical analysis demonstrates the possibility of a satisfactory description of the evolution of the pulse field at large times by means of a three-dimensional generalization of the sine-Gordon and double sine-Gordon equations.


Introduction
Carbon nanotubes-quasi-one-dimensional macromolecules of carbon [1][2][3] -are now considered to be one of the most promising materials with a high potential of applicability in the development of the elemental base for modern electronics.From the point of view of applications in optoelectronics, it is of particular interest to study the properties of nanotubes with respect to the peculiarity of their electronic structure.The non-parabolicity of the dispersion law for the conduction electrons of nanotubes [1][2][3] (i.e. the energy dependency on the quasi-momentum) determines the essential nonlinearity of their response to the action of the electromagnetic field [4].This circumstance provides the principal possibility of observing various unique electromagnetic phenomena in structures based on carbon nanotubes, including nonlinear diffraction and self-focusing of laser beams [5,6].
The possibility of stable propagation of ultrashort electromagnetic pulses in arrays of carbon nanotubes was first theoretically established in [32,33] in a one-dimensional model, and was subsequently studied in a two-dimensional model [34][35][36][37][38] in the framework of the homogeneous field approximation along the axes of nanotubes.The approach used in the above works did not take into account the effects associated with the inhomogeneity of the field along the nanotube axes; therefore, in order to obtain a more realistic description of field evolution, the mathematical model was generalized to the case of two and three spatial dimensions, and was supplemented by equations that take into account the redistribution of the charge density in the system under the influence of the field inhomogeneity along the axis of the nanotubes.As a result, the peculiarities of propagation of solitary electromagnetic waves in arrays of carbon nanotubes were studied both in two-dimensional [39,40] and three-dimensional [41][42][43] models, taking into account the effect of localization of the field along the nanotube axes.In the course of numerical experiments carried out within the framework of these studies, the possibility of stable propagation of bipolar ultrashort electromagnetic pulses in the form of bipolar "breather-like" light bullets at distances significantly exceeding the characteristic dimensions of the pulses along the direction of their motion was confirmed.
To date, a large amount of information has already been accumulated as a result of extensive studies of various aspects of the interaction of extremely short laser pulses with arrays of nanotubes, but there are still numerous questions that require further clarification.In particular, from the point of view of possible practical applications, it is instrumental to study the evolution of the electromagnetic pulse field in an array of nanotubes at time intervals that significantly exceed the characteristic pulse duration, but that are still shorter than the relaxation time in the system.Such a problem of considering the dynamics of a pulse at times several orders of magnitude greater than its duration naturally arises in the case of the propagation of extremely short pulses with a characteristic duration ∆t pulse ∼ 10 −15 − 10 −14 s in the medium with the relaxation time t rel ∼ 10 −12 − 10 −11 s, which is quite achievable with modern technologies for fabricating nanotube structures.A study of the asymptotic dynamics of an ultimately short pulse at large times was carried out earlier in [44] for a one-dimensional model in the homogeneous field approximation along the nanotube axis.However, the results obtained within that simplified framework cannot always be automatically extrapolated to the behavior of a pulse in a model containing more than one spatial dimension.Peculiarities of the behavior of a pulse in a nonone-dimensional geometry can occur owing to the influence of the transverse effects associated with the diffraction spreading of the wave packet.In this connection, it seems expedient to solve the problem of modeling the evolution of parameters of an extremely short laser pulse in an array of carbon semiconductor nanotubes, taking into account the three-dimensional localization of the field.

System configuration and key assumptions
In this paper, we consider the propagation of a solitary electromagnetic wave in the bulk array of single-walled semiconductor carbon nanotubes (CNTs) embedded in a homogeneous dielectric medium.The nanotubes considered here are of the "zigzag" type (m, 0), where the number m (not a multiple of three for semiconductor nanotubes) determines the nanotube radius R = mb √ 3/2π, and b = 1.42 × 10 −8 cm is the distance between neighboring carbon atoms [1][2][3].The CNTs are supposed to be arranged in such a way that their axes are parallel to the common x-axis, and distances between adjacent nanotubes are much larger than their diameter, so that we neglect the interaction between them.This latter assumption allows us to consider the system as an electrically quasi-one-dimensional one, in which electron tunneling between neighboring nanotubes can be neglected, and electrical conductivity is possible only along the axis of nanotubes.
We define the configuration of the system in such a way that the pulse propagates through the nanotube array in a direction perpendicular to their axes (that is, for definiteness along the z-axis), and the electric component of the wave field E = {E, 0, 0} is collinear with the x-axis (see Fig. 1).We note that for a wide range of values of the system parameters, the characteristic distance at which an appreciable change in the field of a bipolar electromagnetic pulse occurs is significantly greater than the distance between neighboring nanotubes, and also the characteristic length of the conduction electron path along the axis of the nanotubes.For example, for a nanotube radius R ≈ 2.74 × 10 −8 cm and m = 7, the characteristic distance between them (sufficient to exclude the overlap of the electron wave functions of neighboring nanotubes)-even if substantially exceeding the value of the radius R-will be negligibly small in comparison with the characteristic wavelength of the electromagnetic radiation in the system.In particular, in the case of infrared electromagnetic radiation this assumption remains valid.With this approximation, the nanotube array, which is a discrete structure at the microscopic level, can nevertheless be considered as a continuous medium at the considered wavelength scale of electromagnetic radiation propagating in it.
Another important assumption that we adopt in this paper concerns the magnitudes of the time duration of the electromagnetic pulse, ∆t pulse , the relaxation time of the conduction current along the nanotube axis, t rel , and also the time interval for observing the evolution of the field in the system, ∆t.Specifically, we assume that the observation time ∆t is much longer than the pulse duration ∆t pulse , but still shorter than the relaxation time t rel , namely ∆t pulse ≪ ∆t < t rel .
Given the orientation of the coordinate system axes relative to the nanotube axis chosen above (see Fig. 1), the electron energy spectrum for carbon nanotubes takes the form where the electron quasimomentum is p = {p x , s}, s being an integer characterizing the momentum quantization along the perimeter of the nanotube, s = 1, 2, . . ., m, m is the number of hexagonal carbon cycles, forming the circumference of a nanotube, γ 0 is the overlap integral, and where

Complete set of equations
As we showed earlier (see, for example [42,43] and references therein), the evolution of the electromagnetic field in an array of nanotubes is described by Maxwell's equations for the vector and scalar potentials, A = {A, 0, 0} and φ, together with the continuity equation [45,46] for the current density j = { j, 0, 0}, which is calculated according to the approach proposed in [47,48].As a result, the set of evolutionary equations consists of three of them: for the vector and scalar potentials, A = {A, 0, 0} and φ, respectively, and for the conduction electron density n.
The equation describing the evolution of the vector potential of a self-consistent electromagnetic field in an array of nanotubes has the form [42] Here, we use the following notation: Ψ = Aed x /(c ) is the projection of the dimensionless vector potential onto the nanotube axis; Φ = φ √ εed x /(c ) is the dimensionless scalar potential; ε is the average relative permittivity of the sample (see, for example [49]); η = n/n bias = η(ξ, υ, ζ, τ) is the reduced (dimensionless) distribution of conduction electron density, with n the instantaneous value of the conduction electrons concentration at an arbitrary point of the sample bulk; n bias is the concentration of conduction electrons in the sample in the absence of an electromagnetic field (that is, essentially a constant component of the concentration).The coefficients G r are the dimensionless quantities that decrease with increasing r [42]; τ = ω 0 t/ √ ε is the scaled time; ξ = xω 0 /c, υ = yω 0 /c and ζ = zω 0 /c are the dimensionless coordinates.The quantity ω 0 = 2|e|d x −1 √ πγ 0 n bias has the meaning of the characteristic frequency of the electron subsystem of nanotubes (e < 0 is the electron charge).
The equation describing the spatio-temporal evolution of the concentration of conduction electrons in an array of nanotubes has the form [39][40][41][42][43] with α = d x γ 0 √ ε/c .This equation determines the change in the electron concentration under the effect of a self-consistent electromagnetic field in the array of nanotubes.We emphasize that the inhomogeneity (localization) of the field only along the x-axis can cause non-stationarity and dynamic spatial inhomogeneity of the electron concentration in the sample, due to the presence of conductivity only along the axis of the nanotubes.In this case, the inhomogeneity of the field along the directions orthogonal to the axes of the nanotubes will not contribute to the redistribution of the electron concentration due to the negligible overlap of the electron wave functions of neighboring nanotubes and the absence of conductivity in the yOz-plane.
As was shown earlier in [39][40][41][42][43], the equation describing the evolution of the scalar potential of a field in an array of nanotubes has the form where β = 1/α = c / d x γ 0 √ ε , and η 0 = n 0 /n bias = η(ξ, υ, ζ, τ 0 ) is the reduced (dimensionless) value of the concentration of conduction electrons at a given point in space at the initial instant of time τ 0 in the absence of the field.In the simplest particular case corresponding to an initially homogeneous sample (n 0 = n bias ), we have η 0 = 1.
Thus, the evolution of the field in the array of nanotubes is described by the set of Eqs. ( 2)-( 4).These equations describe a self-consistent system of physical parameters, Ψ, Φ, and η, the dynamics of which reflect the mutual influence of the electromagnetic field and the electronic subsystem of the nanotubes array (self-consistent light-matter interaction).

Effective equation for the vector potential
When considering the process of propagation of an electromagnetic pulse in an array of nanotubes at times considerably exceeding the pulse duration, special attention should be given to an adequate consideration of the physical factors affecting the evolution of the wave packet parameters.An analysis of the degree of influence of these factors on the pulse dynamics in certain cases can allow us to optimize the mathematical model describing the process.In particular, one of such factors is the inhomogeneity of the electromagnetic field along the nanotube axis.As follows from the set of Eqs. ( 2)-( 4), the limiting (three-dimensional) localization of the field in the array of nanotubes can cause a redistribution of the conduction electron density.Thus, the propagation of the electromagnetic pulse occurs when interacting with the dynamic inhomogeneities of the medium, which are induced by the field of the pulse itself.
As shown by the numerical simulation performed earlier (see [39,41,42]), the electron density differences (dynamic inhomogeneities) formed during the passage of electromagnetic pulses in the sample have magnitude of the order of one percent relative to the initial equilibrium concentration n 0 .In this case, there was no dramatic change in the dynamics of the pulses with respect to the results obtained within the framework of the model limited to the approximation of a homogeneous field along the nanotube axis (see, for example [34]).This can be explained by the fact that the higher the speed of the electromagnetic pulse, the shorter the time during which the pulse affects the electronic subsystem of the array of nanotubes, the less is the formation of dynamic electron density inhomogeneities.At the same time, the higher the pulse speed, and the smaller the time during which the pulse propagates in the region of the space containing the dynamic inhomogeneities of the electron concentration, which cause the pulse field to be adjusted to the changed properties of the medium.Thus, under the conditions considered in this paper, the time of interaction between the field and the perturbations of the medium is not sufficient to generate a noticeable change in the shape of the pulse.Thus, in the case of ultrashort electromagnetic pulses (when the condition ∆t pulse ≪ t rel is satisfied), the non-stationarity of the conduction electron density distribution can be neglected.Proceeding from these considerations, we will assume that the distribution of the conduction electron concentration in the sample remains approximately constant, that is, ∂η/∂τ = 0.If the concentration of conduction electrons in the nanotubes array initially is uniform throughout the sample volume, i.e., there were no regions of high or low electron concentration (that is, n 0 = n bias ), then it can be assumed that η ≈ η 0 = 1 throughout the sample during the entire observation time t < t rel .Under the condition η = η 0 , the right-hand side of Eq. ( 4) vanishes, and without loss of generality (see [45,46]), we can safely impose that Φ(ξ, υ, ζ, τ) = 0.
Thus, taking into account the above arguments, one can exclude Eqs. ( 3) and ( 4) for the concentration of conduction electrons and the scalar potential from further consideration.As a result, taking into account the assumptions made about the short duration of the electromagnetic pulse propagating in the array of nanotubes, with sufficient accuracy, the evolution of the field in a homogeneous array of nanotubes (which does not initially contain regions of increased or decreased electron concentration) can be described by the only equation resulting from Eq. ( 2) with η = 1 and Φ = 0: This equation is an effective equation describing the evolution of the field in an array of nanotubes in the case of propagation of "fast" three-dimensional extremely short electromagnetic pulses.Since Eq. ( 5) contains all the terms of the series G r -where the subscript r corresponds to the mode number in the expansion of the conduction electron energy (1) in the Fourier series (see, for example, Eq. ( 4) in [42]), this equation can be regarded as the "complete" effective equation.
Calculations show that the coefficients G r decrease rapidly with increasing values of r (see [42] and references therein).Therefore, in a number of cases, for a qualitative approximate description of the evolution of the electromagnetic field in the system under consideration, we can restrict ourselves to only two terms with r = {1, 2} under the summation sign in Eq. ( 5), thereby leading to Equation ( 6) can be regarded as a "reduced" effective equation, which happens to be a threedimensional generalization of the Double sine-Gordon equation [50].To the best of our knowledge, this represents the first attempt in using the sine-Gordon equations for the three-dimensional study of such pulse propagation in arrays of CNTs.It is worth adding that using model (6), as opposed to model ( 5), provides significant savings from the computational standpoint.The case of the most radical reduction of the total effective Eq. ( 5) deserves special consideration.In the case of a significant excess of the absolute value of the quantity G 1 over the absolute value of the quantity G 2 , we can restrict the summation to the first term under the summation sign, having obtained a non-one-dimensional generalization of the sine-Gordon equation: where σ = √ G 1 (G 1 > 0).It is worth adding that in the present case, the amplitude of the 3D breathers is not small.Thus, linearizing the sine-Gordon equations is not an option.However, a Taylor expansion with respect to Ψ, up to quintic terms, could be considered.Doing so would enable the use 3D nonlinear equations with cubic-quintic nonlinear terms, in which the classical variational approximation could be applied to 3D breathers.Note that these equations with the cubic-quintic nonlinearity are known to support stable solitons with embedded vorticity [51].
Following the reasoning proposed in [44], the reduced effective equation for the vector potential (7) can be replaced by an equation of the form where σ Σ is determined by the empirical formula [44] The model simplification in the form of the 3D generalization of the sine-Gordon equation conveniently offers the ability to estimate analytically the asympotic dynamics of the electromagnetic pulses (see Appendix).The complete effective Eq. ( 5) and the three versions of the reduced equation, ( 6)-( 8), will be used below as comparative models for studying the dynamics of the ultrashort electromagnetic pulse in an array of semiconductor carbon nanotubes over long times.Laslty, it is worth adding that our effective simplified model in the form of Eq. ( 5) (or alternatively in the form of Eqs. ( 6)-( 8)) is primarily applicable to the description of the central part of the pulse, rather than its tail.

Characteristics of the electromagnetic pulse field
As a quantity illustrating the localization of an electromagnetic pulse in space, we will use the energy characteristic of the field, which is proportional to the bulk density of the field energy.
For such a value, it is convenient to take the square of the electric field strength along the axis of the nanotubes, x .This quantity will be referred to as "intensity" for convenience.Using the well-known formula E = −c −1 ∂A/∂t − ∇φ (see, e.g.[45,46]), and also taking into account that the vector potential has a nonzero component only along the nanotube axis (see the description of the system configuration above), one can simply express the electric field vector E = {E, 0, 0} as where Taking into account Eq. ( 10), the energy characteristic of the field chosen above can be represented as where I 0 = E 2 0 .To illustrate the position of the electromagnetic pulse in space at an arbitrary time instant τ, we define the dimensionless averaged characteristic coordinates of the pulse as the coordinates of the pulse centroid ("center of mass"), ξ pulse , υ pulse , ζ pulse .To this end we calculate the first order of the moments of the pulse intensity, following the approach proposed in [52]: The speed of the electromagnetic pulse is then defined as As quantitative characteristics of the localization of the pulse field in space, we calculate the longitudinal half-width of the pulse λ ζ (along the ζ-axis), as well as the transverse half-widths λ ξ and λ υ along the ξ-axis and υ-axis, correspondingly.We employ the M-squared method for characterizing laser beams [52].We adopt this approach to the situation corresponding to the limiting (three-dimensional) localization of the field, therefore defining the half-widths of the pulse by means of the second moment of the pulse intensity profile as follows: For simplicity, we introduce the generalized transverse characteristic half-width of the pulse, λ ⊥ , as:

Initial condition: shape of the electromagnetic pulse
We will investigate the evolution of the field of an electromagnetic pulse in an array of nanotubes numerically, solving the equation for the vector potential (see Eqs. ( 5)-( 8)).As an initial condition, we take an "instantaneous snapshot" of the distribution of the dimensionless projection of the field vector potential on the ξ-axis (the axis of the nanotubes) at the time instant τ = τ 0 : where the function Ψ (ζ, τ 0 ) determines the distribution (along the Oz direction) of the projection of the dimensionless vector potential on the nanotube axis (Ox); and the function Ψ ⊥ (ξ, υ) represents the initial distribution of the field in the ξOυ-plane orthogonal to the direction of propagation of the pulse.We select the profile Ψ (ζ, τ 0 ) corresponding to the breather solution of the sine-Gordon equation, namely, to the non-topological oscillating soliton [50], where with U = u/v 0 being the ratio between the initial propagation velocity u of the sine-Gordon breather in the one-dimensional approximation along the ζ-axis, and the speed of light in the medium v 0 = c/ √ ε.The quantity ζ 0 is the dimensionless coordinate of the center of mass of breather (21) along the ζ-axis at the time instant τ = τ 0 , and Ω = ω B /ω 0 , where ω B is the self-oscillation frequency of the breather (0 < Ω < 1).As for the value of the quantity σ in Eqs. ( 22) and ( 23), we choose σ = √ G 1 for certainty and comparability of the results of modeling the evolution of the same initial pulse in different models (see Eqs. ( 5)-( 8)); the values of the coefficients G r are calculated using Eq. ( 4) from [42].The justification of the choice of the initial pulse profile ( 21) is given in our previous papers (e.g., see [40,42]).We choose the transverse profile of the intensity of the pulse field with the corresponding Gaussian distribution, which is justified by the high degree of adequacy of this description to real situations widely represented in various fields of physics and technology [5,6,[53][54][55]: where ξ 0 and υ 0 are the dimensionless coordinates of the "center of mass" of the pulse at the initial instant of time τ 0 , and w 0 is the value characterizing the transverse localization of the electromagnetic pulse at the initial instant of time (the value of this characteristic is proportional to the initial value λ ⊥ (0) determined by Eq. ( 19), but not exactly equal to it).
Taking into account the expressions ( 20)-( 24), the projection of the electric field strength (10) of the electromagnetic pulse field onto the axis of nanotubes at the initial instant of time has the form where The electromagnetic pulse corresponding to the breathing solution of the form (21) of the sine-Gordon equation for the vector potential, during the propagation, performs internal oscillations: the shape of its profile.The projection of the electric field strength of such a pulse takes both positive and negative values, so a pulse of this type is said to be "bipolar." The duration of the electromagnetic pulse with a profile whose shape can be approximately described by Eq. ( 21) is determined by the function cosh µ that determines the envelope of a given wave packet.The duration of the electromagnetic pulse ∆t pulse is defined as the time interval during which the values of the "running" envelope (i.e., the function cosh µ) measured at a fixed point will exceed half of its peak amplitude value (i.e., have values above the level 0.5).Given that the width of the function cosh µ at the level of 0.5 (FWHM -full width of half maximum) is ∆µ 0.5 = 2 ln 2 + √ 3 , we get the value of the duration of the electromagnetic pulse:

System parameters
As an environment for modeling the propagation of an extremely short electromagnetic pulse, we chose an array of nanotubes of the zig-zag type (m, 0): m = 7, γ 0 = 2.7 eV, b = 1.42 × 10 −8 cm, d x ≈ 2.13 × 10 −8 cm, n bias = 10 18 cm −3 at temperature T = 293 K.We assume that the array of CNTs is embedded in a dielectric matrix with the effective dielectric constant ε = 4.The dimensionless parameter U (see Eqs. ( 21)-( 23)) was varied within the interval U ∈ (0.5; 0.999).For U < 0.5 and further decrease in the value of this parameter, the longitudinal width of the electromagnetic pulse begins to approach the value of the distance traveled by the pulse over a duration ∼ t rel , which is of no significant practical interest.The case corresponding to the condition U > 0.999 is not described in the present paper because of the limitations of the numerical scheme we used.
The dimensionless parameter of the frequency of the internal oscillations Ω of the initial pulse (21) was varied over the interval Ω ∈ (0.1; 0.9).As this parameter decreases, the characteristic width of the pulse along the ζ-axis decreases, while for Ω < 0.5, the variation in the form of the profile is insignificant.For Ω > 0.9, the width of the pulse along the ζ-axis becomes comparable with the dimensions of the numerical grid.The parameter w 0 of the characteristic transverse width of the pulse was varied in the range from 1.0 to 10.0.
The set of differential Eqs. ( 5)-( 8) does not have an exact analytical solution in the general case.Therefore, we carried out numerical simulations to study the propagation of an electromagnetic pulse in an array of CNTs.To solve each of the equations with initial conditions ( 20)-( 25), we used an explicit finite-difference three-layered scheme of the cross type described in [56][57][58] and adapted by us for the three-dimensional model, using the approach developed and detailed in [40].As a result of the numerical experiment, we have modeled the evolution of the vector potential Ψ(ξ, υ, ζ, τ), and also calculated the distribution of the field strength and intensity using Eqs.( 10) and (11).
We emphasize that the mathematical model used in this work is valid for observation times t shorter than the relaxation time t rel of the electronic subsystem, since under the condition t < t rel , only the evolution of the electromagnetic field in the system can be adequately evaluated, neglecting damping due to collisions of electrons with irregularities in the crystal structure of the nanotube array.For example, with t rel ≈ 10 −11 s-which can be achieved with modern sampling technologies-the limiting distance traveled by light in the medium under consideration is ∆z ≈ ct rel / √ ε = 0.15 cm.We performed a numerical experiment on the propagation of an electromagnetic pulse over the time interval ∆τ = 3 ×10 2 (corresponding to dimensional time ∆t = ∆τ √ ε/ω 0 ≈ 8.4 ×10 −12 s), which is close in the order of magnitude to the value of the relaxation time t rel indicated above.The purpose of the numerical experiment was to find out whether the solutions of equations ( 5)-( 8) retain the properties inherent to the initial condition in the form of a three-dimensional pulse with a longitudinal profile of the breather (see (21)).Specifically, we were interested in whether the individuality of the electromagnetic pulse persists in spite of the phenomena of dispersive and diffractional spreading, whether the wave packet remains bipolar, and whether it retains the properties of the breather with respect to periodic changes in the shape and amplitude over long times.
For definiteness and comparability, the results of modeling the propagation of an electromagnetic pulse in an array of carbon nanotubes is presented below for the following values of the dimensionless parameters of the model: Ω = 0.5, w 0 = 5.0, U = 0.95 (u = 1.425 × 10 10 cm/s).The maximum amplitude of the electric field of the pulse in this case has a value |E x | max ≈ 1.170×10 7 V/cm (see Eq. ( 25)), and its characteristic duration is ∆t pulse ≈ 2.9×10 −14 s (see Eq. ( 26)).The characteristic frequency of the electromagnetic field of a given sub-cycle pulse can be estimated as ω pulse ≈ 2π/∆t pulse ≈ 2.14 × 10 14 s −1 .
The results of the comparison of the numerical solution of the full effective Eq. ( 5) with the solutions of the reduced Eqs.( 6)-( 8), with the same initial condition ( 20)-( 25) for all these four cases, are given below.

The complete effective equation
Figures 2-4 show the results of the numerical simulations of the propagation of a threedimensional electromagnetic pulse within the framework of the model represented by the full effective Eq. ( 5). Figure 2 shows the distribution of the intensity of the field, I(ξ, υ 0 , ζ, τ) = I 0 (∂Ψ/∂τ) 2 , in the ξOζ cross section (at υ = υ 0 ) at various instances of the dimensionless time τ.The color scale is assigned to different values of the ratio I/I 0 : the yellow areas correspond to the maximum values, and the dark-blue areas correspond to the minimum values.The horizontal and vertical axes of the graph represent the dimensionless coordinates ξ = xω 0 /c and ζ = zω 0 /c.Given the system parameter values selected above, the unit along the ξand ζ-axes corresponds to the distance ≈ 4.2 × 10 −4 cm.We note that here we give the distribution of the quantity I/I 0 in the ξOζ-plane only, since the pattern of the distribution of the intensity of the field in the υOζ-plane is qualitatively similar.
We draw attention to the fact that in order to save computation time and to simplify the visual representation, we used a numerical scheme with periodic boundary conditions (see the details in our previous work [40]).For this reason, in Figure 2 (and also in Figure 3), the positions of the wave packet at the later instants ("e"-"h") belong to the same spatial interval along the ζ-axis that corresponds to the positions of the wave packet at the earlier time instants ("a"-"d").It can be seen from Fig. 2 that the wave packet, having overcome a distance substantially exceeding its characteristic size along the direction of propagation (along the positive ζ direction), retains its individuality without undergoing a decay due to diffraction and dispersive spreading.
Figure 3 shows the evolution of the distribution profile of the electric component of the field E x (Eq.( 10)) of a given wave packet and the intensity I along the ζ-axis passing through the point with coordinates ξ = ξ 0 and υ = υ 0 , corresponding to the initial position of the "center of mass" of the pulse in the ξOζ-plane at the same instants of dimensionless time τ as in Fig. 1.In this case, the center of mass of the pulse in the ξOζ-plane is not displaced, i.e. according to the simulation results, ξ pulse = ξ 0 and υ pulse = υ 0 .
It can be seen from Fig. 3 that during the observation time, E x is alternating between positive and negative values, that is, the electromagnetic pulse remains bipolar.The nature of the change in the configuration of the wave packet can be illustrated by the dependency of the scaled amplitude of the intensity I max /I 0 = max (∂Ψ/∂τ) 2 on the dimensionless time τ. Figure 4 shows the time dependency of the quantity I max /I 0 for the same values of the system parameters as in Figs. 2 and 3.
As a result of numerical analysis, it is established that the wave packet preserves its individuality in the course of propagation, without undergoing significant spreading and damping; the electric field remains alternating over a time interval substantially exceeding the characteristic duration of the wave packet.Thus, this three-dimensional electromagnetic pulse has the properties of a nontopological soliton-breather, determined by the initial conditions ( 21) and (25).This circumstance allows one, in a certain sense, to consider a three-dimensional bipolar electro-  magnetic extremely short pulse propagating in an array of carbon nanotubes as being a soliton.
We note that the results presented here correspond to the situation in which the electromagnetic pulse propagates in an essentially nonlinear regime.We emphasize that in this case, the absolute values of the extrema of the projection of the dimensionless vector potential Ψ on the nanotubes axis reach values of the order of unity (that is, the condition | sin Ψ| ≪ 1 is not satisfied in the general case) and, consequently, Eq. ( 5) cannot be linearized with the chosen values of the system parameters.From the physical point of view, this means that the response j (electric current density) depends in an essentially nonlinear manner on the vector potential (and hence on the electric field strength) of the self-consistent field of the electromagnetic pulse (see formula (3) from our previous work [42]).

3D generalization of the double sine-Gordon equation
As noted above, the coefficients G r (see Eq. ( 5)) decrease rapidly with increasing values of the index r.According to Eq. ( 4) from our previous work [42], and with the system parameters chosen above, we obtain the following values for the first five terms of the series G r (up to the second decimal place): G 1 ≈ 0.91, G 2 ≈ 0.33, G 3 ≈ 0.18, G 4 ≈ 0.12, and G 5 ≈ 0.08.Thus, in fact, the third term is already significantly smaller than the first one, which justifies the actual truncation to obtain the full effective equation, so that Eq. ( 6) is a satisfactory approximation of the original Eq. ( 5).
We have numerically solved the truncated effective Eq. ( 6) in the form of a three-dimensional generalized double sine-Gordon equation with initial conditions ( 20)-( 25) on the dimensionless time interval ranging from 0 to τ = 300.As a result, we established that the electromagnetic pulse propagates over a distance substantially exceeding its characteristic size, and does not undergo appreciable diffraction, transverse or dispersive longitudinal spreading.At the same time, the solution of the truncated effective Eq. ( 6) is a bipolar solitary electromagnetic wave with a "breathing" periodically repeating its shape profile.In other words, the evolution of the pulse field in model ( 6) is qualitatively similar to the behavior of the electromagnetic pulse in the model of the full effective Eq. ( 5).Thus, the evolution of the field of a three-dimensional electromagnetic pulse (with the same initial parameters) in an array of nanotubes is described by Eqs. ( 5) and ( 6) in a qualitative way, but some quantitative differences in the values of the parameters of the steady-state solution in the form of a stably propagating bipolar solitary wave.

3D generalization of the sine-Gordon equation
We further simulated the propagation of the electromagnetic pulse within the framework of the model represented by the truncated effective Eq. ( 7), for the same values of the system parameters and with the same initial condition as before.It follows from the results of the calculations that the pulse of model ( 7) at long time scales also possesses the properties of a breather-a bipolar solitary wave-with a periodically changing "breathing" amplitude.Model (8) with the corrected coefficient σ 2 Σ in front of the sine leads to a similar picture of the evolution of the pulse field: after the transient stage (a short-term increase in the degree of transverse field localization with a corresponding increase in amplitude), the pulse reaches a stable propagation mode, qualitatively similar to that in the framework of model (5).

Comparison of the pulse dynamics in the framework of complete and reduced models
To compare the behavior of the solutions obtained in the framework of different models represented by Eqs. ( 5)- (8), we analyze the time dependencies of the scaled amplitude of the maximum field intensity I max /I 0 , and the ratio of the transverse half-width of the pulse λ ⊥ (τ) to its initial value λ ⊥ (0) (see Eq. ( 19)).Figures 5 and 6 show, respectively, the time dependencies of the values of I max /I 0 and λ ⊥ /λ ⊥ (0), obtained as a result of solving the effective Eqs. ( 5)-( 8) at identical values of the system parameters and the same initial condition given by formulas ( 20)- (25).It can be seen from Fig. 5 that a noticeable difference in the behavior of the solution of Eq. ( 8), with the modified coefficient in front of the sine, from the behavior of the solutions within the remaining truncated models ( 6) and (7), is a more pronounced outgrowth of the field amplitude with the field focusing at the initial stage, which is similar to the behavior of the pulse within the framework of the model represented by the full effective Eq. ( 5).The amplitude and transverse width of the pulse reach values close to the initial values at the end of the transient process again, and the pulse dynamics reaches the mode corresponding to the stable propagation of a breather, which is the common result obtained within all the analyzed models ( 5)- (8).We also note that the truncated model in the form of Eq. ( 8) with the modified coefficient σ Σ , in the most accurate manner (with respect to the full Eq.( 5)) approximates the variations of the transverse width of the electromagnetic pulse due to diffraction.This conclusion follows from the analysis of the curves λ ⊥ (τ) in Fig. 6.Thus, comparing the numerical simulation results shown in Figs. 5 and 6, it can be seen  5)) and truncated models represented by Eqs. ( 6)-( 8): (a) red -full effective Eq. ( 5); magenta -truncated model in the form of the three-dimensional generalization of the sine-Gordon Eq. ( 8) containing the coefficient σ Σ ; (b) green -truncated model in the form of the three-dimensional generalization of the double sine-Gordon Eq. ( 6); blue -truncated model in the form of the three-dimensional generalization of the sine-Gordon equation (7).
Fig. 6.Dependency of the scaled transverse half-width of the pulse λ ⊥ /λ ⊥ (0) on the dimensionless time τ during the propagation of the electromagnetic pulse under the same conditions as in Fig. 5: red -full effective Eq. ( 5); green -truncated model in the form of a three-dimensional generalization of the double sine-Gordon equation ( 6); blue and magenta -truncated model in the form of the three-dimensional generalization of the sine-Gordon equation ( 7) and ( 8), respectively.
that the model represented by the three-dimensional generalization of the sine-Gordon and the double sine-Gordon equations can adequately approximate the dynamics of the pulse at long time scales.We note that the truncated model in the form of a generalization of the sine-Gordon equation with the corrected coefficient calculated from the empirical formula (9), is the simplest and satisfactory model for the asymptotic analysis of the electromagnetic pulse dynamics at the long time scales.

Conclusions
Key results of this work may be summarized as follows: i) An effective equation (see Eq. ( 5)) is presented that determines the evolution of the electromagnetic field in an array of semiconductor carbon nanotubes, taking into account the limiting (three-dimensional) field localization.
ii) The evolution of an extremely short three-dimensional bipolar electromagnetic pulse in an array of nanotubes is analyzed by means of numerical simulations over a long time interval, which significantly exceeds the pulse durations, and still shorter than the relaxation time.
The possibility of stable propagation of a bipolar electromagnetic pulse with limiting (three-dimensional) field localization in an array of nanotubes at large time scales is confirmed.
iii) The possibility of approximating the effective Eq. ( 5) by simpler models in the form of three-dimensional generalizations of the double sine-Gordon equation ( 6) and the sine-Gordon equation (( 7)-( 8)) is suggested and numerically verified.
iv) The dynamics of the pulse is compared in the framework of the models associated with the full effective Eq. ( 5) and the truncated effective Eqs. ( 6)- (8).It is established that the truncated effective equations are adequate models for the description of the dynamics of a three-dimensional electromagnetic pulse in an array of carbon nanotubes over long times.

Possibilities of analytical consideration. Asymptotic analysis of the electromagnetic pulse dynamics
The applicability of Eqs. ( 6)-( 8) proposed above opens the possibility for further progress in the analytical study of the problem of the dynamics of extremely short electromagnetic pulses in carbon nanotubes.Consider the illustration of an analytical approach to the asymptotic analysis of the field evolution over long times using the example of the model described by Eq. ( 7) (generalization to Eq. ( 6) is simple).After re-normalization of time τ → σ τ and coordinates ξ → σ ξ, υ → σ υ, and ζ → σ ζ , Eq. ( 7) can be considered as the Euler-Lagrange equation for the Lagrangian density, The change to the cylindrical coordinate system ( υ = r cos θ, ξ = r sin θ), as well as the assumption that the extremely short pulse preserves the cylindrical symmetry in the course of its propagation (dΨ/dθ → 0) allows us to represent the Lagrangian density (27) in the form Next, consider, according to the Whitham approach [59], a solution in the form of a 2π-kink, in which the "fast" and "slow" variables are clearly distinguished: where ρ and Φ are the fast and slow variables, respectively.Passing further to the coordinates of the light cone for τ and ζ , and also averaging over ( τ − ζ ), we obtain the corresponding Lagrangian density: where The set of equations for the quantities ρ and Φ corresponding to the Lagrangian density (30) has a hydrodynamic form [59]: where the notation ϕ = −2Φ is introduced for convenience of presentation.An analogy with the hydrodynamics of an ideal fluid arises when the right-hand side of Eq. ( 31) is equated to zero.In this case, the quantity ρ has the meaning of density, the quantity ϕ corresponds to the potential of the velocity field, and the quantity −2/ρ 2 ≡ H(ρ) plays the role of enthalpy [60].
Considering the relationship between enthalpy and pressure P in the form H(ρ) = ∫ ρ −1 dP, we consider the stability of our solution with respect to long-wavelength transverse perturbations of such a kind that the right-hand side of Eq. ( 31) can be neglected.In this case, taking into account the results from hydrodynamics, it is necessary to satisfy the condition P > 0 for the stability of the solution.Thus, the hydrodynamic analogy for the model described by Eq. ( 7) allows one to establish the stability of its solutions in the form of kinks over arbitrary (long) time scales.Development and application of this analytical approach with respect to solitary waves of other types, including breathers, is also of interest, and we consider it as a subject for further research.

Fig. 1 .
Fig. 1.Schematic diagram of the considered system with the associated coordinate system.

Fig. 4 .
Fig. 4. Dependency of the ratio I max /I 0 on the dimensionless time τ during the propagation of the electromagnetic pulse shown in Figs. 2 and 3.The dimensionless time τ = ω 0 t/ √ ε is plotted along the horizontal axis.

Fig. 5 .
Fig. 5. Dependency of the scaled maximum energy density of the pulse I max /I 0 on the dimensionless time τ along the ζ -axis in the course of the electromagnetic pulse propagation in the full (Eq.(5)) and truncated models represented by Eqs.(6)-(8): (a) red -full effective Eq.(5); magenta -truncated model in the form of the three-dimensional generalization of the sine-Gordon Eq. (8) containing the coefficient σ Σ ; (b) green -truncated model in the form of the three-dimensional generalization of the double sine-Gordon Eq. (6); blue -truncated model in the form of the three-dimensional generalization of the sine-Gordon equation(7).

Funding
A. V. Zhukov and R. Bouffanais are financially supported by the SUTD-MIT International Design Centre (IDC).N. N. Rosanov acknowledges the support from the Russian Foundation for Basic Research, Grant 19-02-00312, and from the Foundation for the Support of Leading Universities of the Russian Federation (Grant 074-U01).M. B. Belonenko acknowledges support from the Russian Foundation for Fundamental Research.