Scalar and vector electromagnetic solitary waves in nonlinear hyperbolic media

TL;DR

This paper explores the propagation of electromagnetic solitary waves in hyperbolic nonlinear media, deriving exact solutions and analyzing how hyperbolic dispersion influences wave stability and unique propagation characteristics.

Contribution

It provides analytical solutions for scalar and vector solitary waves in hyperbolic nonlinear media, highlighting the role of hyperbolic dispersion in wave stability and preventing singularities.

Findings

Hyperbolic dispersion adds a degree of freedom to control wave singularities.

Exact analytical forms of solitary waves are derived using Hirota bilinear method.

Solitary waves exhibit unique propagation properties not seen in conventional media.

Abstract

In this paper, we investigate the problem of electromagnetic wave propagation in hyperbolic nonlinear media. To address this problem, we consider the scalar hyperbolic nonlinear Schr\"odinger system and its coupled version, namely hyperbolic Manakov type equations. These hyperbolic systems are shown to be non-integrable. Then, we examine the propagation properties of both the scalar and vector electromagnetic solitary waves by deriving their exact analytical forms through the Hirota bilinear method. A detailed analysis shows that the presence of hyperbolic transverse dispersion provides an additional degree of freedom to prevent the formation of singularity in both the scalar and vector solitary wave structures in this hyperbolic nonlinear media. Besides this, we realize that the solitary waves in this media possess fascinating propagation properties which cannot be observed in conventional nonlinear media. We believe that the present study will be very useful in analyzing electromagnetic wave propagation in hyperbolic nonlinear metamaterials.

Equations

$$\frac{k_{x}^{2}+k_{y}^{2}}{\varepsilon_{o}}+\frac{k_{z}^{2}}{\varepsilon_{e}}=\bigg{(}\frac{\omega}{c}\bigg{)}^{2},$$

$$i\frac{\partial q_{j}}{\partial z}+\frac{1}{2}\bigg{(}\frac{\partial^{2}q_{j}}{\partial x^{2}}+\beta\frac{\partial^{2}q_{j}}{\partial y^{2}}\bigg{)}+\gamma_{j}(\lvert q_{1}\rvert^{2}+\sigma\lvert q_{2}\rvert^{2})q_{j}=0,~{}q_{j}\equiv q_{j}(x,y,z),~{}j=1,2.$$

$$i\frac{\partial q}{\partial z}+\frac{1}{2}\bigg{(}\frac{\partial^{2}q}{\partial x^{2}}+\beta\frac{\partial^{2}q}{\partial y^{2}}\bigg{)}+\gamma\lvert q\rvert^{2}q=0,~{}q\equiv q(x,y,z).$$

$$(i{D}_z+\frac{1}{2}(D_x^2+\beta D_y^2)-\lambda )g\cdot f=0,(D_x^2+\beta D_y^2-2\lambda )f\cdot f=2\gamma |g{|}^2.$$

$$\displaystyle D_{z}^{m}D_{x}^{n}(a\cdot b)=\bigg{(}\frac{\partial}{\partial z}-\frac{\partial}{\partial z^{\prime}}\bigg{)}^{m}\bigg{(}\frac{\partial}{\partial x}-\frac{\partial}{\partial x^{\prime}}\bigg{)}^{n}a(z,x)b(z^{\prime},x^{\prime})_{\lvert z=z^{\prime},~{}x=x^{\prime}}.$$

$$q(x,y,z)=\frac{ϵg_1}{1+{ϵ}^2{f}_2}=\frac{{\alpha }_1{e}^{{\eta }_1}}{1+e^{{\eta }_1+{\eta }_1^{\ast }+\delta }},$$

$${\eta }_1=k_{1}x+\beta k_{2}y+\frac{i}{2}(k_1^2+{\beta }^3{k}_2^2)z,e^{\delta }=\frac{\gamma |{\alpha }_1{|}^2}{(k_1+k_1^{\ast }{)}^2+{\beta }^3(k_2+k_2^{\ast }{)}^2}.$$

$$q(x,y,z)=A{e}^{i({\eta }_{1I}+\theta )}\text{mbox}sech({\eta }_{1R}+\frac{\delta }{2}).$$

$$q(x,y,z)$$

$$Z_1$$

$$\tau =\pm \sqrt{\frac{(k_1^2+{\beta }^3{k}_2^2)}{\gamma (-2+Z_1+Z_1^{\ast })}}.$$

$$\displaystyle q(x,y,z)=\frac{\tau}{2}e^{i\theta}\big{[}(1+Z_{1})-(1-Z_{1})\tanh\frac{\xi_{1}}{2}\big{]}.$$

$$\lvert q\rvert^{2}=\tau^{2}\bigg{(}1-X\mbox{sech}^{2}\frac{\xi_{1}}{2}\bigg{)},$$

$$X=\frac{(k_1^2+{\beta }^3{k}_2^2{)}^2}{|(k_1^2+{\beta }^3{k}_2^2)+2i(k_1{l}_1+k_2{l}_2{\beta }^3-{\Omega }_1){|}^2}.$$

$$D_1{g}^{(j)}\cdot f=0,D_{2}f\cdot f=2\gamma (|g^{(1)}{|}^2+\sigma |g^{(2)}{|}^2),j=1,2,$$

$$q_j(x,y,z)=\frac{ϵg_1^{(j)}}{1+{ϵ}^2{f}_2}=\frac{{\alpha }_1^{(j)}{e}^{{\eta }_1}}{1+e^{{\eta }_1+{\eta }_1^{\ast }+R_1}},j=1,2,$$

$${\eta }_1=k_{1}x+\beta k_{2}y+\frac{i}{2}(k_1^2+{\beta }^3{k}_2^2)z,e^{R_1}=\frac{\gamma (|{\alpha }_1^{(1)}{|}^2+\sigma |{\alpha }_1^{(2)}{|}^2)}{(k_1+k_1^{\ast }{)}^2+{\beta }^3(k_2+k_2^{\ast }{)}^2}.$$

$$q_j(x,y,z)={\hat{A}}_j(\frac{k_{1R}^2+{\beta }^3{k}_{2R}^2}{\gamma }{)}^{1/2}{e}^{i{\eta }_{1I}}\text{mbox}sech({\eta }_{1R}+\frac{R_1}{2}),$$

$${\hat{A}}_j=\frac{{\alpha }_1^{(j)}}{\sqrt{|{\alpha }_1^{(1)}{|}^2+\sigma |{\alpha }_1^{(2)}{|}^2}},\frac{R_1}{2}=\frac{1}{2}\log \frac{\gamma (|{\alpha }_1^{(1)}{|}^2+\sigma |{\alpha }_1^{(2)}{|}^2)}{(k_1+k_1^{\ast }{)}^2+{\beta }^3(k_2+k_2^{\ast }{)}^2},$$

$$\displaystyle\eta_{1R}=k_{1R}x+\beta k_{2R}y-\big{(}k_{1R}k_{1I}+\beta^{3}k_{2R}k_{2I}\big{)}z,$$

$$\displaystyle\eta_{1I}=k_{1I}x+\beta k_{2I}y+\frac{1}{2}\big{(}(k_{1R}^{2}-k_{1I}^{2})+\beta^{3}(k_{2R}^{2}-k_{2I}^{2})\big{)}z.$$

$$q_1=\frac{ϵg_1^{(1)}}{1+{ϵ}^2{f}_2}=\frac{{\alpha }_1^{(1)}{e}^{{\eta }_1}}{1+e^{{\eta }_1+{\eta }_1^{\ast }+R_2}},$$

