1. Introduction

Spatial soliton solutions in 1D periodic structures have been theoretically predicted [1] and experimentally demonstrated in nonlinear waveguide arrays [2]. More recently, similar structures have been observed in optically induced gratings [3, 4] and, for the first time, in 2D photonic lattices [5]. Modeling has been traditionally performed using the discrete nonlinear Schrödinger equation (NLSE)[1, 6], valid only in the so-called tight-binding approximation. However, the accurate study of nonlinear solutions in optically-induced lattices requires the resolution of the NLSE with a continuous periodic potential model [3, 5]. Discrete optical vortices have been also predicted in the discrete NLSE [7] and in continuous models [8] and, recently, experimentally observed in optically-induced 2D lattices [9].

A different scenario appears in a 2D photonic crystal in the presence of a defect. An example of such a system is a photonic crystal fiber (PCF). A PCF is a silica fiber that possesses a regular array of air-holes extending along the entire fiber length (see Fig. 1(a)). It is an axially invariant 2D photonic crystal with a central defect (the PCF core, where guidance occurs). It has been recently proven that a PCF acting as a 2D nonlinear photonic crystal can support and stabilize fundamental solitons [10]. Unlike fundamental solitons embedded in perfectly periodic potentials [5, 8], a distinguishing feature of these systems is that no power gap is needed to generate the nonlinear localized solution. In this letter we go a step further and give numerical evidence of the existence of vortex solitons in PCF’s.

2. Description of the method

When one considers that silica can have a nonlinear response, electromagnetic propagation in a PCF is given by the following non-paraxial equation (∇·E≈0):

${\nabla }_t^2$ being the 2D-transverse Laplacian operator and k ₀=ω ₀/c the vacuum wavenumber. The linear refractive index profile function n ₀(x,y) is 1 in the air-holes and equals nsilica in silica, whereas the nonlinear index profile function n ₂(x,y) is different from zero only in silica ($n_{2\left(\text{silica}\right)}^2$≡3χ ⁽³⁾(_silica)/(2ε ₀ cn _0(silica))). We search for monochromatic (or quasi-monochromatic) stationary electric field solutions with well-defined constant polarization: $\mathbf{E}(x,y,z,t)\approx \mathbf{u}\varphi (x,y)e^{i\left(\beta z-{\omega }_{0}t\right)}$. That is, we analyze the following equation:

where L ₀=${\nabla }_t^2$+$k_0^2$ $n_0^2$(x,y) and L_NL (|ϕ|)=$k_0^2$ $n_2^2$(x,y)|ϕ|² stand for the linear and nonlinear parts of the differential operator L(|ϕ|), respectively. Due to the geometry of the air-holes distribution, both n ₀(x,y) and n ₂(x,y) are invariant under the action of the 𝓒 _6v point-symmetry group. This group is constituted by discrete π/3-rotations (R _π/3) plus specular reflections with respect to the x and y axes: $y\stackrel{R_x}{\to }-y$ and $x\stackrel{R_y}{\to }-x$, or, equivalently, $\theta \stackrel{R_x}{\to }-\theta$ and $\theta \stackrel{R_y}{\to }\pi -\theta$, in polar coordinates.

Our aim is to look for nonlinear solutions of Eq. (2) beyond the fundamental solitons reported in Ref. [10]. Our search strategy will be based on group theory arguments, which will be enlightening in the classification and understanding of new solutions according to the sixth-fold symmetry of the system. In the linear case, this approach has been adopted to study the degeneracy of solutions in a PCF [11]. Here, we will use it to obtain analytical expressions for the angular dependency of solutions. According to group theory, since L ₀ is invariant under the 𝓒 _6v group, all its eigenfunctions have to lie on discrete representations of this group [12, 13]. Solutions of the vortex type belong to one of the two two-dimensional representation of 𝓒 _6v. A pair of functions (ϕ_l,ϕ*_l) (l=1, 2) belonging to these representations have the following transformation properties: $({\varphi }_l,{\varphi }_l^{*})\stackrel{R_{\pi ⁄3}}{\to }({\epsilon }^l{\varphi }_l,{\epsilon }^{*l}{\varphi }_l^{*})\epsilon =e^{i\pi ⁄3}$, $({\varphi }_l,{\varphi }_l^{*})\stackrel{R_{\pi }}{\to }({\varphi }_l^{*},{\varphi }_l)$ and $({\varphi }_l,{\varphi }_l^{*})\stackrel{R_y}{\to }{\left(-1\right)}^l({\varphi }_l^{*},{\varphi }_l)$. These symmetry constraints impose the form of these functions:

where ${\varphi }_l^s$ is a scalar function, that is, a function invariant under finite rotations and specular reflections, whereas ${\varphi }_l^p$ is a pseudo-scalar function, invariant under rotations but sign-changing under reflections. Mathematically, ${\varphi }_l^s$ (r,θ+π/3)=${\varphi }_l^s$ (r,θ) and ${\varphi }_l^s$ (r,-θ)=${\varphi }_l^s$ (r,π-θ)=${\varphi }_l^s$ (r,θ), whereas ${\varphi }_l^p$ (r,θ+π/3)=${\varphi }_l^p$ (r,θ) but ${\varphi }_l^p$ (r,-θ)=${\varphi }_l^p$ (r,π-θ)=-${\varphi }_l^p$ (r,θ). Both ${\varphi }_l^s$ and ${\varphi }_l^p$ are periodic functions of θ, so they can be expanded in Fourier series in cos(6nθ) and sin(6nθ) (n∈N). Reflection symmetry forces ϕ^s (ϕ^p ) to depend on the cosine (sine) terms only. Thus, ${\varphi }_l^s$ (r,θ)=∑ _na_ln (r)cos(6nθ) and ${\varphi }_l^p$ (r,θ)=∑ _nb_ln (r) sin(6nθ).

 

Fig. 1. (a) Schematic representation of a PCF. (b)–(d) Amplitudes of several vortices for increasing values of γ in a PCF with Λ=31µm and a=6µm (λ=1064nm): (b) γ=0,45×10^-3; (c) γ=0,95×10^-3; (d) γ=1,75×10^-3.

Download Full Size | PDF

In the nonlinear case, functions ϕ_l of the form given by Eq. (3) can satisfy self-consistency, according to group theory. The full operator L(|ϕ_l |)=L ₀(n ₀)+L_NL (r^l${\varphi }_l^s$;n ₂) is invariant under the 𝓒 _6v group because L ₀, r^l${\mathit{,}\varphi }_l^s$ and n ₂ are all 𝓒 _6v invariants. Thus L(|ϕ_l |), like in the linear case, provides the representation where ϕ_l lies on. However, group self-consistency is not sufficient to ensure that the resolution of the nonlinear equation (2) with the ansatz (3) provides a nontrivial solution in all cases, since Eq. (2) always admits the ϕ=0 solution. We solve Eq. (2) by means of the Fourier iterative method previously used in Ref. [10] to find soliton solutions in PCF’s. These solutions also fulfill the group self-consistency condition and they belong to the fundamental representation of 𝓒 _6v. The important difference now is that we restrict ourselves to a different representation space; we search for nonlinear solutions in the l=1,2 representations of 𝓒 _6v. The ϕ_l solution and its conjugate ϕ*_l represent a vortex and an anti-vortex soliton of order l, respectively.

 

Fig. 2. Effective index of a family of vortex solutions as a function of γ (solid line). Same for a family of fundamental solitons (dashed line), as in Ref. [10]. The shadow region corresponds to the conduction band, constituted by Bloch modes, of the 2D photonic cladding with Λ=31µm and a=10µm (λ=1064nm).

Download Full Size | PDF

3. Results

We have simulated large-scale PCF’s (in order to prevent silica breakdown [10]) with a lattice period, or pitch, Λ=31µm, and several air-hole radius at a fixed wavelength λ=1064nm (a convenient quasi-continuum source). We define a dimensionless nonlinear parameter as γ=$n_{2\left(\text{silica}\right)}^2$ P/A ₀ (A ₀ is an area parameter: A ₀=π (Λ/2)²; P is the total power carried by the solution) and perform calculations for increasing values of γ. We discover that a family of optical vortex solutions of the type given by the ansatz (3) is found. The amplitudes of several solutions (with l=1) for increasing values of γ are represented in Fig. 1(b)–(d). According to group theory (Eq. (3)), these amplitudes have to be scalar functions, thus showing full invariance under 𝓒 _6v, and they have to vanish at the origin, as seen in Fig. 1(b)–(d). Besides, they become gradually narrower as the nonlinearity increases. It is interesting to plot the effective index of these solutions (n _eff=β/k ₀) versus γ, as shown in Fig. 2. The value of neff increases as γ increases, accordingly to the narrowing of solutions depicted in Fig. 1(b)–(d). For comparison, we have also included the n _eff(γ) curve of the fundamental soliton family for the same PCF structure, as in Ref. [10]. Both lie on the upper forbidden band of the perfectly periodic cladding structure. Here it can be clearly envisaged the role of the central defect. Notice that there is no threshold γ (power) to generate a nonlinear vortex soliton, i.e, there is a continuum of solutions in γ starting from the linear (TE or TM) mode. Unlike in perfectly periodic structures, the defect can eliminate the presence of a threshold power [8]. If we consider structures where the linear vortex-like mode is not present [14], a threshold γ is necessary to generate the nonlinear vortex solution in the forbidden band. By playing with the geometric parameters of the PCF, it is possible to tune this threshold value. In all cases, this value is much lesser than in perfectly periodic structures, which can be important for experimental purposes.

 

