1. Introduction

Four-wave mixing (FWM) is a fundamental nonlinear optical process, in which four photons interact by means of the third order nonlinear susceptibility of the medium. FWM can be observed in a wide range of materials, including optical fibers [1]. Modulational instability (MI), also known as Benjamin-Feir instability [2], is a special kind of FWM instability, and it is typical of nonlinear dispersive systems. It is known to exist in different branches of physics such as fluid-dynamics, plasma physics, solid-state physics and nonlinear optics [2, 3, 4]. It manifests itself as an exponential growth of two spectrally symmetric sidebands when launching a continuous wave (CW), by virtue of the interplay between Kerr nonlinearity and group velocity dispersion (GVD). In optical fibers it leads to the breakup of a CW or quasi-CW beam into a train of ultrashort pulses [5], and in the spatial case to the formation of patterns [6]. The long-term dynamics of MI can be very complicated, leading for instance to recurrence phenomena [7], and it is strictly connected to the existence of solitary waves [8].

In the context of nonlinear optics MI can be interpreted as degenerate FWM: two pump photons of identical frequencyω are converted into two new photons of different frequencies ω – δ and ω + δ, symmetric with respect toω, called respectively Stokes and anti-Stokes photons [1]. The existence and efficiency of frequency conversion by FWM crucially depends on the so-called phase-matching conditions [1], which are associated to the momentum conservation of the four photons involved in the process. Such conditions are determined in large part by the linear characteristics of the system, i.e. its GVD, even though the entire process is nonlinear in its nature.

Traditionally, FWM instabilities are analyzed in nonlinear optics in the framework of the nonlinear Schrödinger equation (NLSE) [1]. This equation (and its various extensions) provides the simplest, and in many cases even an accurate model for predicting the unstable frequency regions. In this work we show that, although NLSE (and its various generalizations, known as GNLSE, see [1, 9]) works extremely well in the context of glass optical fibers [10], there are physical situations in which, in addition to the term containing the second order derivative in time (which models the temporal dispersion), it is necessary to include a spatial dispersion coefficient, proportional to the Laplacian operator. This ‘spatiotemporal’ dispersion may be due, for example, to materials that exhibit the presence of elementary matter excitations, such as excitons [11]. The transverse components of the spatial dispersion term due to excitons contribute to the diffraction of the light beam propagating in the medium. This diffraction can be compensated for example by using common semiconductor waveguides [14]. However, the longitudinal component of the spatial dispersion cannot be compensated by the waveguide, and it introduces a new element in the dynamics of CWs and pulses, that, as we show here, strongly affects the nonlinear optical instabilities.

The structure of the paper is as follows. In section 2 we propose a propagation equation, by extending the (1+1)-NLSE to include both spatial and temporal dispersion coefficients, and we analyze the linear and nonlinear behavior of carrier plane waves in detail. Our equation must be distinguished from the non-paraxial NLSE [12, 13], due to a crucial difference in the sign of the spatial dispersion coefficient, which provides a totally different physical interpretation of the equation. Appendix A gives a physical justification of Eq. (1), derived in a system pumped in the proximity of an exciton resonance in a specially designed superlattice structure described in the literature. In section 3 we analyze the nonlinear parametric instabilities of the master equation, and we find some rather unconventional and unexpected properties. In section 4 we show the existence of an analytical family of bright and dark soliton solutions of the master equation, and we study their existence regions in the parameter space. Conclusions and final remarks are given in section 5.

2. Nonlinear evolution equation

In the present work we study the following evolution equation:

Eq. (1) is a NLSE with a spatially dispersive term proportional to the coefficient D, which acts along the longitudinal coordinate z. Parameter s = ±1 depending on whether we have anomalous or normal temporal dispersion, respectively. ψ is a field proportional to the slowly varying amplitude of the electric field. The last term models nonlinear cubic (Kerr) nonlinearity. We wish not to distract the reader from the main purpose of this work, thus the full physical justification of the model described by Eq. (1) shall be given in the Appendix A, where we study a waveguide of semiconductor material populated by exciton-polaritons possessing a negative effective mass, when optically excited slightly off-resonance. Eq. (1) is formally equivalent to the equation for non-paraxial solitons [12, 13], with the fundamental difference that parameter D can assume both signs, and that this parameter does not depend on the transverse profile of the field, see also discussions in Appendix A.