$$q_2=g_0\frac{(1+{ϵ}^2{g}_2^{(2)})}{1+{ϵ}^2{f}_2}=\tau e^{i\gamma \sigma {\tau }^{2}z}\frac{(1-e^{{\eta }_1+{\eta }_1^{\ast }+R_2+{\chi }_1})}{1+e^{{\eta }_1+{\eta }_1^{\ast }+R_2}},$$

$${\eta }_1=k_{1}x+\beta k_{2}y+\frac{i}{2}(k_1^2+{\beta }^3{k}_2^2)z-i\lambda z,\lambda =-\gamma {\tau }^2\sigma ,$$

$$e^{{\chi }_1}=\frac{k_1(k_1+k_1^{\ast })+{\beta }^3{k}_2(k_2+k_2^{\ast })}{k_1^{\ast }(k_1+k_1^{\ast })+{\beta }^3{k}_2^{\ast }(k_2+k_2^{\ast })},e^{R_2}=\frac{\gamma |{\alpha }_1^{(1)}{|}^2{\chi }_2}{(k_1+k_1^{\ast }{)}^2+{\beta }^3(k_2+k_2^{\ast }{)}^2},$$

$$\displaystyle\chi_{2}=\bigg{(}1+\frac{\sigma\tau^{2}\gamma\big{(}(k_{1}+k_{1}^{*})^{2}+\beta^{3}(k_{2}+k_{2}^{*})^{2}\big{)}}{(k_{1}(k_{1}+k_{1}^{*})+\beta^{3}k_{2}(k_{2}+k_{2}^{*}))(k_{1}^{*}(k_{1}+k_{1}^{*})+\beta^{3}k_{2}^{*}(k_{2}+k_{2}^{*}))}\bigg{)}^{-1}.$$

$$q_1(x,y,z)=A(\frac{k_{1R}^2+{\beta }^3{k}_{2R}^2}{\gamma }{)}^{1/2}{e}^{i({\eta }_{1I}+\theta )}\text{mbox}sech({\eta }_{1R}+\frac{R_2}{2}),$$

$$\displaystyle q_{2}(x,y,z)=\tau e^{i(\gamma\sigma\tau^{2}z+\phi+\pi)}\bigg{[}\cos\phi\tanh(\eta_{1R}+\frac{R_{2}}{2})+i\sin\phi\bigg{]},$$

$$\displaystyle\eta_{1R}=k_{1R}x+\beta k_{2R}y-\big{(}k_{1R}k_{1I}+\beta^{3}k_{2R}k_{2I}\big{)}z,$$

$$\displaystyle\eta_{1I}=k_{1I}x+\beta k_{2I}y+\frac{1}{2}\big{[}k_{1R}^{2}-k_{1I}^{2}+\beta^{3}(k_{2R}^{2}-k_{2I}^{2})-\lambda\big{]}z,$$

$$\displaystyle\phi=\frac{1}{2}\tan^{-1}\bigg{(}\frac{2(k_{1R}^{2}+\beta^{3}k_{2R}^{2})(k_{1R}k_{1I}+\beta^{3}k_{2R}k_{2I})}{(k_{1R}^{2}+\beta^{3}k_{2R}^{2})^{2}-(k_{1R}k_{1I}+\beta^{3}k_{2R}k_{2I})^{2}}\bigg{)},~{}~{}e^{i\theta}=\sqrt{\frac{\alpha_{1}}{\alpha_{1}^{*}}},$$

$$\displaystyle A=\bigg{[}1+\frac{\sigma\tau^{2}\gamma\big{(}(k_{1}+k_{1}^{*})^{2}+\beta^{3}(k_{2}+k_{2}^{*})^{2}\big{)}}{(k_{1}(k_{1}+k_{1}^{*})+\beta^{3}k_{2}(k_{2}+k_{2}^{*}))(k_{1}^{*}(k_{1}+k_{1}^{*})+\beta^{3}k_{2}^{*}(k_{2}+k_{2}^{*}))}\bigg{]}^{1/2}.$$

$$q_j={\tau }_j{e}^{i{\theta }_j}\frac{(1+Z_j{e}^{{\xi }_1})}{1+e^{{\xi }_1}},j=1,2,{\xi }_1=k_{1}x+\beta k_{2}y-{\Omega }_{1}z,$$

$$Z_1=-\frac{i({\Omega }_1-l_1{k}_1-{\beta }^3{l}_2{k}_2)+\frac{1}{2}(k_1^2+{\beta }^3{k}_2^2)}{-i({\Omega }_1-l_1{k}_1-{\beta }^3{l}_2{k}_2)+\frac{1}{2}(k_1^2+{\beta }^3{k}_2^2)},$$

$$Z_2=-\frac{i({\Omega }_1-l_3{k}_1-{\beta }^3{l}_4{k}_2)+\frac{1}{2}(k_1^2+{\beta }^3{k}_2^2)}{-i({\Omega }_1-l_3{k}_1-{\beta }^3{l}_4{k}_2)+\frac{1}{2}(k_1^2+{\beta }^3{k}_2^2)}.$$

$${\tau }_1=\pm \sqrt{\frac{\gamma \sigma {\tau }_2^2(-2+Z_2+Z_2^{\ast })-(k_1^2+{\beta }^3{k}_2^2)}{\gamma (2-Z_1-Z_1^{\ast })}}.$$

$$\displaystyle q_{1}=\frac{\tau_{1}}{2}e^{i\theta_{1}}\bigg{(}(1+Z_{1})-(1-Z_{1})\tanh\frac{\xi_{1}}{2}\bigg{)},$$

Full text

††journal: Chaos, Solitons & Fractals

1 Introduction

Recently, engineering the propagation of light using artifical material is one of the hot areas of research in nonlinear optics [1]. In particular, metamaterials or left-handed materials, having double-negative refractive indices [2, 3], are highly beneficial to construct the meta-devices with desirable properties. An interesting aspect of this material is that one can tune its properties depending on a particular need or application. Such materials exhibit properties that are never or rarely observed in nature. To enhance these properties further, in recent years, highly anisotropic media with hyperbolic dispersion, which is characterized by the effective electric and magnetic tensors, have been introduced to the metamaterials [1]. As a result, the resultant composite materials exhibit unique properties including enhanced superlensing effects [1], strong enhancement of spontaneous emission [1], and a large Purcell factor [4, 5, 6]. These synthetic materials have also applications in sub-wavelength optics [1], heat transport [7], acoustics [8], and cosmology [9, 10]. The latter fascinating properties emerge from by engineering the following dispersion relation of extraordinary transverse magnetic polarized waves. It is represented by

[TABLE]

where, $k_{x}$, $k_{y}$, and $k_{z}$ are the $x$, $y$ and $z$ components of the wave vector $k$, respectively, $\omega$ is the frequency of the wave, $c$ is the speed of light, $\varepsilon_{o}$, and $\varepsilon_{e}$ are the principal values of the permittivity of tensors. In the conventional dielectric medium, both of the latter quantities take positive values whereas in hyperbolic metamaterials they take opposite values. For instance, they can be either $\varepsilon_{o}>0$ and $\varepsilon_{e}<0$ or $\varepsilon_{o}<0$ and $\varepsilon_{e}>0$. It is easy to produce the structures, having such principal values of permittivity tensors, that were first experimentally illustrated in magnetized plasma [11]. Then, layered hyperbolic metamaterials have been created to demonstrate negative refraction and hyperbolic dispersion [12, 13, 14, 15]. These layered hyperbolic materials have been extensively studied by theoretically to explain the wave dispersion, refraction and propagation in such materials [16, 17, 18]. It is noted that hyperbolic metamaterials are also created in the form of metallic nanowire arrays [19]. It is interesting to point out that waveguides made from the hyperbolic media can act as non-magnetic left-handed media [20]. Very recently, it has been predicted that a periodic arrangement of graphene layers behave as a hyperbolic metamaterials by tuning them from elliptic to hyperbolic dispersion through external voltage [21]. This is because graphene based hyperbolic metamaterial exhibits a strong Purcel effect. In addition to these several other factors, the presence of hyperbolic dispersion and intrinsic nonlinearity makes the hyperbolic metamaterials as very interesting to investigate the problem of electromagnetic wave propagation in such a bulk nonlinear media [22, 23, 24, 25, 26, 27].

