Solitary waves in three-dimensional crystal-like structures

The motion of three-dimensional (3D) solitary waves and solitons in nonlinear crystal-like structures, such as photonic materials, is studied. It is demonstrated that collective excitations in these systems can be tailored to move in particular directions of the 3D system. The effect of modulation instability is studied showing that in some cases it can be delayed by using a lensing factor. Analytical results supported by numerical simulations are presented.


I. INTRODUCTION
Recently, great progress has been made in the experimental fabrication of three-dimensional (3D) crystal-like optical nano structures, such as photonic crystals [1] and photonic metamaterials [2]. Light in 3D crystal-like optical nano structures can interact with the regular pattern of the structure setting up resonances. These resonances can cause beams and pulses to be deflected in unconventional directions and even to slow down the speed of light.
Engineered optical metamaterials with unique electromagnetic properties, have become in recent years a hot research topic because of their interesting physics and exciting potential applications as, e.g cloaking [3,4], or negative refractive index [5], among others.
The theoretical analysis of these systems have been mostly performed with the help of numerical simulations, where the 3D spatial distribution of the effective electromagnetic properties of the material medium (such as the permittivity and permeability) are tailored to have specific electromagnetic properties in the continuum limit [3]. Other discrete effects of the system on the propagation of light waves are usually neglected. Similar approximations have been also adopted for studying cloaking of matter waves [4]. So far, most of the studies have been done for linear planewaves, so nonlinear effects have been neglected. There is, however, strong evidence from low-dimensional systems [6,7,8,9] that discrete effects in combination with nonlinearities may play an important role in describing realistic 3D crystallike structures. So, nonlinear excitations such as moving solitary waves (in short solitons) and vortices can be expected.
With respect to the theoretical models, it has been shown that electromagnetic waves interacting with metamaterials in the tight binding limit can be described as a 3D lattice of microresonators modeled by the 3D discrete nonlinear Schrödinger equation (3D-DNLSE) [5]. In the case of 3D photonic crystals, so far, no discrete model has been proposed. However, it is well known that the dynamics of light beams in photonic waveguides can be described by the 2D-DNLSE [8,9], where the discreteness is transversal to the light propagation. Moreover, it has been also shown that trapped light waves travelling along chains of optical resonators can be described by the 1D DNLSE [6,7], where the discreteness is longitudinal to the light propagation. These two low-dimensional models strongly suggest that 3D photonic crystals with high index of refraction [1,2] may be effectively modeled in the tight binding limit by a 3D lattice of optical resonators governed by the 3D-DNLSE Other parameters are U = −1, J x = J y = J z = 1, γ x = γ y = γ z = 1/3. Here, β x = β y = β z = +1, so only sech-type solutions in Eq. (8) are considered. [10,11,12].
Notice that the 3D-DNLSE is an ubiquitous dynamical-lattice model which may emerge from a variety of other important problems and has its direct physical realization in Bose-Einstein condensates (BECs) trapped in strong optical lattices [10,11,12].
The aim of the present work is to study the dynamics of moving solitons and vortices in nonlinear 3D-crystal-like structures described by the 3D-DNLSE.

