## Abstract

The discovery of a new type of soliton occurring in periodic systems is reported. This type of nonlinear excitation exists at a Dirac point of a photonic band structure, and features an oscillating tail that damps algebraically. Solitons in periodic systems are localized states traditionally supported by photonic bandgaps. Here, it is found that besides photonic bandgaps, a Dirac point in the band structure of triangular optical lattices can also sustain solitons. Apart from their theoretical impact within the soliton theory, they have many potential uses because such solitons are possible in both Kerr material and photorefractive crystals that possess self-focusing and self-defocusing nonlinearities. The findings enrich the soliton family and provide information for studies of nonlinear waves in many branches of physics.

© 2017 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

It is well-known that the photonic bandgap of a periodic system can support linear localized modes at a defect. It is also well-known that, if the defect is induced by nonlinearity of the medium, the corresponding self-localized mode is called gap soliton. Confinement of waves within a finite area is the basis of all information processing [1–4]. Traditionally, such wave trapping is achieved by cavities and waveguides that rely on total internal reflection, or a photonic bandgap, to suppress radiation losses [5, 6]. Cavities and waveguides can be formed by a high-index core surrounded by a cladding with a lower refractive index so that total internal reflection can take place. Alternatively, cavities and waveguides can be formed in discrete, or periodic, systems as defects. Due to the presence of allowed bands and forbidden gaps, radiation losses of a wave packet accommodated by the defect are suppressed by any photonic bandgap of the systems. There are various other pathways towards localization that uniquely involve nonlinearity. Indeed, localized modes due to nonlinearity are commonly called solitons [7, 8]. In homogeneous media, nonlinearity raises the refractive index of the media so that light creates its own high-index core. In this sense, light of a conventional soliton is still essentially guided by total internal reflection. In periodic systems, nonlinearity creates its own defect by changing its onsite refractive index, so that light is trapped to the defect by photonic bandgaps of the periodic lattices, forming a gap soliton in this case. Solitons, as nonlinear self-trapped wave packets, have been extensively studied in many branches of physics including optics [9–11]. Among them, algebraic solitons with power-law asymptotic behavior arise mathematically from a particular limit of hyperbolic solitons [12].

In parallel with the above advances, research on graphene has made great progress [13]. The electronic band structure of graphene contains Dirac cones at the six corners of the hexagonal Brillouin zone. The Dirac point is a conical singularity of a band structure where the associated energy-wavenumber relation resembles the two-dimensional massless Dirac equation −*iv*(*σ _{x}*∂

*+*

_{x}*σ*∂

_{y}*)Ψ = (*

_{y}*ω*−

*ω*)Ψ for relativistic electrons in a vacuum, where

_{D}*v*is the group velocity,

*ω*/2π is the Dirac frequency,

_{D}*σ*and

_{x}*σ*are Pauli matrices, and |Ψ|

_{y}^{2}is the probability of finding the electrons in space. The band structure of a photonic crystal formed by a triangular lattice also possesses Dirac cones at the corners of the Brillouin zones [14]. At these high-symmetry points Maxwell’s equations can be replaced by the Dirac equation, with Ψ being the wave functions of two degenerate Bloch states. Due to similarity between band structures of a solid and of a photonic crystal, recent years have seen an increased interest in simulations of relativistic quantum effects using photonic structures [15]. Wave behavior at the Dirac frequency has been studied extensively in photonic crystals [16–20], and localized modes have also been found recently at the Dirac frequency [21, 22]. It has been shown that the Dirac point in the band structures of these lattices can take the role of a bandgap to form localized modes at a defect. A soliton in a photonic crystal is, essentially, a localized mode at a nonlinearly-induced defect. Since the Dirac point can suppress radiation loss like a band gap, it is natural to imagine the existence of Dirac-point soliton, in parallel to gap soliton. As long as a defect is formed by nonlinearity, it is understandable that a self-localized mode can be supported around the defect by the Dirac point.

In this work we report the discovery of a new type of soliton, the Dirac-point soliton that exists at the Dirac point in a nonlinear triangular photonic lattice. We begin by studying the case of a linear defect and the associated localized Dirac mode. The transition to the nonlinear case is effected by replacing the defect by the auto-induced defect created by nonlinearity, then a soliton localized mode is found to exist. This new specific entity is sustained by the Dirac point rather than photonic bandgaps, so is designated here as Dirac-point soliton. It should be noted that the name, Dirac soliton, has been designated to a soliton of the nonlinear relativistic Dirac equation [23, 24]. We show that the Dirac-point solitons are possible in both Kerr and photorefractive crystals with self-focusing and self-defocusing nonlinearities. Characteristics of the Dirac-point solitons are revealed and their stability is analyzed with linear stability analysis. It is found that the Dirac-point soliton is unstable in a self-focusing lattice, and stable in a self-defocusing lattice, obeys the so-called Vakhitov-Kolokolov (V-K) stability criterion [25]. We verify the stability criterion by direct numerical simulations.