To describe the electromagnetic wave propagation in an anisotropic nonlinear hyperbolic media, the following ($2+1$)-dimensional coupled hyperbolic nonlinear Schrödinger (CHNLS) equations have been derived from the Maxwell’s equations under slowly varying envelope approximation [28]. They are written as follows:

[TABLE]

Here, $q_{j}$’s are complex electric field envelopes of a transverse magnetic wave, $z$ is the normalized propagation distance, $x$ and $y$ are the transverse coordinates. In water wave theory, $z$ represents the physical time, $x$ and $y$ denote the spatial coordinates. The dispersion coefficient $\beta$ is equal to $-1$ for a hyperbolic medium whereas it is equal to one for a conventional nonlinear medium. The nonlinearity coefficients $\gamma_{1}=\frac{4\pi\chi^{(3)}}{\varepsilon_{e}}$, and $\gamma_{2}=\frac{4\pi\chi^{(3)}}{\varepsilon_{o}}$ along with $\sigma$ determine the self-focusing and self-defocusing nature of both the ordinary and hyperbolic nonlinear media. This can be decided by the signs of $\gamma_{j}$’s, which essentially depend on the values of $\varepsilon_{e}$ and $\varepsilon_{o}$. As we pointed out earlier, for the conventional media, $\varepsilon_{e}$ and $\varepsilon_{o}$ have the same sign, that is $\varepsilon_{e,o}>0$ or $\varepsilon_{e,o}<0$. However, for a hyperbolic medium, the signs of $\varepsilon_{e}$ and $\varepsilon_{o}$ are opposite. The latter two choices imply that the CHNLS equations or hyperbolic Manakov type system of equations (2) cover the problem of electromagnetic wave propagation in both the hyperbolic media and conventional anisotropic media. This consideration is very much essential in order to compare the solitary wave properties one from the other. However, we point out that the problem of dynamics of electromagnetic solitary waves in the hyperbolic media and their underlying mathematical structures have been focused to a lesser extent [29, 30, 31]. Therefore, the main objective of this paper is to unveil the electromagnetic solitary wave structures associated with the CHNLS Eq. (2) by deriving their analytical forms through the Hirota’s bilinear method and also analyze their propagation properties. We note that if the fields $q_{1}$ and $q_{2}$ do not vary with respect to $y$ or $x$ (or when $\beta=0$), Eq. (2) becomes the well known Manakov type equations [32] describing the evolution of orthogonally polarized light waves in an ordinary anisotropic media or birefringent media. We also note that the real constant $\sigma=\pm 1$ help us to bring out the mixed focusing and defocusing Kerr nonlinearities in a hyperbolic medium in analogy with the mixed coupled NLS system. Due to this, the system (2) with $\sigma=-1$ and $\gamma_{j}>0$, $j=1,2$, or with $\sigma=+1$ and $\gamma_{j}<0$, $j=1,2$, are referred as the mixed CHNLS system. It is interesting to point out that several soliton structures are identified in the nonlinear Schrödinger family of systems [33, 34, 35, 36].

Further, we also analyze a physical situation, in which an electromagnetic field has only one component. In other words, if the field depends on only one independent variable then under this situation, the wave propagation is described by the following scalar hyperbolic nonlinear Schrödinger (SHNLS) equation [37],

[TABLE]

The above SHNLS equation can be considered as a special case of the ($2+1$)-dimensional CHNLS Eq. (2). In Eq. (3), the cubic Kerr nonlinearity coefficient $\gamma$ is defined in physical parameter as $\displaystyle{\gamma=\frac{4\pi\chi^{(3)}}{\varepsilon_{e}}}$, where $\varepsilon_{e}$ is a value of permittivity tensor. For completeness, in the present paper, we also analyze the existence of electromagnetic solitary wave structures in the underlying SHNLS Eq. (3) and reveal their properties. We note that the existence of travelling wave solution for Eq. (3) with $\beta=-1$ was proved in Ref. [38]. We also note that the dynamics of spatial soliton inside a conventional Kerr nonlinear medium [39]. It is also interesting to point out that Eq. (3), with $\beta=-1$, is also appears in the dynamics of gravity waves in deep-water [40].

The rest of the paper is organized as follows. In Section 2, we derive the scalar bright and dark solitary wave solutions and discuss their propagation properties. Then, the different types of vector solitary wave solutions of the CHNLS system (2), such as bright-bright, bright-dark and dark-dark solitary wave solutions, are bring out in Section 3 and explore their properties. Further, we also indicate the stability analysis of solitary wave solutions in Section 4. Finally, we summarize our results in Section 5. In the Appendices, we present the Painlevé singularity structure analysis of the SHNLS Eq. (3) and its coupled version (2) and show that they are non-integrable.

2 Scalar bright and dark electromagnetic solitary waves

We start with the SHNLS Eq. (3) and its bright and dark solitary wave solutions. To derive these solutions, we bilinearize Eq. (3) by adopting the well known Hirota’s bilinear method [41]. We achieve this, through the bilinear transformation, $q(x,y,z)=g(x,y,z)/f(x,y,z)$, where $g(x,y,z)$ is a complex function and $f(x,y,z)$ is a real function. The following bilinear forms of the (2+1)-dimensional SHNLS Eq. (3) are obtained by introducing the latter transformation in it. The resultant bilinear forms are given by

[TABLE]

In the above, $D_{z}$, $D_{x}$ and $D_{y}$ are the standard Hirota’s $D$-operators. According to Hirota [41], these bilinear operators are defined as follows:

[TABLE]

The above Hirota’s notation is very useful to write the bilinear equation in a more compact form. It is noted that the non-integrability nature of the scalar system (3) and its coupled version (2) is established by us and the details are given in Appendix A and B, respectively. Therefore, in the following sections, we refer the wave solutions associated with these systems as solitary wave solutions rather than calling them as ‘soliton solutions’.

2.1 Scalar bright solitary wave

We consider a vanishing boundary condition $q\rightarrow 0$ as $x$ and $y$ go to $\pm\infty$ in order to derive the bright solitary wave solution of Eq. (3). To incorporate this condition in the solution construction process, we fix $\lambda=0$ in the bilinear Eq. (4). Solving the resultant bilinear forms using the truncated series expansion of the unknown functions $g$ and $f$, $g=\epsilon g_{1}$ and $f=1+\epsilon^{2}f_{2}$, where $\epsilon$ is a formal series expansion parameter, we arrive at the fundamental bright solitary wave solution of Eq. (3). It is written as

[TABLE]

From the above solution, one can observe that the structure of bright solitary wave in a hyperbolic medium is determined by the complex constants $\alpha_{1}$ and $k_{j}$, $j=1,2$, and a real constant $\gamma$, apart from a hyperbolic transverse dispersion coefficient $\beta=-1$. We point out that the same number of complex and real constants describe the structure of bright solitary wave in a conventional isotropic nonlinear medium but with the positive transverse dispersion ($\beta=1$) in $y$-direction. We also wish to point out that by repeating the same solution construction procedure with appropriate forms of seed solutions one cannot obtain the two-and general higher-order solitary wave solutions of Eq. (3) with a sufficient number of free parameters because the series expansions do not terminate to a finite order thereby ensuring the non-integrable nature of Eq. (3).