To explore the parametric instabilities of Eq. (1), the knowledge of the stationary states of the system is necessary. For plane waves, this is accomplished in the Fourier domain by the expansion ψ = ψ ₀ e ^ikz-_iωt, where ψ ₀ is a real constant field. This gives a dispersion relation 𝓓(k,ω,I) = 0, where I≡ψ ² ₀ is the dimensionless electric field intensity. In scattering problems usually one wants to keep ω real by solving the dispersion relation for complex k.

For D = 0, 𝓓 reduces to the well-known NLSE dispersion: 𝓓(k,ω,I)=-k-sω ²/2+I = 0, which has only one solution in terms of k = k(ω). Thus, only one steady state is physically excitable for NLSE. In the case D ≠ 0, we have 𝓓(k,ω,I)=-k+Dk ²-sω ²/2+I = 0, which has the two solutions k _α(ω,I)=[1+α(1-4D[I-sω ²/2])^1/2]/(2D), where α = ±1 defines the upper and lower branches respectively.

Two cases must be distinguished: D > 0 and D < 0. The former represents the most interesting situation for this paper, due to the unconventional instability content of Eq. (1). When D < 0 the modulational instabilities of Eq. (1) show no significant difference with the case of the non-paraxial NLSE [12, 13]. Therefore here we only treat the case D > 0. A physical interpretation of these two cases in terms of the sign of effective excitonic mass is considered in Appendix A.

Let us first examine the linear excitation case (I = 0). In anomalous dispersion conditions, i.e. for s = +1, the linear wave number k _L≡k(ω,I = 0) is always real and positive. Dispersion is made of two branches, k ₊ and k _-, separated by a k-bandgap, the extension of which is Δk = 1/D, see Fig. 1(a). For s = -1, i.e. in normal dispersion conditions, we have an inverted ω-bandgap (IBG), i.e. a range of frequencies in which k _α are purely real, see Fig. 1(b). The extension of this gap is Δω = [2/D]^1/2. Outside this bandgap, i.e. for |ω| > [2D]^-1/2, wave numbers become purely imaginary, and stationary states are not possible, due to strong field absorption or gain (depending on whether Im{k} > 0 or Im{k}<0 respectively). Therefore, in contrast to conventional bandgaps, the IBG is the only region which allows indefinite propagation of linear plane waves.

 

Fig. 1. Dispersion k _± of the two branches of the carrier wave in the linear excitation case (I = 0) for (a) anomalous dispersion regime (s = +1) and (b) normal dispersion regime (s = -1). D = 1 in both cases. Solid (dashed) lines indicate the real (imaginary) parts of k _±. Blue (red) lines refer to k ₊ (k _-). Gray shaded regions in (a) and (b) indicate the k-bandgap and the ω-inverted bandgap, respectively.

Download Full Size | PDF

Let us now examine the dispersion relations in the nonlinear regime (I ≠ 0). For s = +1, dispersion does not change in an essential way for I < I _t≡1/(4D). At this threshold intensity, the two branches touch each other. For I > I _t we have the situation shown in Fig. 2(a). A true nonlinearly-induced bandgap in ω appears in the dispersion (with extension Δω = [2/D-8I]^1/2), where eigenwaves possessing the full nonlinear wave number cannot propagate without being either suppressed or amplified. For s = -1, dispersion also has a qualitative change for I > I _t, see Fig. 2(b), where the closure of the IBG is shown, and field absorption (or gain) is dominant at all frequencies.

It is thus possible to go from a dispersion which exhibits a k-bandgap [Fig. 1(a)] to a dispersion which shows an ω-bandgap [Fig. 2(a)], or from a dispersion which exhibits an inverted ω-bandgap [Fig. 1(b)] to a dispersion showing total absorption/gain [Fig. 2(b)], by just increasing the incident intensity of the carrier plane wave. This effect is made possible by the simultaneous presence of both spatial and temporal dispersion terms in Eq. (1).

3. Parametric instabilities

Analysis of parametric wave generation in the system described by Eq. (1) follows the wellknown Benjamin-Feir scheme [2]. The stationary states found above are disturbed by a small perturbation field f (z,t), and the eigenvalues of the operator describing the small fluctuations are monitored in the frequency domain. Thus, in Eq. (1) we make the substitution ψ = (ψ ₀+ε f)e ^ik(ω,I)z-^iωt, where ε is a small parameter. The leading term in the perturbation theory leads to the following equations:

where with the bar we always indicate complex conjugation. The plus (minus) sign of α in Eq. (2) refers to the perturbation of the steady state in the upper (lower) dispersion branch. Defining the new vector v≡[f,f̄], we can write Eq. (2) in the short form

where Â≡diag(a,ā), a≡[1-4D(I-sω ²/2)]^1/2 is the complex discriminant, σ̂_z≡diag(1,-1) is the third Pauli matrix (which commutes with Â), and matrix operator M。 is given by

In Eq. (4), the dispersion operator induced by the perturbation is defined by D̂(i∂ _t)≡-isω∂ _t+(1/2)s∂ ² _t. We are now interested in the spatial evolution of the small perturbational plane waves (propagating on the carrier plane wave with frequency ω and intensity I) governed by Eq. (3). In order to extract the stability content of Eq. (3), we express v as $\mathbf{v}(z,t)=e^{i{\hat{\sigma }}_z\left[\kappa z-\delta t\right]}{\mathbf{w}}_{AS}(\kappa ,\delta )+e^{-i{\hat{\sigma }}_z\left[\overline{\kappa }z-\delta t\right]}{\mathbf{w}}_S(\kappa ,\delta )$ , where w _AS and w _S are respectively the anti-Stokes (high frequency) and the Stokes (low frequency) modes. By directly substituting into Eq. (3) we immediately obtain the secular equation for the complex wave number detuning of the perturbational wave κ, as a function of its (real) frequency detuning δ:

In the following we shall always assume that we work in regions of ω that make a real and positive, so that  = a1̂, where 1̂ is the identity matrix. For given D, s and ω, this gives a constraint on the maximum intensity, I < 1/(4D)+sω ²/2. This constraint is due to the fact that we are only interested in carrier waves which allow stationary states, and therefore cannot have any gain or loss in the dispersion. This condition will be clarified in section 4 when solitary wave solutions will be discussed. Equation (5) is a set of two quartic equations in κ (one equation for the upper branch, α = +1 and one for the lower branch, α = -1). However, Eq. (5) is unchanged under the discrete symmetry α→-α, a→ā, δ→-δ, and due to the fact that instability spectra are always symmetric with respect to the δ = 0 axis, it is sufficient to consider the case α = +1 without any loss of generality. The imaginary parts of the four solutions κ _j(δ) (j = {1,…,4}) of Eq. (5) determine the linear stability of the system. More specifically, if Im{κ _j} is negative, then frequency ω + δ will be unstable.

 

Fig. 2. Dispersion k _± of the two branches of the carrier wave in the nonlinear excitation case (I = 0.3 > I _t) for (a) anomalous dispersion regime (s = +1) and (b) normal dispersion regime (s = -1). D = 1 in both cases. Solid (dashed) lines indicate the real (imaginary) parts of k _±. Blue (red) lines refer to k ₊ (k-).

Download Full Size | PDF

Let us now review all possible unstable scenarios shown by the solutions of Eq. (5) for D > 0. First of all, let us explore the anomalous dispersion case, s = +1. In the following, we always pose for simplicity ω = 0. As we have seen, for light intensities I < I _t, the two dispersion branches of the carrier wave are always real, see Fig. 3(a,c). Absolute values of the imaginary parts of all unstable roots (i.e. the ones for which Im {κ _j(δ)}<0) as functions of δ are shown in Fig. 3(b,d), for two representative values of the spatial dispersion parameters D. For D < D _t≡1/(6I), one recovers the usual two lobes of modulational instability which is known for optical fibers [1], see Fig. 3(c). For increasing values of D, such that D > D _t, novel effects due to spatial dispersion term come into play, and the zero-detuning part of the roots becomes unstable, and the instability plot resembles the shape shown in Fig. 3(d).

 

Fig. 3. Parametric instabilities for anomalous dispersion regime, s = +1. (a,c) Nonlinear dispersion branches of the carrier wave, for D = 0.7<D _t < D _c [(a)], and for D = 0.85 > D _t < D _c [(c)]. Solid (dashed) lines are real (imaginary) parts of k _±. Blue (red) lines refer to k ₊ (k-). (b,d) Nonlinear instability gain Im{κ _j(δ)}<0 for the parameters described in (a) and (c) respectively. The carrier wave intensity is I = 0.22 in both cases. The values of D _t and D _c are 0.7576 and 1.1364 respectively.