## 2. Band structure and the Dirac point

The propagation of optical beams in a photorefractive material with an optically induced nonlinear photonic lattices is described by the nonlinear Helmholtz equation [10, 11, 26, 27]:

*n*=

*n*

_{e}^{4}

*r*

_{33}

*E*

_{0}/[1 +

*I*(

*x*,

*y*) + |

*E*|

^{2}] represents an index change due to linear modulation and saturable nonlinearity,

*n*and

_{e}*r*

_{33}are the refractive index and electro-optic coefficient of the extraordinary wave,

*E*

_{0}is the applied field, and

*k*

_{0}= 2π/λ

_{0}. In the derivation of Eq. (1), the scalar approximation has been invoked on the basis that both the index modulation and nonlinearity are so small that they can be treated as perturbations to the homogeneous medium. For a highly anisotropic Strontium Barium Niobate crystal SBN:75,

*n*≈2.299,

_{e}*r*

_{33}≈1340pmV

^{−1}. When a static electric field of

*E*

_{0}= 102

*V*/

*mm*is applied, Δ

*n*≤ 3.82 × 10

^{−3}, verifies that Δ

*n*<<

*n*

_{e}^{2}. The linear lattice function Δ

*n*=

*n*

_{e}^{4}

*r*

_{33}

*E*

_{0}/(1 +

*I*) is optically induced by the linear periodic diffraction-free wave with the intensity distribution given by

*I*(

*x*,

*y*). For this study the intensity function is assumed

*I*=

*I*

_{0}[χ + cos(

*b*

_{1}⋅

*r*) + cos(

*b*

_{2}⋅

*r*) + cos(

*b*

_{3}⋅

*r*)]

^{2}, which is a triangular photonic lattice with

*a*being the lattice spacing and

*a*

_{1}=

*aẋ*,

*a*

_{2}=

*a*/2(

*ẋ*+ √3

*ẏ*) being the lattice basis vectors, where

*ẋ*and

*ẏ*are respectively the unit vectors in the

*x*and

*y*directions. The reciprocal basis vectors are

*b*

_{1}= 2π/

*a*(

*ẋ*-√3

*ẏ*/3),

*b*

_{2}= 2π/

*a*(2√3

*ẏ*/3). Such pattern can be generated experimentally, on a stationary background illumination (

*I*

_{0}χ), from the interference of three plane waves possessing an intensity

*I*

_{0}and transverse wave vectors

*b*, where

_{i}*i*= 1,2,3 and

*b*

_{3}is defined as

*b*

_{3}= -

*b*

_{1}-

*b*

_{2}. In fact, the three wave vectors of the plane waves form a triangle, i.e.,

*b*

_{1}+

*b*

_{2}+

*b*

_{3}= 0. The waveguide is uniform in the

*z*direction. Suppose that the wave propagates along the

*z*-direction predominantly, Eq. (1) can be further reduced by factoring out the fast phase variation

*E*=

*U*exp(

*ik*

_{0}

*n*) and dropping the second derivative of the slowly varying envelope function

_{e}z*U*:

*X*=

*x*/

*a*,

*Y*=

*y*/

*a*,

*Z*=

*z*/

*L*has been adopted,

_{d}*L*= 2

_{d}*k*

_{0}

*n*

_{e}a^{2}is the diffraction length of a light beam of width

*a*,

*V*

_{0}=

*k*

_{0}

^{2}

*a*

^{2}

*n*

_{e}^{4}

*r*

_{33}

*E*

_{0},

*V*=

*V*

_{0}/[1 +

*I*(

*X*,

*Y*)] is the linear lattice potential, and

*V*=

_{NL}*V*

_{0}/(

*V*

_{0}/

*V*+ |

*U*|

^{2}) is the nonlinear lattice potential. The new intensity function

*I*(

*X*,

*Y*) has a period normalized to 1. The lattice constant

*a*and the diffraction length

*L*define a basic length scale for the system. If

_{d}*V*

_{0}>0 (<0) the nonlinearity considered here has a self-focusing (self-defocusing) nature. In addition to the saturable nonlinearity, the Kerr nonlinearity can also be incorporated into the above model by resetting

*V*to

_{NL}*V*−

*σ*|

*U*|

^{2}, where

*σ*= 1 (−1) corresponds to a nonlinear self-focusing (self-defocusing) Kerr medium, and

*V*is the linear lattice potential function that is defined as

*V*=

*V*

_{0}[

*χ*+ cos(

*b*

_{1}⋅

*r*) + cos(

*b*

_{2}⋅

*r*) + cos(

*b*

_{3}⋅

*r*)]

^{2}for this study. Four typical examples of the linear lattice potential

*V*are depicted in Fig. 1.

Wave propagation in such a linear periodic lattice *V*(*X,Y*) is known to exhibit unique features that arise from the presence of allowed bands and forbidden gaps. The associated band structure of the lattice can be found by the plane wave expansion method [1]. By substituting a solution of the form *U*(*X,Y,Z*) = *Φ*(*X,Y*)exp(-*iqZ*) into Eq. (2), in the linear limit the following eigenvalue equation is resulted