To bring out the propagation properties associated with the bright solitary wave solution (6a)-(6b), we rewrite it in the following hyperbolic form:

[TABLE]

Here, A=\bigg{(}\frac{k_{1R}^{2}+\beta^{3}k_{2R}^{2}}{\gamma}\bigg{)}^{1/2}, $e^{i\theta}=\sqrt{\alpha_{1}/\alpha^{*}_{1}}$, \eta_{1R}=k_{1R}x+\beta k_{2R}y-\big{(}k_{1R}k_{1I}+\beta^{3}k_{2R}k_{2I}\big{)}z, \eta_{1I}=k_{1I}x+\beta k_{2I}y+\frac{1}{2}\big{(}(k_{1R}^{2}-k_{1I}^{2})+\beta^{3}(k_{2R}^{2}-k_{2I}^{2})\big{)}z. The subscripts $R$ and $I$ appearing in the above and in the following represent the real and imaginary parts of that particular complex parameter. The amplitude of the bright solitary wave is \displaystyle{\big{(}\frac{k_{1R}^{2}+\beta^{3}k_{2R}^{2}}{\gamma}\big{)}^{1/2}}. For $\beta=0$, the solution (7) matches with the scalar bright soliton solution of the standard NLS equation [42]. In the scalar Eq. (3), the so-called line or stripe-solitary wave can travel in three different planes, namely ($x-y$)-plane, ($x-z$)-plane, and ($y-z$)-plane, for fixed $z$, $y$, and $x$, respectively. In ($x-y$)-plane , the solitary wave travels with velocity $\beta k_{2R}/k_{1R}$ whereas it propagates in ($x-z$)-plane with velocity $\displaystyle{k_{1I}+\frac{\beta^{3}k_{2R}k_{2I}}{k_{1R}}}$. On the other hand, the solitary wave moves in ($y-z$)-plane with velocity $\displaystyle{\frac{k_{1R}k_{1I}}{\beta k_{2R}}+\beta^{2}k_{2I}}$. From the latter, one can observe that the velocity of solitary wave in both planes depends on the real and imaginary parts of wave numbers $k_{1}$ and $k_{2}$. This implies that one cannot tune the velocity of solitary wave in a particular plane (say ($x-z$)-plane) without affecting the velocity of solitary wave in the other plane (say ($y-z$)-plane). This fact is always true in both the conventional and hyperbolic nonlinear media.

Further, an interesting fact about the hyperbolic nonlinear media is that one can observe the bright solitary wave formation in focusing as well as defocusing Kerr nonlinear regimes. This fact can be realized by carefully looking at the solution (6a)-(6b) of Eq. (3) since it satisfies Eq. (3) in both the focusing ($\gamma>0$) and defocusing ($\gamma<0$) nonlinear regions with $\beta=-1$. As a result, one can observe a non-singular bright solitary wave profile in such a hyperbolic medium by imposing either $k_{1R}^{2}>k_{2R}^{2}$, $\gamma>0$ or $k_{1R}^{2}<k_{2R}^{2}$, $\gamma<0$. In both situations, the quantity $e^{\delta}$ is always positive in Eq. (6b). However, the above fact is not true in the case of conventional nonlinear medium, where a regular bright solitary wave exist only in the focusing nonlinear region. From the solution (6a)-(6b), we realize that it is always singular for $\beta=1$ and $\gamma<0$. It physically implies that the presence of hyperbolic dispersion in the nonlinear hyperbolic metamaterials provides an additional degree of freedom to avoid the formation of singularity in bright solitary wave structures in both the self-focusing and self-defocusing Kerr nonlinear regimes. Such a property is not possible in the conventional Kerr nonlinear medium. A typical propagation of non-singular bright solitary wave in a hyperbolic focusing nonlinear medium is illustrated in Figs. 1a-1c. We have omitted a graphical demonstration of the propagation of the bright solitary wave in a hyperbolic defocusing Kerr nonlinear medium for brevity.

2.2 Scalar dark solitary wave

Next, to derive the dark-solitary wave solution of the system (3): we wish to incorporate the non-vanishing boundary condition $q\rightarrow\tau e^{i\theta}$, $\theta=l_{1}x+\beta l_{2}y-\omega z$, as $x$ and $y$ go to $\pm\infty$ in the solution derivation procedure by setting $\lambda\neq 0$, which can be determined later, in the bilinear forms (4). Then, solving Eq. (4) along with the truncated series expansions, $g=g_{0}(1+\epsilon g_{1})$, $f=1+\epsilon f_{1}$, for the unknown functions, $g$ and $f$, we obtained the fundamental dark solitary wave solution of Eq. (3) as

[TABLE]

where $\theta=l_{1}x+\beta l_{2}y-\frac{1}{2}(l_{1}^{2}+\beta^{3}l_{2}^{2})z-\lambda z$, $\lambda=-\gamma\tau^{2}$, along with a constraint equation for $\tau$:

[TABLE]

In the above, the wave parameters, $\tau$, $k_{j}$, $l_{j}$, $j=1,2$, $\Omega_{1}$ and $\xi_{1}^{(0)}$ are real constants. The amplitude of the background plane wave $\tau$ can be a complex constant but in the present case we treat it as a real one. By rewriting the solution (8a), we obtain the following form of the fundamental dark solitary wave solution. It reads as

[TABLE]

Here, an interesting fact associated with dark-solitary wave solution (10) is that it appears also in the focusing nonlinear region. This possibility is allowed in a hyperbolic nonlinear medium. The intensity of the scalar dark solitary wave can be calculated from Eq. (10) as [39],

[TABLE]

where

[TABLE]

The darkness of the solitary wave (10) can be controlled by tuning the value of $X$ appropriately using the wave parameters involve in it. If $X=1$, we get a dark solitary wave profile whereas a gray solitary wave profile appears for $X<1$. In the present scalar hyperbolic medium as well as in the conventional nonlinear medium, the profile of dark-solitary wave is always regular. The propagation of a such non-singular dark solitary wave, for $\beta=-1$, in three different planes is illustrated in Fig. 2. For instance, the dark solitary wave propagation in ($x-y$) plane is demonstrated in Fig. 2a for $z=0$, where it moves with velocity $\displaystyle{\beta k_{2}/k_{1}}$. Then, the dark-solitary wave propagates in ($x-z$)-plane, with velocity $\Omega_{1}/k_{1}$, is depicted in Fig. 2b. A typical dark-solitary wave dynamics in ($y-z$)-plane is illustrated in Fig. 2c, where it propagates with velocity $\Omega_{1}/\beta k_{2}$. From the above, we observe that one can tune the velocity of dark solitary wave in ($x-z$)-plane through wave numbers without affecting the velocity in ($y-z$)-plane for fixed value of $\Omega_{1}$. The hyperbolic dispersion $\beta$ only influences the grayness of the dark solitary wave. Very recently, the grayness of the dark-solion in the generalized long-wave short-wave resonance interaction system has been studied in [43] by tuning the system parameter. It is noted that one can study the interaction properties of these solitary wave solutions numerically by considering our solutions as initial conditions [44]. We point out that the dark-soliton solution of the defocusing NLS equation can be recovered by imposing $\beta=0$ in Eq. (10).

3 Vector solitary waves of the CHNLS system (2)

Now, we intend to study the propagation properties of the various vector electromagnetic solitary waves in a hyperbolic nonlinear medium in which the electric field has two components. To bring out these properties, we now consider a coupled version of the SHNLS Eq. (3). An important feature of the CHNLS system (2) is that it admits a rich structure of vector solitary wave solutions [45], such as bright-bright, bright-dark or dark-bright, dark-dark type solitary wave solutions depending on the choice of boundary conditions. To construct these solutions, we again consider the bilinearization method. The CHNLS Eq. (2) can be bilinearized through the transformation, $q_{j}(x,y,z)=g^{(j)}(x,y,z)/f(x,y,z)$, $j=1,2$, where $g^{(j)}(x,y,z)$’s are in general complex functions and $f(x,y,z)$ is a real function. By substituting the latter dependent variable transformation in Eq. (2), we arrive at the following forms of bilinear equations. They are written as