II. THEORY
The general form of the 3D-DNLSE with coupling constants J x , J y and J z is where ρ m,n,p = |ψ m,n,p | 2 can be interpreted as a intensity (e.g. in crystals built of microresonators) or as a probabilitydensity (in BEC arrays). In Eq. (1) the nonlinear coefficient U is a real constant and t is the time coordinate.
In order to proceed we shall consider an expansion into harmonics of a travelling wave ansatz for an envelope complex function, i.e ψ m,n,p = ∞ µ,ν,ξ=1 where and and Ω µ,ν,ξ = Ω 0 (µω x + νω y + ξω z ) is the frequency (chemical potential in BEC arrays).
Since the coordinates {x, y, z} are mutually independent, Eq. (6) can be integrated. So, finally, approximate soliton solutions of Eq. (1) read as where δ βw,±1 is the Kronecker delta with β w = sign(γ w α 0 /α w ) and w = x, y, z. In Eq. (8) the soliton widths read as and where the relation γ x + γ y + γ z = 1 should be satisfied. Besides and In Fig. 1 it is shown an example of the motion of solitons following from initial conditions given by Eq. (8).
Finally we obtain the solution where the embedded vorticity is given by the order parameter µ. In Eq.
In Figs. 1 and 3 we observe that solitons and vortices when propagating undergo a selfdefocusing instability. This inestability follows from the interplay between nonlinearity and discreteness. Notice that nonlinearity tends to localize the wave while the discrete diffraction tends to spread it out. For moving solitons the self-defocusing instability corresponds to the case when diffraction effect is stronger than the nonlinear localization effect. This selfdefocusing process manifest itself as a soliton-shape broadening accompanied by exponentiallike amplitude decay. An example of this amplitude decay can be observed in Fig. 2(a), where the evolution of the absolute ρ m,n,p maximum (ρ max ) of the collision in Fig. 1 is plotted.
Notice that in Fig. 2(a) ρ max peaks when the collision in Fig. 1 [label D (green color)] occurs. After collision the solitons further broaden in shape, as can be observed in Fig.   1 [labels A ′ , B ′ , and C ′ (blue color)]. Besides, for comparison in Fig. 2(a) the amplitude ρ max of the individual solitons in the absence of collision is plotted. We can observe that before and after collision the ρ max value coincides with that of individual solitons. Moreover, we can observe that the peak amplitude during collision is much higher than three times the amplitude of the individual solitons. It is because during collision the nonlinear effect becomes strong, so a transitory self-focusing process takes place.
The self-defocusing process in solitons depends on their form and direction of motion, which in turn are governed by k. As in two dimensional systems [8,9], solitons moving along the diagonal directions, are less prone to self-defocusing than those moving along the main-axes directions, Notice that it is possible to control the defocusing process so that solitons moving in different directions and velocities can have the same decay rate. For example, solitons in Fig. 1 have similar decay rates, as shown in Fig. 2(a). This can be easily achieved by first choosing with the help of k the wanted directions and velocities of motion and then with the help of Ω 0 , included in α 0 in Eqs. (7), the same initial amplitudes. Here we use the fact that the amplitude of a soliton is proportional to its decay rate [9,14]. Notice that for similar decay rates we obtain solitons with different sizes and velocities, as shown Figs. 1 and 3.
In Fig. 3 we plot the initial soliton and vortex solutions [labels A, B, and C (color red)] following from Eq. (21). In particular, the motion of a vortex shell (A → D → A ′ in Fig. 3) along its azimuthal axis (p axis) with vorticity µ = 1 is shown. The soliton dynamics is similar as in Fig. 1, however small perturbations of the soliton shape are observed after the collision (labels B ′ , and C ′ in Fig. 3). The broadening and breaking-apart of the vortex in Fig. 3 is due to the modulation instability and can be also observed in the absence of collisions. This instability effect on vortices is more pronounced for higher values of the vorticity (µ ≥ 2).
We remark that the center of mass of both solitons and vortices move with a constant velocity whose magnitude measured in simulations is very well predicted by Eq. (13). On the other hand, though not seen in Figs. 1 and 3, the presence of low-intensity radiation tails can be detected for both soliton and moving vortices.
A question that emerges from the results above is how to mitigate the self-defocusing effect. In order to tackle this problem we investigate the effect of a magnification in the form of a 3D thin-lens phase, i.e. ψ m,n,p → ψ m,n,p exp(−i r 2 /(4T )).
In the present analysis ψ m,n,p solutions follow from Eq. (21), where the radial coordinate r is given by Eq. (16). In Eq. (25), due to the discreteness of the system, the parameter T does not correspond exactly to the focal length, as it happens in continuous systems.
In Fig. 4(a) we consider three solitons moving in the positive direction of the m axis, and their individual ρ max values are plotted in Fig. 4(b). In Fig. 4(a) the transversal separation distance between the solitons is large enough to avoid soliton-soliton interactions. The initial soliton forms, ψ [Labels C, A, B (red color) in Fig. 4(a)], are identical to each other but distinct from the T value in the imposed thin-lens phase. In order to compare the soliton behavior, the level value of the isosurfaces in Fig. 4(a) was defined to be fixed to some initial value, as shown in 4(b) (red line). So, solitons with ρ max value above this level appear plotted in Fig. 4(a). Notice, e.g., that in the soliton evolution A → A ′ → A ′′ in Fig. 4(a) the soliton shape vanishes at A ′′ because its ρ max value [black line in Fig. 4(b)] decays below the level value [red line in Fig. 4(b)] of the isosurface for t > 12.
Three different T values have been considered in Fig. 4 However, it is to remark that for larger time scales a self-defocusing process is observed and unavoidable. The value T = 0.1 (C → C ′ → C ′′ ) corresponds to a case where a focusing due to the thin-lens phase is immediately followed by a self-focusing effect, as can be observed in Fig. 4(b) (purple dotted line). Besides, a strong radiation tail can be also observed in For further comparison in Fig. 4(b) the case of T = 0.05 (green dot-dashed line) is also plotted. This case is similar to the case T = 0.1 where the first peak corresponds to a focusing due to the imposed thin-lens phase and the second peak is due to a self-focusing effect. Notice that after the second peak the soliton amplitude strongly decays and a long radiation tail, more in the fashion of Fig. 4(a) (C ′′ ), can be also observed.