Eigenfunction *Φ* can be written as a summation of plane waves *Φ*(*r*) = Σ*h*(*G*)exp[*i*(*G* + *k*)*·r*], where *k* is the Bloch momentum, *G* = *b*_{1}*P*_{1} + *b*_{2}*P*_{2}, (*P*_{1},*P*_{2}) are any integers and (*b*_{1}, *b*_{2}) are the reciprocal basis vectors of the lattice. Potential *V* can be expanded as *V*(*r*) = Σ*f*(*G*)exp(*iG·r*), where *f*(*G*) = 1/*S _{cell}* ∫

*V*(

*r*)exp(-

*iG·r*)

*ds*is the expansion coefficient and

*S*= √3/2 is the area of a unit cell. Substituting these expressions into Eq. (3) leads to

_{cell}Given the potential *V*(*r*), the eigenvalue Eq. (4) can be solved numerically to obtain the lattice band structure. The results are shown in Fig. 2 for the potentials of Fig. 1. As can be seen, the band structure of the lattice exhibits Dirac point at the six corners of the Brillouin zone. The eigenvalues at these Dirac points are respectively *q _{D}* = 115.64,

*q*= −8.67,

_{D}*q*= 29.45, and

_{D}*q*= −41.81 for the four different potentials.

_{D}The Dirac point can support guided modes [21]. A defect, which is introduced by resetting *V*(*X,Y*) to *V _{d}* for √(

*X*

^{2}+

*Y*

^{2})≤

*R*, works as a circular waveguide in the configuration. Based on the linear version of Eq. (2), guided modes are discovered by the finite difference beam propagation method (BPM). A source beam

*S*(

*X,Y,Z*) = exp[-(

*X*

^{2}+

*Y*

^{2})]exp(−

*iq*) with a fixed propagation constant

_{D}Z*q*is launched in the waveguide. By changing parameters

_{D}*R*and

*V*of the defect, propagation constant of the waveguide eigenmode varies. When propagation constant of the eigenmode coincides with propagation constant of the source beam, a situation arises in which power monotonously grows with distance, as illustrated in Fig. 3(a). This happens when synchronization of the eigenmode with the source beam is established, so that the eigenmode is always in-phase with the source, and energy is extracted from the source at each step of propagation. This indicates the excitation of an eigenmode of the waveguide with a propagation constant

_{d}*q*=

*q*. In this way, we find the defect mode of the system by carefully designed numerical tests. Field profile of this defect-guided Dirac mode is shown in Fig. 3(b). In Fig. 3(c) |

_{D}*U*| of the mode is multiplied by

*r*

^{3/2}. The almost constant level of oscillation in the tail of the product reveals a 1/

*r*

^{3/2}, asymptotic, algebraic-decay feature of the mode, similar to the Dirac mode found in a triangular photonic crystal [21]. In the second stage of evolution shown in Fig. 3(a), the source is switched off and amplitude of the eigenmode starts a propagation decline. In the free evolution stage, electromagnetic power within the mode decays exponentially according to the format

*P*=

*P*

_{0}exp(−

*αZ*). This law shows that the instantaneous decay rate is

*α*= −1/

*P*d

*P*/d

*Z*and this can be calculated from slope of the numerically obtained power evolution curve. The calculated instantaneous

*α*value is shown in Fig. 3(d) as a function of distance. The average

*α*value is taken for power loss rate of the mode. The BPM yields the space domain response

*U*(

*X,Y,Z*) directly. The spectral response

*u*(

*q*) is subsequently obtained by a discretized Fourier transform from the spatial series:

*u*(

*q*) = 1/N Σ

*U*(0,0,

*n*Δ

*Z*)exp(

*iqn*Δ

*Z*), where Δ

*Z*is step-size, N is the number of steps,

*q*is the propagation constant. The obtained spectrum of the eigenmode is shown in Fig. 3(e), which verifies that the propagation constant of the guided mode is truly centered at

*q*= −41.81.

_{D}Existence of the Dirac mode is attributed to the vanishing density of radiation states at the Dirac point [28]. This means that travelling waves are forbidden in the surrounding triangular lattice at *q* = *q _{D}*. Because of this feature, field concentration around a defect becomes possible and the optical lattices can therefore support defect modes. Such localized states can also exist for a range of propagation constant

*q*around

*q*. When

_{D}*q*≠

*q*the state is no longer a truly localized mode but becomes a guided resonance similar to that observed in [21]. Eigenmode of propagation constant

_{D}*q*other than

*q*can be found in a similar way by setting the source beam to

_{D}*S*(

*X,Y,Z*) = exp[-(

*X*

^{2}+

*Y*

^{2})]exp(−

*iqZ*) while adjusting the parameters of the defect for synchronization. The results are shown in Figs. 4(a)-4(d), corresponding to the lattice potentials, respectively, of Figs. 1(a)-1(d). They show parameter