[TABLE]

where $D_{1}\equiv(iD_{z}+\frac{1}{2}(D_{x}^{2}+\beta D_{y}^{2})-\lambda)$, and $D_{2}\equiv(D_{x}^{2}+\beta D_{y}^{2}-2\lambda)$. We note here that the bilinearization of the CHNLS system (2) enforces us to fix the nonlinear coefficients $\gamma_{1}$ and $\gamma_{2}$ as $\gamma_{1}=\gamma_{2}=\gamma$. By solving the above bilinear forms (11) along with different forms of series expansions for the unknown functions, $g^{(j)}$, $j=1,2$, and $f$, we are able to derive the various types of vector solitary wave solutions of the CHNLS system (2) as given below.

3.1 Bright-bright solitary wave

To obtain the fundamental vector bright solitary wave solution of the hyperbolic coupled NLS system (2), we assume the truncated series expansions, $g^{(j)}(x,y,z)=\epsilon g_{1}^{(j)}$, $j=1,2$, $f(x,y,z)=1+\epsilon^{2}f_{2}$, and $\lambda=0$. By solving the system of equations, which arise while considering the latter series expansions in Eq. (11), we obtain the bright-bright vector solitary wave solution of Eq. (2) as

[TABLE]

In the above solution, the parameters, $\alpha_{1}^{(j)}$, and $k_{j}$, $j=1,2$, are complex and the system parameters, $\gamma$, $\beta$, and $\sigma$ are real constants. If the complex constant $\alpha_{1}^{(1)}=\alpha_{1}^{(2)}$, then the component $q_{1}=q_{2}$. Correspondingly, the solution (12) becomes a solution of the SHNLS (3) with a modified form of the nonlinearity. We note that to derive the solution (12), the seed solutions, $g^{(j)}_{1}=\alpha_{1}^{(j)}e^{\eta_{1}}$, have been considered during the solution construction process. Besides this, we have also tried to construct the recently identified nondegenerate vector bright soliton/solitary wave solutions [46, 47, 48, 49] for Eq. (2) by assuming two distinct starting solutions, $g^{(j)}_{1}=\alpha_{1}^{(j)}e^{\eta_{j}}$, $\eta_{j}=k_{j}x+p_{j}y-\omega_{j}z$, $j=1,2$. However, the series expansions for the unknown functions $g^{(j)}$, $j=1,2$, and $f$ do not converge to a finite order so we are unable to find such nondegenerate vector bright solitary wave solution for Eq. (2) through the Hirota bilinear method. We remark that one cannot find the higher-order solitary wave solutions of kind (12), with a sufficient number of wave parameters, since the series expansions that we consider for the unknown functions $g^{(j)}$, $j=1,2$, and $f$ do not terminate to a finite order. This is due to the non-integrable nature of the CHNLS system (2) as we proved it in Appendix B.

In order to reveal the propagation properties associated with the vector bright solitary wave, we rewrite the solution (12) as follows:

[TABLE]

where

[TABLE]

The amplitudes of the bright solitary wave in both the components $q_{1}$ and $q_{2}$ are given by \displaystyle{\hat{A}_{1}\big{(}\frac{k_{1R}^{2}+\beta^{3}k_{2R}^{2}}{\gamma}\big{)}^{1/2}}, and \displaystyle{\hat{A}_{2}\big{(}\frac{k_{1R}^{2}+\beta^{3}k_{2R}^{2}}{\gamma}\big{)}^{1/2}}, respectively. Here, ${A}=\begin{pmatrix}\hat{A}_{1}&\hat{A}_{2}\end{pmatrix}^{T}$ is a unit polarization vector, which obeys $A^{\dagger}BA=1$, $B=\begin{pmatrix}1&0\\ 0&\sigma\end{pmatrix}$. The complex constants $\alpha_{1}^{(j)}$’s play a crucial role in fixing the intensity of the vector bright line-solitary waves. In the vector Eq. (2), the line solitary wave can propagates in three different planes, like in the case of scalar Eq. (3), namely ($x-y$)-plane, ($x-z$)-plane, and ($y-z$)-plane, for fixed $z$, $y$, and $x$, respectively. The velocity of vector solitary wave in ($x-z$)-plane is $\displaystyle{k_{1I}+\frac{\beta^{3}k_{2R}k_{2I}}{k_{1R}}}$, while the bright solitary wave travels in ($y-z$)-plane with velocity $\displaystyle{\frac{k_{1R}k_{1I}}{\beta k_{2R}}+\beta^{2}k_{2I}}$. We displayed these propagation scenarios of vector bright solitary waves in Fig. 3. We note that the constant, $\beta$, which explicitly appears on the velocity of solitary wave.

Further, the singularity nature of the vector bright solitary wave (13) can be understood from the quantity $e^{R_{1}}$, which mainly depends on the value of $\alpha_{1}^{(j)}$, $k_{jR}$, $j=1,2$, $\beta$, $\gamma$, and $\sigma$. In the hyperbolic nonlinear media, a regular bright solitary wave profile appears in different nonlinear regions by fixing the conditions as follows: (i) Focusing nonlinear case ($\gamma,\sigma>0$): $k_{1R}^{2}>k_{2R}^{2}$, (ii) Defocusing nonlinear case ($\gamma,\sigma<0$): $k_{1R}^{2}<k_{2R}^{2}$, $\lvert\alpha_{1}^{(1)}\rvert>\lvert\alpha_{1}^{(2)}\rvert$, (iii) Mixed nonlinear case ($\gamma>0$, $\sigma=-1$ or $\gamma<0$, $\sigma=1$): $k_{1R}^{2}>k_{2R}^{2}$, $\lvert\alpha_{1}^{(1)}\rvert>\lvert\alpha_{1}^{(2)}\rvert$ (or $k_{1R}^{2}<k_{2R}^{2}$). However, this is not true in the case of conventional anisotropic nonlinear medium, where a non-singular vector bright solitary wave exists in the focusing nonlinear regime whereas a bright solitary wave becomes singular in the defocusing nonlinear regime. Then, it is possible to observe a regular bright solitary wave profile in the conventional medium with mixed nonlinearity. For instance, the quantity $e^{R_{1}}$ becomes positive definite for $\gamma>0$ and $\sigma=-1$ (or $\gamma<0$ and $\sigma=-1$) when $\lvert\alpha_{1}^{(1)}\rvert>\lvert\alpha_{1}^{(2)}\rvert$ (or $\lvert\alpha_{1}^{(1)}\rvert<\lvert\alpha_{1}^{(2)}\rvert$). From the above, one can clearly confirm that by setting the wave numbers $k_{1}$ and $k_{2}$, which determine the propagation behaviors of vector bright solitary waves along $x$ and $y$ directions, in such a way one can avoid the formation of singularity in all the nonlinear regions of the hyperbolic media.

3.2 Bright-dark vector solitary wave

To investigate the role of dark solitary wave on the bright solitary wave component, we intend to construct the fundamental mixed (bright-dark) vector solitary wave solution of the system (2). For this purpose, we consider the series expansions as $g^{(1)}=\epsilon g_{1}^{(1)}$, $g^{(2)}=g_{0}(1+\epsilon^{2}g_{2}^{(2)})$, and $f=1+\epsilon^{2}f_{2}$ along with the bilinear Eq. (11) and $\lambda\neq 0$. In this case, the boundary conditions are assumed as $q_{1}\rightarrow 0$, and $q_{2}\rightarrow\tau e^{i\gamma\sigma\tau^{2}z}$, as $x$ and $y$ go to $\pm\infty$. By following the same solution construction procedure, the fundamental bright-dark solitary wave solution is obtained as