Fig. 3. Phase of a vortex with l=1 at r=21µm. We represent both the total phase arg(ϕ_l ) (solid line) and the group phase arg(ϕ_l )-lθ (dashed line).

Download Full Size | PDF

Even more interesting is the phase behavior of these solutions. We show in Fig. 3 a typical phase profile of a PCF vortex calculated at a fixed radius. It shows strict accordance with the group theory prediction for the phase: arg(ϕ_l )=lθ+${\varphi }_l^p$ (r,θ)=lθ+${\sum }_{n=1}^{\infty }$ b_ln (r) sin(6nθ) and it points out that n≈1. The phase of these vortex solutions differs from that of an ordinary vortex by the presence of the pseudo-scalar function extra term: the group phase. Besides the existence of the standard linear behavior in θ, the group phase provides an additional sinusoidal dependence on the angle (with period determined by the group order: π/3 for order 6) not present in ordinary vortices.

In order to check the stability of a vortex soliton solution ϕ_l , we have to consider z-dependent perturbations. This implies solving the non-paraxial equation (1) for the perturbed field ϕ′=ϕ_l +δϕ. In terms of group theory, (ϕ_l,ϕ*_l) constitute an orthogonal basis of the l representation of 𝓒 _6v, so we can consider two types of perturbations around a vortex solution ϕ_l : diagonal perturbations, in which ϕ′ remains in the one-dimensional subspace spawned by ϕ_l , that is, 〈ϕ′|ϕ*_l〉=0, and non-diagonal, where the perturbation takes ϕ′ out of this subspace, 〈ϕ′|ϕ*_l〉=0.

A typical initial condition of the form ϕ′=re^iθ exp[-(r/r ₀)²] is an example of diagonal perturbation for a vortex (l=1). Another example is a perturbation of the form ϕ′=(1+ε)ϕ_l ,ε being a small real constant number. In both cases, 〈ϕ′|ϕ*_l〉=0. We have simulated the evolution of a perturbation of the latter type launched into the fiber core for different PCF configurations and values of ε (ε~0.025-0.1). In Fig. 4 we present an animation that displays the evolution of the transverse field amplitude along the fiber. The animation shows the transient from the initial amplitude towards an asymptotically stationary profile corresponding to a vortex. This transient is a few centimeters long and has the characteristic form of a damped oscillation. Non-paraxial evolution shows that a vortex is stable under diagonal perturbations. Diagonal perturbations generalize the concept of radial perturbations to the discrete symmetry case. We have also simulated non-diagonal perturbations by introducing random perturbations (|δϕ(x,y)|≤ε, ε~0.05-0.1) that modify both the phase and amplitude of the vortex solution in such a way that 〈ϕ′|ϕ*_l〉≠0. Non-diagonal perturbations correspond to azimuthal ones in the radially symmetric case. In this case, an oscillatory instability occurs, as seen in the simulation presented in Fig. 5. This oscillatory instability initially appears as a rotating flux along the vortex ring. After a short evolution, the perturbed vortex breaks into a two strongly interacting single structures of zero vorticity with a cumbersome dynamics, as observed in our simulation. It is remarkable that, unlike vortices in an homogeneous medium [15], these two structures do not fly off tangential to the vortex ring but remain confined in the PCF defect due to the inhibition of transverse flux produced by the PCF cladding. It is also noticeable that, in spite of this instability, no collapse of the structures is detected due to the non-paraxial nature of evolution [16].

 

Fig. 4. (997 KB) Evolution of the field amplitude in z under a diagonal perturbation (ε=0.095) for a large-scale PCF with Λ=31µm, a=10µm and λ=1064nm. We show the transient from an initial profile towards an asymptotic vortex solution.

Download Full Size | PDF

 

Fig. 5. (965 KB) Evolution of the field amplitude in z under a non-diagonal perturbation (ε=0.05) for a large-scale PCF with Λ=31µm, a=10µm and λ=1064nm.

Download Full Size | PDF

4. Conclusions

Summarizing, we propose a novel way of generating vortex solitons, based on PCF’s, in which the presence of the defect plays a crucial role. Our analytical approach, based on group theory, is general and permits to predict the angular dependence of vortex solutions in a nonlinear system owning a discrete invariance, no matter it is periodic or non-periodic. Finally, the use of a non-paraxial equation permits to analyze new phenomena beyond the paraxial approximation.

We thank H. Michinel for useful discussions. This work was financially supported by the Plan Nacional I+D+I (grant TIC2002-04527-C02-02), Ministerio de Ciencia y Tecnología (Spain) and FEDER funds. M. Zacarés acknowledges Fundación Ramón Areces grant.