*V*of the defect versus the propagation constant

_{d}*q*of the eigenmode for fixed

*R*. Leakage of the Dirac mode can be studied numerically by the BPM conveniently. The average loss rates

*α*, of the localized modes, extracted from free evolution of the power are shown also in Figs. 4(a)-4(d) as a function of the propagation constant

*q*using, respectively, the four different lattices. In the absence of material absorption, there still exists two loss mechanisms for the Dirac mode, namely, leakage

*α*to the surrounding medium due to the finite lattice size (called losses associated with field penetration across boundary into the surrounding), and leakage

_{p}*α*due to scattering into the continuum of states (called losses associated with coupling into radiation modes). The net decay rate is the sum

_{s}*α*=

*α*+

_{p}*α*. At the Dirac point losses

_{s}*α*associated with scattering into radiation modes vanish. Any residual values of losses at this point,

_{s}*α*, are entirely due to the finite lattice sizes but they can be made as small as desired by increasing the boundary surrounding the defect. As the propagation constant moves away from the Dirac point, the density of radiation states increases in a linear fashion. Then, in addition to leakage

_{p}*α*to the surrounding medium, there can be leakage

_{p}*α*to the continuum of states. The net decay rate

_{s}*α*increases with |

*q*-

*q*|. Minima of the loss rates

_{D}*α*in Fig. 4 occur around the corresponding

*q*, thus confirming the positions of the Dirac points. As revealed in Fig. 4, there exists a small range

_{D}*q*around

*q*where decay rate is low and waveguiding by a defect is practically realizable. A linear damping length

_{D}*L*can be defined for the localized mode through exp(−

_{lin}*αL*) = 0.9, which represents the propagation distance of the mode when it is damped to 90% of its initial amplitude. In the literature damping length is often defined at 1/

_{lin}*e*or 50% of the original amplitude. We use a threshold value of 90% here to leave more margin for experiments. This value also facilitates comparison with soliton propagation distance

*L*introduced later in Fig. 13, where amplitude drops to about 90% as the Dirac point soliton degrades into a gap soliton. By saying the decay rate is small and waveguiding by the defect is possible, the wave is required to propagate for more than 5 diffraction lengths [7]. Using the condition

_{sol}*L*≥5 in the definition of damping length exp(−

_{lin}*αL*) = 0.9 leads to a critical decay rate

_{lin}*α*= −ln0.9/5≈0.021. This critical value of

_{c}*α*is indicated in Fig. 4 by red horizontal lines, and results in the region

*q*around

*q*where

_{D}*α*is below

*α*, is marked red in the figure.

_{c}## 3. The Dirac-point soliton

The defect that accommodates the above discussed Dirac modes could be created by nonlinearly-induced, onsite, index change. If the defects are self-induced optically by nonlinearity, the corresponding, self-localized nonlinear eigenmodes are referred to as solitons. In other words, the presence of a defect mode at the Dirac point in the linear limit suggests the existence of a Dirac-point soliton in the corresponding nonlinear regime. This is, indeed, the case. To find the Dirac-point solitons of Eq. (2), a solution is sought with the form *U*(*X,Y,Z*) = *Φ*(*X,Y*)exp(-*iqZ*). Using this reduce Eq. (2) to

Equation (5) is the static nonlinear Schrödinger equation, and yields solitary wave solutions as the outcomes of a, numerical, modified squared-operator method [29]. For a fundamental soliton *Φ* is real and *Φ*(∞) = 0. Apart from the fundamental soliton, Eq. (5) also supports a vortex soliton of the form: *Φ* = *f*(*r*)*e ^{imθ}*, where the integer

*m*stands for a phase twist around the intensity ring. Some typical Dirac-point solitons of Eq. (5), including fundamental and vortex solitons, are shown in Figs. 5(a)-8(a) and Figs. 5(b)-8(b) respectively for the focusing and defocusing saturable/Kerr nonlinearities. Characteristics show that the Dirac-point solitons belong to a group labelled as algebraic solitons [12, 30, 31]. As shown in Figs. 5(c)-8(c) and Figs. 5(d)-8(d), tail of the Dirac-point soliton oscillates while decays asymptotically and algebraically according to a power-law (roughly

*r*

^{−3/2}) at large distances. It is understandable because towards the end of the tail soliton field is so weak that it resumes a linear behavior, and in the linear limit an optical lattice can support waves with power-law asymptotics at the Dirac point, as demonstrated in section 2. Total power conveyed by the soliton

*P*= ∫∫|

*Φ*|

^{2}

*r*d

*r*d

*θ*versus the eigenvalue

*q*is shown in Figs. 5(e)-8(e).

The soliton results build upon earlier results in section 2 that show that localized states may exist at the Dirac point in a triangular lattice in the presence of defects. Here the defects are replaced with self-induced defects, achieved through the presence of nonlinearity, so, soliton is somewhat expected. Such states constitute a new addition to the soliton family. It should be stressed that an optical defect is necessary to anchor the localized mode. After all, in the absence of a defect, the system is periodic and all the modes are extended Bloch modes. In the nonlinear case nonlinearity provides an effective defect.