[TABLE]

where

[TABLE]

By rewriting the above solution, we arrive at the following compact form of the mixed vector solitary wave solution of Eq. (2). It reads as

[TABLE]

where

[TABLE]

The above mixed vector solution (15a)-(15b) is characterized by three arbitrary complex constants, $\alpha_{1}$, $k_{1}$ and $k_{2}$, and a real constant $\tau$, apart from three system parameters, $\gamma$, $\sigma$, and $\beta$. The amplitude of the bright solitary wave is $\displaystyle{A(\frac{k_{1R}^{2}+\beta^{3}k_{2R}^{2}}{\gamma})^{1/2}}$. Here, $A$ is related to the polarization vector of the bright solitary wave and is independent on the complex parameter $\alpha_{1}$, which appears as the complex phase $e^{i\theta}$. The bright and dark mixed vector solitary waves travel with velocity $\displaystyle{k_{1I}+\frac{\beta^{3}k_{2R}k_{2I}}{k_{1R}}}$ in ($x-z$)-plane, while they propagate in ($y-z$)-plane with velocity $\displaystyle{\frac{k_{1R}k_{1I}}{\beta k_{2R}}+\beta^{2}k_{2I}}$. The propagation of bright and dark-solitary waves in three different planes is displayed in Fig. 4.

The intensity of dark-solitary wave component is calculated as \lvert q_{2}\rvert^{2}=\tau^{2}\big{(}1-\cos^{2}\phi~{}\mbox{sech}^{2}(\eta_{1R}+\frac{R_{2}}{2})\big{)}. From the latter, one can observe that the value of $\phi$ determines the degree of grayness of the dark-solitary wave in $q_{2}$ component. Besides this, an interesting feature of the mixed vector solitary wave solution is that the presence of dark-solitary wave in $q_{2}$ component influences the structure of bright solitary wave in the other $q_{1}$ component. For instance, in the focusing and defocusing hyperbolic nonlinear media, the intensities of bright and dark solitary waves are increasing with respect to the background wave parameter $\tau$ as it is illustrated in Fig. 5. In addition to this, in the mixed nonlinear case, we observe that the intensity of the bright solitary wave is decreasing while increasing the value of $\tau$. Such a possibility is illustrated in Fig. 6 for different values of $\tau$. From the solution (15a)-(15b), it is obvious that the hyperbolic nonlinear media allows to observe the mixed bright-dark solitary waves.

3.3 Dark-dark vector solitary wave

Finally, to show the existence of the fundamental vector dark solitary wave in a hyperbolic medium, we derive the corresponding solution by considering the boundary condition: $\lvert q_{j}\rvert^{2}\rightarrow\tau_{j}^{2}$, $j=1,2$, as $x$ and $y$ go to $\pm\infty$. We derive this solitary wave solution by solving the bilinear Eq. (11) along with the following truncated series expansions, $g^{(1)}=g_{0}^{(1)}(1+\epsilon g_{1}^{(1)})$, $g^{(2)}=g_{0}^{(2)}(1+\epsilon g_{1}^{(2)})$, and $f=1+\epsilon f_{1}$. This action yields the fundamental vector dark-solitary wave solution of Eq. (2). It reads as

[TABLE]

Here, $\theta_{1}=l_{1}x+\beta l_{2}y-\omega_{1}z$, $\theta_{2}=l_{3}x+\beta l_{4}y-\omega_{2}z$, $\omega_{1}=\frac{1}{2}(l_{1}^{2}+\beta^{3}l_{2}^{2})+\lambda$, $\omega_{2}=\frac{1}{2}(l_{3}^{2}+\beta^{3}l_{4}^{2})+\lambda$, $\lambda=-\gamma(\tau_{1}^{2}+\sigma\tau_{2}^{2})$, and the constraint equation is given by

[TABLE]

We rewrite the above dark-dark solitary wave solution (16) as

[TABLE]

The above vector dark solitary wave is characterized by the real parameters $k_{j}$, $l_{j}$, $\tau_{j}$, $j=1,2$, $l_{3}$, $l_{4}$, $\gamma$, $\beta$, and $\sigma$. The grayness of the dark-solitary wave in two components are determined by calculating intensity from the expressions (17a) and (17b). By doing so, we find

[TABLE]

where

[TABLE]

One can get a complete dark-solitary wave profile in both components, if the above quantities $X_{1}$ and $X_{2}$ are equal to one. On the other hand, we get a gray-solitary wave profiles if $X_{1,2}<1$. The various propagation of vector dark solitary wave, for $\beta=-1$, is illustrated in Fig. (7). For example, the propagation of vector dark solitary wave in ($x-y$)-plane is demonstrated in Fig. (7)a1-a2, where it propagates with velocity $\displaystyle{\beta k_{2}/k_{1}}$. Then, the dark-solitary wave propagate in ($x-z$) plane for $y=0$, with velocity $\Omega_{1}/k_{1}$, is depicted in Fig. (7)b1-b2. A typical dark-solitary wave propagation in ($y-z$)-plane is displayed in Fig. (7)c1-c2, where it travels with velocity $\Omega_{1}/\beta k_{2}$. One can control the velocity of dark solitary wave, in both modes, in ($x-z$)-plane by tuning the values of wave numbers $k_{1}$ and $k_{2}$ without affecting the velocity in ($y-z$)-plane for fixed $\Omega_{1}$. The hyperbolic dispersion $\beta=-1$ is very useful to control the grayness of the vector dark solitary waves. We note that the vector dark solitary wave solution (16) becomes a solution of the scalar Eq. (3) if the restriction, $l_{3}=l_{1}$ and $l_{4}=l_{2}$, is imposed on it. Under this condition, $Z_{2}=Z_{1}$, the solution (16) turns out to be q_{j}=\frac{\tau_{j}}{2}e^{i\theta}\bigg{(}1+Z_{1}-(1-Z_{1})\tanh\frac{\xi_{1}}{2}\bigg{)}, where $\theta_{1}=\theta_{2}=\theta=l_{1}x+l_{2}y-\omega_{1}z$, $\omega_{1}=\frac{1}{2}(l_{1}^{2}+\beta l_{2}^{2})+\lambda$, and $\lambda=-\gamma(\tau_{1}^{2}+\sigma\tau_{2}^{2})$.

4 Transverse instability analysis

We examine the stability of obtained scalar solitary wave solutions by treating the transverse dispersion term $\beta q_{yy}$ as perturbation to the $1$-dimensional scalar NLS soliton. To start with, we consider a very special case of the solution (7) with $\beta=0$ and $z=0$ ($1$-dimensional stationary scalar NLS soliton). In other words, in order to examine the stability of

[TABLE]

which satisfies

[TABLE]

with respect to the transverse perturbations, we consider

[TABLE]

In the above, $\lambda$ is a complex constant which is related to the spectral problem (Eq. (23)) given below, $k$ is the transverse wavenumber (which is real), and $v(x)$ and $w(x)<<1$ are small transverse perturbations. Here, we note that Eq. (21) is obtained by substituting Eq. (20), after taking $z\rightarrow 2z^{\prime}$ and choosing $\gamma=1$, in Eq. (3). Then, substituting (22) in (3), with again $z=2z^{\prime}$ and $\gamma=1$, and linearize the resulting equation for $v(x)$ and $w(x)$, we get the following eigenvalue problem:

[TABLE]