Download Full Size | PDF

The unstable behavior in normal dispersion conditions (s = -1) is somewhat more interesting. As we have seen before, the full nonlinear carrier wave shows an IBG of extension Δω = [2/D-8I]^1/2, centered around ω = 0, see Fig. 4(a,b). In the limit D→0, this extension goes to infinity, together with the critical intensity I _t. For a fixed value of I, and for D∈[0,D _t], there are two unstable frequency regions near the boundaries of the IBG window. Outside this window, the imaginary parts of the carrier wave number is complex, and as a consequence there cannot be any steady state, due to photon absorption/gain. This is shown in Fig. 4(a,b). The two unstable regions tend to merge at D = D _t. For D∈[D _t,_{D c}], where D _c≡1/(4I) > D _t, a qualitative change occurs, and we have a situation similar to that shown in Fig. 4(d), where it is interesting to observe a multi-valued function for the gain plot. In the range of frequency detunings where the instability gain is multi-valued, the most negative imaginary part will dominate the exponential growth of perturbations.

For D > D _c the IBG window closes completely, and steady states over which unstable noise can grow cannot exist.

 

Fig. 4. Parametric instabilities for normal dispersion regime, s = -1. (a,c) Nonlinear dispersion branches of the carrier wave, for D = 0.7<D _t < D _c [(a)], and for D = 0.85 > D _t < D _c [(c)]. Solid (dashed) lines are real (imaginary) parts of k _±. Blue (red) lines refer to k ₊ (k-). (b,d) Nonlinear instability gain Im{κ _j(δ)κ<0 for the parameters described in (a) and (c) respectively. The carrier wave intensity is I = 0.22 in both cases. D _t and D _c are the same as in Fig. 3.

Download Full Size | PDF

4. Bright and dark solitons

The linear stability analysis performed in the previous section shows a non-trivial unstable content for plane waves travelling according to Eq. (1). It is well-known that parametric instabilities are connected to the existence of solitary waves in nonlinear systems. In this section we clarify this connection for Eq. (1), and we give explicit analytical solutions for bright solitons, which vanish at the space-time boundaries. Of course an essential requirement for bright solitary wave solutions of Eq. (1) to exist, is that they must be localized in both space (z) and time (t), because Kerr nonlinearity must balance both spatial and temporal dispersion terms in the equation.

In order to find bright soliton solutions for Eq. (1), we have used the following ansatz:

where A ₀ is the soliton amplitude, ν is the soliton velocity and z ₀ is the soliton width. After plugging expression (6) into Eq. (1), and imposing it to vanish for any value of z and t, one obtains

Again, α = ±1 in Eqs. (7)–(9) refers to upper and lower dispersion branches respectively. It is possible to see that solution (7)–(9) represents a family of 2-parameter solitons, parameterized by (ν,ω). Moreover, coefficients A ₀ and z ₀ in Eqs. (7)–(9) must be real, and this imposes restrictions on the values that v ad ω can assume. In addition, it is possible to prove with simple algebra that the following extra condition must hold:

 

Fig. 5. (a) Region of existence on the parameter space (ν,ω) for bright solitons (s = +1, D = 1), calculated by imposing reality conditions on Eqs. (7)–(9), plus the extra condition Eq. (10). (b) Regions of existence on the parameter space (ν,ω) for dark (black + gray) solitons (s = -1, D = 1), calculated by imposing reality conditions on Eqs. (12–14), plus the extra condition sνω = ±α|ν||ω|, valid for the black and the gray solitons respectively.

Download Full Size | PDF

Constraints given above are quite restrictive. Importantly, it can be demonstrated with a bit of elementary algebra that, analogously to the case of nonlinear fiber optics [1, 9], bright solitons can never exist for s = -1, irrespective of the values of the parameters (ν,ω), in any branch of the dispersion relation. The regions of existence for the family of bright solitons found above is displayed in Fig. 5 for the upper branch of the dispersion k ₊. This region of the parameter space is, however, profoundly different from the analogous region found for fiber solitons [1, 9]. Black regions in Fig. 5 are those regions in the parameter space where reality of (7)–(9) and constraint (10) are simultaneously satisfied. For solitons living on the lower branch k-, the region of existence (not shown here) is obtained by taking the mirror image of Fig. 5 with respect of the ν = 0 axis. A noteworthy feature is the onset of a velocity gap, which forbids the existence of bright solitons inside the velocity gap -[2D]^1/2<ν < [2D]^1/2. This gap is filled with dark solitons, as we shall see shortly. Solitons possessing a smaller velocity v (but such that |v| > [2D]^1/2), exist for a wide range of frequencies (positive or negative, depending on the sign of v and on the specific dispersion branch k _α), while large velocity solutions can only exist near a narrow region around ω = 0.

We have also found dark soliton solutions for Eq. (1). We chose the following ansatz:

where θ is an angular parameter which regulates the depth of the dark soliton. For θ = 0 one has a so-called ‘black soliton’, which has a vanishing minimum of the amplitude, to be distinguished from the ‘gray soliton’, which has θ ≠ 0. This distinction is made not only for classification purposes, but because in general black and gray solitons may have different dynamical and stability properties [9]. We have found by direct substitution of Eq. (11) into Eq. (1) that parameters of solution Eq. (11) are given by the following expressions:

Note that, given ν and ω, the value of θ is not arbitrary, and is determined by Eq. (15). Also, it is clear from Eq. (15) that the black soliton solutions (θ = 0) and gray soliton solutions (θ ≠ 0) have different existence regions in the parameter space (ν,ω). In fact the black soliton exists under the condition sνω = α|ν||ω|, while the gray soliton exists if sνω = -α|ν||ω|, where α = ±1 correspond to the upper and lower branch of the dispersion, respectively. This can be also seen in Fig. 5(b), where the existence region for the ‘gray’ soliton (in the upper dispersion branch) is shown by the light gray area, while the existence region for the ‘black’ soliton is shown with a black area. Note that the above dark solitons fill the velocity gap which appeared in the parameter space of the bright solitons family, see Fig. 5(a).

5. Conclusions

In conclusion, we have analyzed a model of Nonlinear Schrödinger equation which includes both spatially and temporally dispersive terms. This model is mathematically equivalent to the non-paraxial model described in [12, 13], but it is physically very different, mainly due to the difference in sign for the spatial dispersion parameter D. It was found that qualitatively new parametric processes are possible in presence of excitons with a negative effective mass. In particular, unconventional parametric instabilities are present in the normal dispersion regime for a certain range of values of the input intensities or the spatial dispersion parameter. Analytical expressions have been found for bright and dark soliton solutions of Eq. (1). The stability of such objects, which is a highly non-trivial issue, will be rigorously analyzed in a separate publication. Finally we have demonstrated in Appendix A that semiconductor waveguides made of excitonic materials possessing negative exciton masses might realize the conditions for the model described by Eq. (1) to be valid with a good degree of approximation.

A. Derivation of Eq. (1) for an off-resonance excitonic transition

In this Appendix we demonstrate that the dynamics of intense light propagation near an excitonic resonance in a 1D semiconductor waveguide in presence of Kerr nonlinearity, is regulated by Eq. (1). The starting point of our derivation is the dielectric function of an isotropic nonmagnetic material under excitation of monochromatic light of frequency ω, which, in the proximity of an excitonic resonance reads:

In the above expression, ε _b is the background dielectric constant, ω̃₀≡[ω ² ₀-γ ²]^1/2, with ω ₀ and γ the exciton resonance and damping, respectively, Δ≡Ω ² _c/(2ω̃₀), Ω_c being a frequency proportional to the oscillator strength for the coherent exciton-photon interaction. The parameter Γ≡h̄/(2M*_x), where M*_x is the effective exciton mass, takes into account the excitonic spatial dispersion through the wavevector dependence of the excitonic transition, ω = ω ₀+Γ|k|². Note that in expression (16) only the resonant part of the dielectric function is taken into account, while the non-resonant part is generally totally negligible [15]. Solving Maxwell equations with the constitutive relation D=ε _k,_ω E, with ε _k,ω given by Eq. (16), yields the following equation for transverse modes:

which represents the dispersion of exciton-polariton waves, i.e., the true mixed modes propagating in the medium, resulting from the interaction between excitons and the retarded electromagnetic field. As first pointed out by Pekar [16], the inclusion of the spatial dispersion in the excitonic dielectric function Eq. (16), leads in general to a twofold solution of Eq. (17) at a given ω: an upper transverse polariton branch ω _UT=ω _UT(|k|) and a lower transverse one ω _LT=ω _LT(|k|). Polaritons in bulk semiconductors (as well as in confined structures) are well-established concepts, and the existence of the above branches has been demonstrated for various crystals in the proximity of the Wannier exciton resonances (see e.g. Ref. [17]). In semiconductor crystals, a Wannier exciton has typically a positive effective mass (M*_x > 0), and the ω _LT(|k|) branch is thus characterized by a positive group velocity ν _gr≡∂ω _LT/∂|k| > 0. In some special cases, however, the lower branch ω _LT(|k|) is characterized by a negative group velocity, ν _gr<0, or, in other words, a negative effective mass, M*_x <0.

Let us now assume that we pump an intense laser beam into a semiconductor waveguide which provides an effective transverse confinement, therefore making the problem effectively one-dimensional. k = [0,0,k] is the wave vector for the effective 1D problem, the only non-vanishing component of which points along the longitudinal propagation direction Z. Then the linear constitutive relation between the Fourier components of polarization P _k,_ω and electric field E _k,_ω can be written as

where χ ^X,^W are the exciton and waveguide contributions to the susceptibility, respectively. Note that χ ^W does not depend on k, due to the absence of spatial non-locality in conventional semiconductor single-mode waveguides [9].

The linear propagation of Fourier components of the electric field (due to the exciton part of the susceptibility only) is regulated by the following equation:

where E ^(±) _k,_ω are respectively the forward and backward components of the electric field (which is supposed to be linearly polarized along one of the waveguide’s principal axis), and c is the speed of light in vacuum. Taking equations only for the forward component corresponds to performing the well-known slowly-varying envelope approximation (SVEA), which is one of the cornerstones of nonlinear optics [1].

Substitution of expression (16) into Eq. (19) yields the following equation for the forward propagation:

where δω≡ω-ω̃₀. We now assume that the central frequency of the input laser is tuned slightly off-resonance, but such that |δω|≫|-Γk ²+i γ| and |δω|≫Δ, i.e. on the ‘tail’ of lower exciton-polariton branch. In this physical situation we do not need to take into account the coupled nature of the electric and polarization fields, and the knowledge of the dynamics of the electric field is sufficient to accurately model light propagation. Also, the oscillator damping constant will be negligible in this regime, and in the following we set it equal to zero for simplicity. In this limit Eq. (20) can be approximated by

By taking the inverse Fourier transform in k and ω of Eq. (21), and including the waveguide dispersion contribution as given by Eq. (18), we obtain

where Z is the longitudinal propagation direction, T is time, A≡A(Z,T) is the slow-varying envelope of the linearly polarized electric field, ω ₀ is the central pulse frequency, and γ_NL is the Kerr nonlinear coefficient exhibited by the semiconductor waveguide. The operator

models the waveguide temporal dispersion, and contains terms such as the group velocity [proportional to the β ₁ coefficient in Eq. (23)] and the group velocity dispersion (GVD, proportional to the β ₂ coefficient), plus higher-order dispersion terms (β _j≥3) if necessary [18]. For our purposes, it is sufficient to make the expansion D̂^W(i∂ _T)≃iβ ₁ ∂ _T-(1/2)β ₂ ∂ ² _T, where β ₁,₂ are the inverse of the group velocity and the GVD parameters respectively. β ₁ can always be eliminated by a proper choice of moving coordinates without loss of generality (mixed derivatives that can arise are typically very small due to the fact that the exciton transverse dispersion is much smaller than the waveguide dispersion), and we therefore consider only β ₂ in the calculation. Eq. (22) is the conventional generalized nonlinear Schroedinger equation (GNLSE, [1, 9, 18]), but with the additional contribution of the longitudinal spatial dispersion of the exciton, proportional to the second order derivative in space.

By introducing dimensionless variables z≡Z/z ₀, t≡T/t ₀, and ψ≡A/A ₀, by imposing the scalings z ₀=t ² ₀/|β ₂|, A ₀=[γ _NL z ₀]^-1/2, t ₀ being the input pulse duration, and defining $s\equiv \mathrm{sgn}(-{\beta }_2)$ and $D\equiv -\frac{\left[\sqrt{{\epsilon }_b}\Gamma \Delta {\tilde{\omega }}_0\right]}{\left[2\delta {\omega }^{2}c{z}_0\right]}$, we finally obtain our master Eq. (1).

Two physically different situations can be realized for Eq. (1). The first regime occurs for D < 0, i.e. when the effective exciton mass M*_x is positive. We name this case normal spatial dispersion regime. This is a standard case which can be realized, for instance, for the very common semiconductor materials exhibiting conventional bulk Wannier exciton resonances. This case does not present particular interest for us, due to its equivalence to the non-paraxial NLSE case which can be found in the literature [12, 13], and it will be not analyzed in the present paper. The second regime, when D > 0 (which we name anomalous spatial dispersion regime), is characterized by a negative effective mass M*_x < 0 and is the most interesting one for this paper, due to the new and unconventional unstable content which can result from Eq. (1). This case is physically realizable, for instance, in some organic molecular crystals [19, 20, 21]. However, here we consider as our representative example a structure made of a strained ZnCdSe/ZnSe superlattice, investigated in a series of recent works, where the authors clearly show that the experimental data (taken at room temperature) are fully consistent with a negative exciton-polariton dispersion, induced by high compressive strain [22]. The excitonic feature for the type-II superlattice of [22] is located in the visible range with a central value of ω ₀≃2.51 eV, i.e. λ ₀≃493 nm. The background dielectric constant for this material is ε _b≃8.66, the coupling frequency is Ω_c≃85.4 meV, and the exciton mass assumes the value M*_x = -0.65m ₀, where m ₀ is the free electron mass. With these values of parameters, one can calculate Γ≃-8.9 × 10^-5 m²/sec, and Δ≃1.5 meV. By choosing t ₀≃200 fsec, and |β ₂| = 0.1 psec²/m, and pumping the input pulse with a detuning equal to δω≃ω0/250, we obtain a value of the dimensionless spatial dispersion parameter equal to D≃1.18 × 10^-2, which is relatively large and should be observable by means of the nonlinear propagation regimes described in the previous sections. The Kerr nonlinear coefficient for ZnSe can be estimated to assume the large value γNL~0.8 m^-1W^-1 for a pulse transverse effective area of about A _eff~100 µm² [23].

One final thing to be noticed (which we have already mentioned) is the formal similarity of Eq. (1) with the governing equation of non-paraxial solitons intensively investigated by Chamorro-Posada, see Eq. (3) of [12]. However, one must notice the profound physical difference between the non-paraxial term, proportional to the second order derivative in the longitudinal coordinate, and our analogous term in Eq. (1). In fact, the non-paraxial spatial dispersion coefficient is given in our notation by D _nonpar≡-(λ ₀/w ₀)²/(4π ² ε _b), where w ₀ is the transverse pulse width (A _eff≡w ² ₀). D _nonpar is always negative and becomes appreciable only for w ₀→λ ₀. This is in striking contrast with the physical origin of coefficient D in Eq. (1), attributed to exciton-polaritons, which can have both signs and is independent from the transverse width of the pulse. This means that the non-paraxial contribution is negligible when compared with the exciton-polariton contribution (i.e. |D _nonpar|≪|D|) whenever w ² ₀≫λ ² ₀/[4π ² εb D]^1/2, a condition that in the example described above translates to the relation w ² ₀≫0.06 µm², which is easily satisfied for our choice of A _eff. It is interesting to note that non-paraxiality can be provided also by other physical situations, other than narrow beams, like for instance oblique propagation of the beam [13]. However, in the present work these interpretations do not apply, due to the fact that our model has one spatial coordinate only. Moreover, solutions Eqs. (6) and (11) are qualitatively different from the solutions found in [12, 13]. In these latter works, for D _nonpar→0, solutions collapse to the expressions for the NLSE soliton, while this is not the case for D→0 in Eqs. (6) and (11).

Our model is also suitable of further extensions, such as the inclusion of dispersive terms of order higher than the second in Eq. (23) (that become important for short pulse durations), and of more realistic nonlinear terms, such as Raman and saturable nonlinearities.