Strictly speaking, the Dirac-point soliton exists only at a single point in the [*q*, *P*] space. In practice, however, it can be excited within a small range around *q _{D}*. When

*q*≠

*q*the soliton becomes a nonlinear guided resonance similar to the linear guided resonance at these

_{D}*q*values. The nonlinear guided resonance can leak at a small rate into the continuum of states that makes up the band. Nevertheless, the mode still concentrates much of its field energy near the defect. Loss of the nonlinear guided resonance is estimated at about the same level of loss of the linear guided resonance in the same potential, so the range of existence is allocated to the valley of the linear loss shown in Fig. 4. Outside the valley, loss increases quickly and disturbs the solitons. To indicate the level of linear losses, the net decay rate

*α*of the linear wave is superimposed in Figs. 5(e)-8(e) and Figs. 5(f)-8(f) as background. The Dirac-point solitons are expected to exist roughly within the ranges shaded light grey, where level of the linear decay rate

*α*is below

*α*. Beyond the shaded light grey range heavy losses may seriously restrict propagation distance of the Dirac-point solitons.

_{c}## 4. Stability analysis

The Dirac-point solitons are nonlinear stationary solutions of Eq. (2) so that preserve their shapes upon evolution in the media, but their stability is not guaranteed because of the non-integrable nature of the underlying equation. In fact, their stability against small perturbations is a crucial issue because only stable (or weakly unstable) self-trapped beams can be observed experimentally. To study the stability of these Dirac-point solitons, a perturbation, with the form *U* = *e*^{-}* ^{iqZ}*[

*Φ*+ (

*v*-

*w*)

*e*+ (

^{λZ}*v*+

*w*)*

*e*], is invoked, where

^{λ*Z}*v*,

*w*<<1. Substituting this expression into Eq. (2) and then linearizing, results in an eigenvalue problem of the form

*LΨ*=

*λΨ*, here

*Ψ*= [

*v, w*]

^{T}is the eigenfunction,

*λ*is the eigenvalue. Expressions of the operator

*L*for Kerr and saturable nonlinearities are respectively

*L*

_{0}= ∂

^{2}/∂

*x*

^{2}+ ∂

^{2}/∂

*y*

^{2}-

*V*+

_{NL}*q*. The eigenvalue equation

*LΨ*=

*λΨ*is solved by a numerical iteration method [32] for growth rate. The real part of the soliton perturbation growth rate is Re(

*λ*) and this is plotted in Figs. 5(f)-8(f) as function of the propagation constant of the soliton for respectively the four different lattices shown in Figs. 1(a)-1(d). It can be seen that in a saturable self-focusing lattice (Fig. 5), growth rates exceed 10 and both the fundamental and the first vortex solitons are unstable. In a saturable self-defocusing lattice (Fig. 6), growth rates are of the order of 1 so the fundamental and the first vortex solitons are weakly unstable. In Kerr lattices (Figs. 7 and 8), the fundamental solitons for self-focusing nonlinearity have excessively large growth rates (>40) so are surely unstable. By contrast, the fundamental solitons in self-defocusing nonlinearity have much smaller growth rates (∼0.1) so are referred to as being quasi-stable. The power curves in Figs. 5(e)-8(e) give information about the stability of the solitons also. According to the V-K stability criterion [25], a soliton is stable (unstable) if the slope of the corresponding power curve is positive (negative). The V-K stability criterion is derived for homogenous nonlinear medium but a periodic potential is present here, so it is not directly applicable. However, the depth of index modulation is very shallow. As such, the linear stability analysis comes up with results that agree well with the V-K stability criterion. The Dirac-point solitons in self-defocusing lattices [Figs. 6(a), 6(b), Fig. 7(b), Fig. 8(b)] exhibits d

*P*/d

*q*≥0 on the power curve and are quasi-stable (or weakly unstable), whereas the Dirac-point solitons in self-focusing lattices [Figs. 5(a), 5(b), Fig. 7(a), Fig. 8(a)] exhibits d

*P*/d

*q*<0 and are unstable. The Dirac-point solitons are confirmed to obey the simple V-K stability criterion despite the potential.