where $\mathcal{L}_{0}=\partial_{x}^{2}-1+2q_{0}^{2}-\beta k^{2}$, $\mathcal{L}_{1}=\partial_{x}^{2}-1+6q_{0}^{2}-\beta k^{2}$, and $q_{0}=\mbox{sech}x$. The above eigenvalue problem describes the stability of $1$-dimensional NLS soliton. One cannot find the exact analytical solutions for the linear eigenvalue problem (23). However, this eigenvalue problem was extensively analysed in the literature both by analytically and numerically and reported the approximate growth rate for unstable eigenvalue $\lambda(k)$ and a threshold condition for the onset of the transverse instability [50, 51, 52, 53, 54, 40, 55, 56, 57, 58, 59, 60, 61, 62]. One may refer the paper (Ref. [62]) by Deconinck et al. for a complete analysis on the transverse instability of NLS soliton and also the previously achieved important results along this line are summarized there. Therefore, in order to avoid the repetition, we provide here only the expressions for the growth rates $\lambda(k)$, which was first obtained by Zakharov-Rubenchik [50], for the long ($k\rightarrow 0$) and short ($k\rightarrow k_{c}$, here $k_{c}$ refers a critical value) transverse waves. Based on these expressions and their corresponding plots, we conclude the stability of our obtained solitary wave solutions. To achieve the expressions for approximate growth rates the asymptotic expansion method was adopted. This process yields the growth rate for the elliptic NLS case ($\beta=1$) is given by

[TABLE]

The above growth rate corresponds to the long transverse waves. However, for the short transverse waves one may refer [62], where a cut-off condition for the onset of instability is given. We draw Eq. (24) by varying the transverse wavenumber $k$ in Fig. 8(a), where the growth rate reaches its maximum value for $k\approx 0.6$.

For the hyperbolic case ($\beta=-1$), the approximate expressions for the four branches of $\lambda(k)$ have been obtained in [50, 62] by considering and long transverse waves as perturbations. Out of which two of them are imaginary whereas the other two is a real negative and a real positive value of $\lambda$. They are given as follows:

[TABLE]

and

[TABLE]

In the above, $\lambda_{3,4H}$ denotes a pair of imaginary eigenvalues, $\lambda_{1H}$ and $\lambda_{2H}$ refer to the unstable (positive) and the stable (negative) eigenvalues, respectively, and $H$ represents the hyperbolic case. We draw the growth rate $\lambda(k)$ corresponds to the unstable eigenvalue in Fig. 8(b) for certain values of $k$. However, there is a cut-off in instability curve after $k=1$. One can understand the reason behind this cut-off from Ref. [62]. Similarly, we also draw the imaginary part of growth rate $\lambda_{3H}$ in Fig. 8(c) for certain values of $k$. From these figures (and also from their corresponding growth rates expressions (25) and (26)), and from Eq. (22), we infer that the scalar NLS soliton $q_{0}=\mbox{sech}~{}x$ is unstable if eigenvalue $\lambda$ is real positive ($\lambda_{1H}$). This is due to the presence of exponential factor $e^{\lambda z}$ in Eq. (22). Because of the latter factor the amplification will occur in the corresponding perturbation profile. On the other hand, if $\lambda$ is complex the effect of instability is determined by the real part of $\lambda_{3,4H}$ leading to the amplification of perturbation along with oscillatory instability which is due to the imaginary part of $\lambda_{3,4H}$. If the eigenvalue is a real negative value, then the NLS soliton is stable against the transverse perturbation. In a similar way, one can do the stability analysis for the rest of our solitary wave solutions. Thus, we conclude that the scalar and coupled solitary waves presented in this paper are also prone to perturbations in transverse directions leading to the transverse instability.

.

5 Conclusion

In this paper, we have considered the scalar and vector hyperbolic nonlinear Schrödinger equations to study the propagation of electromagnetic solitary waves in nonlinear hyperbolic metamaterials. These equations are shown to be non-integrable since they fail to pass the Painlevé test. Then, we have investigated the propagation properties of both scalar and vector electromagnetic solitary waves by deriving their exact analytical forms through the Hirota bilinear method. In the scalar case, we have obtained the bright and dark solitary waves solutions and revealed their interesting propagation properties. On the other hand, in the vector case, we have presented the various types of vector solitary wave solutions, such as bright-bright, bright-dark and dark-dark solitary wave solutions of the CHNLS system (2). Out of which, in the mixed vector solitary wave solution, the existence of dark solitary wave in one component strongly influences the structure of bright solitary wave in another component. Another interesting main result of our study is that the presence of negative/hyperbolic transverse dispersion in hyperbolic nonlinear media prevents the formation of wave collapsing phenomenon by avoiding the singularity formation. We observe that this is essentially due to the presence of two distinct transverse wavenumbers in the obtained solitary wave solutions. This striking property makes our solitary wave solutions as special. Then, the hyperbolic transverse dispersion also allows to observe bright and dark solitary waves in both focusing and defocusing Kerr nonlinear regions. We believe that the results reported in this paper will be useful for studying the electromagnetic wave propagation in the bulk nonlinear hyperbolic media.

Declaration of competing interest

The authors declare that they have no conflict of interest.

Acknowledgement

The works of Mokhtar Kirane, and Stalin Seenimuthu, are supported by the Khalifa University of Science and Technology (KUST), Abu-Dhabi, UAE, under the Project Grant No. 8474000355.

Appendix A Painlevé singularity structure analysis of SHNLS system (3)

In this section, we intend to analyze the integrability of the system (3) through the WTC (Weiss, Tabor and Carnevale) algorithm [63]. To investigate integrability of the scalar hyperbolic NLS system (3), we rename the nonlinear Schrödinger field $q(x,y,z)$ its complex conjugate field $q^{*}(x,y,z)$ as $a(x,y,z)$ and $b(x,y,z)$, respectively. Then Eq. (3) and its complex conjugate equation can be rewritten as follows:

[TABLE]

The Painlevé singularity structure analysis for the above coupled equations can be carried out by expanding their solutions in a standard way in the neighborhood of a non-characteristic manifold $\phi(x,y,z)=0$, with $\phi_{x},\phi_{y},\phi_{z}\neq 0$. The Laurent series expansion for the solutions $a(x,y,z)$ and $b(x,y,z)$ are

[TABLE]

where $\nu_{k}$, $k=1,2$, are integers, which are to be determined later. In the Painlevé test, the first step is to analyze the leading order behaviour of the solutions of Eqs. (27a) and (27b). To do this, we restrict the series expansion (28a) and (28b) as $a=a_{0}\phi^{\nu_{1}}$ and $b=b_{0}\phi^{\nu_{2}}$ by assuming $j=0$. The values of $\nu_{1}$ and $\nu_{2}$ are identified by substituting the latter leading order solutions in Eqs. (27a) and (27b) and balancing the most negative dominant terms. By doing so, we find $\nu_{k}=-1$, $k=1,2$, and the following condition

[TABLE]

From the above equation, one can observe that a single equation is obtained for the two unknowns $a_{0}$ and $b_{0}$, so one is arbitrary. The next step is to obtain the resonance values ($j$), at which the arbitrary functions appear in the Laurent series (28a)-(28b). This can be determined by substituting Eqs. (28a)-(28b) in Eqs. (27a) and (27b) and collect the coefficients of $\phi^{j-3}$, which gives

[TABLE]

From the above equation, one can get the quartic polynomial in $j$ and solving such a polynomial leads to the resonance values: $j=-1,0,3$ and $4$. Here, the resonance value $j=-1$ represents the arbitrariness of the singular manifold $\phi(x,y,z)$ and $j=0$ denotes the arbitrariness of either $a_{0}$ or $b_{0}$, which can be confirmed from Eq. (29).