The stability of Dirac-point solitons can be checked numerically using the split-step Fourier method based on Eq. (2). Some simulated evolution scenarios of typical stable and unstable solitons are shown in Figs. 9-12. The dimension of evolutions can be understood by recalling that *X* = *x*/*a*, *Y* = *y*/*a*, *Z* = *z*/*L _{d}*. Using data from a real experiment [10, 26]:

*a*= 11

*μm*,

*n*≈2.299,

_{e}*λ*

_{0}= 0.5

*μm*, the diffraction length

*L*= 2

_{d}*k*

_{0}

*n*

_{e}a^{2}is about 6.99

*mm*in real space. In Fig. 9 it is shown that the first vortex Dirac-point soliton in a saturable self-focusing lattice is unstable. It breaks down after

*Z*>1. In Fig. 10 it is shown that the first vortex Dirac-point soliton in a saturable self-defocusing lattice is weakly unstable. The breakdown distance of the weakly unstable soliton (Fig. 10) is about an order of magnitude larger than that of the unstable soliton (Fig. 9). In Kerr lattices the fundamental soliton is unstable for self-focusing nonlinearity (Fig. 11) and quasi-stable for self-defocusing nonlinearity (Fig. 12). The breakdown distance of the quasi-stable soliton (Fig. 12) is about three orders of magnitude larger than that of the unstable soliton (Fig. 11). These results confirm the stability analysis that has been deployed in this investigation. The broadening of propagation constant spectrum is attributed to a phenomenon called self-propagation-constant shift. Unlike the eigenmode of a linear waveguide that has a fixed propagation constant, soliton, which is an eigenmode of a nonlinearly-induced waveguide, has a changeable propagation constant. As the Dirac-point soliton loses power its amplitude reduces. This alters the parameters of the nonlinearly-induced defect, and, in turn, shifts its propagation constant. This process is an analogue of the self-frequency shift of a temporal soliton. For the case of d

*P*/d

*q*>0 (<0),

*q*is smaller (larger) for soliton with smaller amplitude. So, as amplitude drops in propagation the corresponding

*q*continuously shifts downwards (upwards). This explains why in Figs. 9-12 the propagation constant of the soliton sweeps through a wide range of values. For the particular lattice of Fig. 9 (Fig. 12), there is a photonic band gap above (below) the Dirac point, as shown in Fig. 2(a) [Fig. 2(d)]. As the propagation constant

*q*shifts upwards (downwards) into the band gap, the Dirac-point soliton degrades into a gap soliton and settles down there. This is confirmed by the second peak in the spectrum as shown in Fig. 9(b) [Fig. 12(b)]. When the Dirac point is surrounded by bands of extended Bloch modes, as in the case of Figs. 10 and 11, the Dirac-point solitons break up into radiation waves in propagation.

## 5. Discussion

It is revealed by Figs. 9-12 that a stable soliton is not fully stable but has a finite propagation distance. For example, it is shown in Fig. 12 that the quasi-stable Dirac-point soliton degrades after about 40 diffraction lengths. Indeed, real part of the perturbation growth rate Re(*λ*) shown in Fig. 8(f) for this soliton does not drop to exact zero at *q _{D}* = −41.81 but has a small residue Re(

*λ*) = 3.81 × 10

^{−2}. As discussed in section 2, the linear defect mode at the Dirac point experiences a propagation decline due to losses. Since soliton can be considered as a self-induced defect state, from this analogue picture, it cannot be fully stable while the corresponding linear defect state is weakly unstable. Similar to the linear Dirac modes, leakage of the Dirac-point soliton consists of leakage into the surrounding medium and leakage into radiation states. The vanishing density of radiation states at the Dirac point can suppress the decay rate leading to long-lived, soliton-like propagation. If window of cladding is made larger, losses associated with field penetration across boundary into the surrounding will be made smaller, and propagation distance of the quasi-stable Dirac-point soliton can be further extended to almost as long as desired.

Recently, a type of nonlinear leaky mode called soleakon has been studied [33, 34]. The soleakon features oscillating delocalized tails, indicating the outgoing flow of energy to the cladding. Such undamped oscillating tail also leads to a diverging total power of the mode, i.e., ∫|*Φ*|^{2}*r*d*r*d*θ*→∞. However, our finding is a new phenomenon and does not fit the description of the so-called soleakon. The Dirac-point soliton occurs at a conical singularity, where the density of radiation states vanishes so, coupling between a Dirac-point soliton and radiation states is actually suppressed at the Dirac frequency. The residue damping is due to the finite lattice size and can be reduced by enlarging thickness of cladding. Because of this reason, leakage of the Dirac-point soliton is accidental and can be removed, whereas the soleakon is intrinsically leaky. Furthermore, oscillating tail of the Dirac-point soliton is localized as it damps asymptotically according to an algebraic law 1/*r*^{3/2}. As a matter of fact, tail of the Dirac-point soliton decays faster than 1/*r*, the quantity |*Φ* | is therefore square integrable, i.e., *P* = ∫∫|*Φ* |^{2}*r*d*r*d*θ* is finite. This leads to a converging total power, and confirms that the Dirac-point soliton at *q*≠*q _{D}* is a well-behaved, true-bound state.

Degradation of the quasi-stable Dirac-point soliton in propagation is understandable as loss is in existence. The losses will sooner, or later, breakup the balance between contraction and diffraction, and eventually diminish the soliton. By examining numerically the evolution scenarios, the finite propagation distance *L _{sol}* of the Dirac-point soliton that shown in Fig. 8(b) is estimated and the obtained value is plotted in Fig. 13 as a function of the initial eigenvalue

*q*(initial propagation constant). The corresponding linear damping length

*L*for the Dirac mode in the same lattice is plotted in Fig. 13 also, alongside the soliton propagation distance

_{lin}*L*. It is found that

_{sol}*L*matches

_{sol}*L*on the left hand side of the Dirac point. However, on the right hand side it does not follow the trend of

_{lin}*L*and, surprisingly, remains at a rather high level, implies that the Dirac-point soliton in this case can survive for longer distance than its linear counterpart. The extended propagation distance of the Dirac-point soliton on the right hand side of the Dirac point is attributed to the phenomenon of self-propagation-constant shift. The correlation of the linear damping length

_{lin}*L*and the propagation distance

_{lin}*L*of the Dirac-point soliton as illustrated in Fig. 13 can be understood as following: For the Dirac-point soliton starts on the left-hand side of the Dirac point, its

_{sol}*q*departs from

*q*straight away so the soliton decays into a gap soliton quickly. In this process |

_{D}*q*-

*q*| increases, and so does the loss rate. This is why the propagation distance

_{D}*L*of the soliton is below the linear damping length

_{sol}*L*on this side. On the other hand, for the Dirac-point soliton starts from the right hand side of the Dirac point, its

_{lin}*q*approaches

*q*from the right so |

_{D}*q*-

*q*| reduces in the initial stage. As |

_{D}*q*-

*q*| reduces, losses associated with coupling into radiation modes decreases, so the soliton encounters a weaker damping than that does the corresponding linear mode encounter, and consequently, it lasts longer in propagation than its linear counterpart. This explains why the propagation distance

_{D}*L*of the soliton is above the linear damping length

_{sol}*L*on the right hand side of the Dirac point.

_{lin}## 6. Conclusions

A new type of soliton is discovered for photonic lattices. This new soliton relies upon a Dirac point, rather than photonic bandgaps, in order to establish field localization. We believe the outcomes reported here will have an immediate impact because of the unique mechanism underpinning the soliton. The Dirac equation is a special symbol of relativistic quantum mechanics, the investigations of the Dirac-point soliton, given here, will lead to new findings in the area of relativistic quantum effects on the transport of photons, phonons, and electrons.

## Funding

Leading Talents of Guangdong Province Program; National Natural Science Foundation of China (11574070, 11404087, 51771186); Fundamental Research Funds for the Central Universities of China; Postdoctoral Science Foundation (2015M571918, 2017T100442); The European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement (744817).

## References and links

**1. **M. Skorobogatiy and J. Yang, *Fundamentals of Photonic Crystal Guiding* (Cambridge University, 2009).

**2. **G. Li, M. Zeisberger, and M. A. Schmidt, “Guiding light in a water core all-solid cladding photonic band gap fiber - an innovative platform for fiber-based optofluidics,” Opt. Express **25**(19), 22467–22479 (2017). [PubMed]

**3. **C. H. Lai and C. I. Hsieh, “Power coupling between gain-guided index-antiguided planar waveguides,” Opt. Express **24**(24), 28026–28038 (2016). [PubMed]

**4. **J. Buencuerpo, J. M. Llorens, P. Zilio, W. Raja, J. Cunha, A. Alabastri, R. P. Zaccaria, A. Martí, and T. Versloot, “Light-trapping in photon enhanced thermionic emitters,” Opt. Express **23**(19), A1220–A1235 (2015). [PubMed]

**5. **J. C. Knight, J. Broeng, T. A. Birks, and P. S. J. Russell, “Photonic band gap guidance in optical fibers,” Science **282**(5393), 1476–1478 (1998). [PubMed]

**6. **K. C. Kao and G. A. Hockham, “Dielectric–fibre surface waveguides for optical frequencies,” IEE Proceedings 113, 1151–1158 (1966).

**7. **Y. S. Kivshar and G. P. Agrawal, *Optical solitons: from fibers to photonic crystals*, (Academic, 2003).

**8. **M. Peccianti, C. Conti, G. Assanto, A. De Luca, and C. Umeton, “Routing of anisotropic spatial solitons and modulational instability in liquid crystals,” Nature **432**(7018), 733–737 (2004). [PubMed]

**9. **D. N. Christodoulides, F. Lederer, and Y. Silberberg, “Discretizing light behaviour in linear and nonlinear waveguide lattices,” Nature **424**(6950), 817–823 (2003). [PubMed]

**10. **J. W. Fleischer, M. Segev, N. K. Efremidis, and D. N. Christodoulides, “Observation of two-dimensional discrete solitons in optically induced nonlinear photonic lattices,” Nature **422**(6928), 147–150 (2003). [PubMed]

**11. **N. K. Efremidis, S. Sears, D. N. Christodoulides, J. W. Fleischer, and M. Segev, “Discrete solitons in photorefractive optically induced photonic lattices,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **66**(4 Pt 2), 046602 (2002). [PubMed]

**12. **C. Brunhuber, F. G. Mertens, and Y. Gaididei, “Long-range effects on superdiffusive algebraic solitons in anharmonic chains,” Eur. Phys. J. B **57**, 57–65 (2007).

**13. **A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, “The electronic properties of graphene,” Rev. Mod. Phys. **81**, 109–162 (2009).

**14. **S. Raghu and F. D. M. Haldane, “Analogs of quantum-Hall-effect edge states in photonic crystals,” Phys. Rev. A **78**, 033834 (2008).

**15. **S. Longhi, “Quantum-optical analogies using photonic structures,” Laser Photonics Rev. **3**, 243–261 (2009).

**16. **S. R. Zandbergen and M. J. A. de Dood, “Experimental observation of strong edge effects on the pseudodiffusive transport of light in photonic graphene,” Phys. Rev. Lett. **104**(4), 043903 (2010). [PubMed]

**17. **X. Zhang, “Observing Zitterbewegung for photons near the Dirac point of a two-dimensional photonic crystal,” Phys. Rev. Lett. **100**(11), 113903 (2008). [PubMed]

**18. **M. J. Ablowitz and Y. Zhu, “Nonlinear Waves in Shallow Honeycomb Lattices,” SIAM J. Appl. Math. **72**, 240–260 (2012).

**19. **O. Bahat-Treidel, O. Peleg, and M. Segev, “Symmetry breaking in honeycomb photonic lattices,” Opt. Lett. **33**(19), 2251–2253 (2008). [PubMed]

**20. **O. Peleg, G. Bartal, B. Freedman, O. Manela, M. Segev, and D. N. Christodoulides, “Conical diffraction and gap solitons in honeycomb photonic lattices,” Phys. Rev. Lett. **98**(10), 103901 (2007). [PubMed]

**21. **K. Xie, H. M. Jiang, A. D. Boardman, Y. Liu, Z. H. Wu, M. Xie, P. Jiang, Q. Xu, M. Yu, and L. E. Davis, “Trapped photons at a Dirac point: a new horizon for photonic crystals,” Laser Photonics Rev. **8**, 583–589 (2014).

**22. **K. Xie, W. Zhang, A. D. Boardman, H. Jiang, Z. Hu, Y. Liu, M. Xie, Q. Mao, L. Hu, Q. Li, T. Yang, F. Wen, and E. Wang, “Fiber guiding at the Dirac frequency beyond photonic bandgaps,” Light Sci. Appl. **4**, e304 (2015).

**23. **Y. Nogami, F. M. Toyama, and Z. Zhao, “Nonlinear Dirac soliton in an external field,” J. Phys. Math. Gen. **28**, 1413–1424 (1995).

**24. **X. Tran, S. Longhi, and F. Biancalana, “Optical analogue of relativistic Dirac solitons in binary waveguide arrays,” Ann. Phys. **340**, 179–187 (2014).

**25. **N. G. Vakhitov and A. A. Kolokolov, “Stationary solutions of the wave equation in a medium with nonlinearity saturation,” Radiophys. Quantum Electron. **16**, 783–789 (1973).

**26. **N. K. Efremidis, J. Hudock, D. N. Christodoulides, J. W. Fleischer, O. Cohen, and M. Segev, “Two-dimensional optical lattice solitons,” Phys. Rev. Lett. **91**(21), 213906 (2003). [PubMed]

**27. **T. J. Alexander, A. S. Desyatnikov, and Y. S. Kivshar, “Multivortex solitons in triangular photonic lattices,” Opt. Lett. **32**(10), 1293–1295 (2007). [PubMed]

**28. **M. Plihal and A. A. Maradudin, “Photonic band structure of two-dimensional systems: The triangular lattice,” Phys. Rev. B Condens. Matter **44**(16), 8565–8571 (1991). [PubMed]

**29. **J. Yang and T. I. Lakoba, “Universally-convergent squared-operator iteration methods for solitary waves in general nonlinear wave equations,” Stud. Appl. Math. **118**, 153–197 (2007).

**30. **H. Alatas, A. A. Iskandar, M. O. Tjia, and T. P. Valkering, “Rational solitons in deep nonlinear optical Bragg grating,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **73**(6 Pt 2), 066606 (2006). [PubMed]

**31. **J. M. Christian, G. S. McDonald, and P. Chamorro-Posada, “Helmholtz algebraic solitons,” J. Phys. A Math. Theor. **43**, 085212 (2010).

**32. **J. Yang, “Iteration methods for stability spectra of solitary waves,” J. Comput. Phys. **227**, 6862–6876 (2008).

**33. **O. Peleg, Y. Plotnik, N. Moiseyev, O. Cohen, and M. Segev, “Self-trapped leaky waves and their interactions,” Phys. Rev. A **80**, 041801 (2009).

**34. **M. Kozlov, O. Kfir, and O. Cohen, “Self-trapped leaky waves in lattices: discrete and Bragg soleakons,” Opt. Express **21**(17), 19690–19700 (2013). [PubMed]