Next, we check the existence of sufficient number of arbitrary functions at the resonance values $j=3$ and $j=4$ without introduction of any movable singular manifold. To verify this, we first collect the coefficients of $\phi^{-2}$, which corresponds to the resonance value $j=1$. By doing so, we get the following equations which are expressed in matrix form as

[TABLE]

where $A_{1}=2\gamma a_{0}b_{0}$, $A_{2}=\gamma b_{0}^{2}$, $B_{1}=\gamma a_{0}^{2}$, c_{1}=ia_{0}\phi_{z}+\frac{1}{2}\big{(}2a_{0,x}\phi_{x}+a_{0}\phi_{xx}\big{)}+\frac{\beta}{2}\big{(}2a_{0,y}\phi_{y}+a_{0}\phi_{yy}\big{)}, and c_{2}=-ib_{0}\phi_{z}+\frac{1}{2}\big{(}2b_{0,x}\phi_{x}+b_{0}\phi_{xx}\big{)}+\frac{\beta}{2}\big{(}2b_{0,y}\phi_{y}+b_{0}\phi_{yy}\big{)}. By solving, Eq. (30), we obtain the expressions for $a_{1}$ and $b_{1}$. which are given below:

[TABLE]

From the above, we confirm that there is no arbitrary function at this order. Similarly, we also confirm no arbitrary function enters into the Laurent series (28a)-(28b) at the resonance value $j=2$. It can be verified by collecting the coefficients of $\phi^{-1}$. This action gives the matrix equation

[TABLE]

where d_{1}=-ia_{0,z}-\frac{1}{2}a_{0,xx}-\frac{\beta}{2}a_{0,yy}-\gamma\big{(}a_{1}^{2}b_{0}+2a_{0}a_{1}b_{1}\big{)}, d_{2}=ib_{0,z}-\frac{1}{2}b_{0,xx}-\frac{\beta}{2}b_{0,yy}-\gamma\big{(}a_{0}b_{1}^{2}+2b_{0}a_{1}b_{1}\big{)}, and other elements are already defined above. By solving the above equation we obtain

[TABLE]

This shows that there is no arbitrary function at this stage. However, we observe the breaking of arbitrariness at the resonance value $j=3$. To ensure this, collect the coefficients of $\phi^{0}$. This yields the following equations, which are defined in matrix form:

[TABLE]

Here, c_{1}^{\prime}=\frac{1}{\gamma a_{0}}\big{(}-ia_{1,z}-ia_{2}\phi_{z}-\frac{1}{2}(a_{1,xx}+2a_{2,x}\phi_{x}+a_{2}\phi_{xx})-\frac{\beta}{2}(a_{1,yy}+2a_{2,y}\phi_{y}+a_{2}\phi_{yy})-\gamma(2a_{1}a_{2}b_{0}+2a_{0}a_{2}b_{1}+a_{1}^{2}b_{1}+2a_{0}a_{1}b_{2})\big{)}, and c_{2}^{\prime}=\frac{1}{\gamma b_{0}}\big{(}ib_{1,z}+ib_{2}\phi_{z}-\frac{1}{2}(b_{1,xx}+2b_{2,x}\phi_{x}+b_{2}\phi_{xx})-\frac{\beta}{2}(b_{1,yy}+2b_{2,y}\phi_{y}+b_{2}\phi_{yy})-\gamma(2b_{1}b_{2}a_{0}+2b_{0}b_{2}a_{1}+b_{1}^{2}a_{1}+2b_{0}b_{1}a_{2})\big{)}. From the above one can easily identify that the elements $c_{1}^{\prime}$ and $c_{2}^{\prime}$ are distinct and they become identical if and only if $\beta=0$, which corresponds to the standard NLS equation. This clearly suggest that the scalar hyperbolic NLS equation (3) fails to satisfy the Painlevé test at this order.

Now, we collect the coefficients of $\phi^{1}$ to check the arbitrariness of the functions $a_{4}(x,y,z)$ and $b_{4}(x,y,z)$ at the resonance value $j=4$. By doing so, we obtain

[TABLE]

where d_{1}^{\prime}=\frac{1}{a_{0}}\big{(}ia_{2,z}+2ia_{3}\phi_{z}+\frac{1}{2}(a_{2,xx}+4a_{3,x}\phi_{x}+2a_{3}\phi_{xx})+\frac{\beta}{2}(a_{2,yy}+4a_{3,y}\phi_{y}+2a_{3}\phi_{yy})+\gamma(2a_{1}a_{3}b_{0}+2a_{0}a_{3}b_{1}+2a_{0}a_{1}b_{3}+b_{0}a_{2}^{2}+2a_{0}a_{2}b_{2}+a_{1}^{2}b_{2}+2a_{1}a_{2}b_{1})\big{)}, and d_{2}^{\prime}=\frac{1}{b_{0}}\big{(}ib_{2,z}+2ib_{3}\phi_{z}-\frac{1}{2}(b_{2,xx}+4b_{3,x}\phi_{x}+2b_{3}\phi_{xx})-\frac{\beta}{2}(b_{2,yy}+4b_{3,y}\phi_{y}+2b_{3}\phi_{yy})-\gamma(2b_{1}b_{3}a_{0}+2b_{0}b_{3}a_{1}+2b_{0}a_{3}b_{1}+a_{0}b_{2}^{2}+2b_{0}a_{2}b_{2}+a_{2}b_{1}^{2}+2a_{1}b_{1}b_{2})\big{)}. This shows that the elements $d_{1}^{\prime}$ and $d_{2}^{\prime}$ are non-identical except for $\beta=0$. Thus, we confirm that the SHNLS system (3) fails to pass the Painlevé analysis due to the non-existence of sufficient number of arbitrary functions in the Laurent series. Therefore, the SHNLS equation (3) is not integrable in the Painlevé sense.

Appendix B Singularity structure analysis of CHNLS system (2)

Next, we move on proving the non-integrable nature of the coupled hyperbolic nonlinear Schrödinger system (2). However, we omit the most of the mathematical details related to this singularity analysis due to the cumbersome nature of equations. To begin with, we rewrite the coupled system (2) in terms of four real functions $P(x,y,z)$, $Q(x,y,z)$, $R(x,y,z)$, and $S(x,y,z)$, which are defined by $q_{1}=P+iQ$, and $q_{2}=R+iS$. Consequently, we get the following four coupled equations:

[TABLE]

As we have described above for the SHNLS equation, we have to perform the singularity structure analysis for the above set of coupled equations. To start with, we consider the leading orders of the form

[TABLE]

where $\nu_{j}$, $j=1,2,3,4$ are integers to be determined later. Substituting the above leading order forms in Eqs. (36a)-(36d) and equating the most dominant terms we obtain the following choice:

[TABLE]

and $\nu_{1}=\nu_{2}=\nu_{3}=\nu_{4}=-1$. The remaining set of equations resulting from Eqs. (36a)-(36d) are identical to the above set of equations. At this stage, one can say that out of four functions, $P_{0}$, $Q_{0}$, $R_{0}$, and $S_{0}$, two of them are arbitrary. However, the obtained resonance values, $j=-1,0,0,0,3,3,3,4$, show that $Q_{0}$ (or $S_{0}$) is arbitrary besides $P_{0}$ and $R_{0}$. Thus, the obtained Eqs. (38a) and (38b) reduce to a single equation, \gamma\big{(}P_{0}^{2}+Q_{0}^{2}+\sigma(R_{0}^{2}+S_{0}^{2})\big{)}=-(\phi_{x}^{2}+\beta\phi_{y}^{2}) along with a condition $\gamma_{1}=\gamma_{2}=\gamma$. Then, we found that the non-existence of sufficient number of arbitrary functions at the resonance values $j=3$ and $j=4$. This indicates that the CHNLS system (2) is also non-integrable in the Painlevé sense.