1. Introduction

The existence of self-trapped light beam solutions, or nonlinear localized guided waves, in homogeneous media is a classical result in nonlinear optics [1, 2]. These rotationally symmetric solutions are, however, unstable and only exist at a critical power. More recently, it has been proven that stable modes can be localized in two-dimensional (2D) nonlinear photonic crystals constituted by periodic arrays of poles of nonlinear material both for infinite-sized [3, 4] and finite-sized structures [5]. These solutions are non propagating in the direction of invariance of the photonic crystal (axial direction) and, physically, they correspond to the situation of in-plane illumination. Guided localized nonlinear solutions have also been found in 2D linear photonic crystals of dielectric poles with embedded nonlinear arrays [6] (also for in-plane illumination) and, experimentally, in 1D waveguide arrays, in the form of the so-called discrete solitons propagating along the waveguide direction [7]. Modeling has been performed using an effective “quasi-free photon approximation”[3] (valid for lowindex contrasts) or a discretized version of the scalar nonlinear equation [4].

In a different context, much attention has recently been paid to a different type of photonic crystal structures. Photonic crystal fibers are silica fibers that have a regular array of air-holes extending along the entire fiber length (Fig. 1(a)). Guidance is produced by transverse localization of light at a defect (the core region) due to a complex mechanism of interference with the periodic air-hole array that form the photonic crystal cladding. This mechanism is analogous to that occurring in electron crystals in the presence of donor or acceptor impurities [8]. Unlike in dielectric pole lattices, most power is not confined in lattice sites (air-holes) but in silica. Since air has negligible nonlinear response, the possibility of non-linearly confining light in PCF’s has to be restricted to the silica region as well. This excludes the possibility of generating nonlinear localized solutions of the discrete-soliton type, with amplitudes supported mainly on lattice sites. Moreover, the typical air-dielectric filling fraction in these structures is very large as compared to the relative small fractions achievable in some PCF’s. From this point of view, nonlinear PCF’s are closer to a continuous medium (with discrete symmetry) and, therefore, different properties are expected to be satisfied by their solutions. For the same reason, the modeling technique cannot be based on a discrete nonlinear Schrödinger equation scheme. “Quasi-free photon approximation” techniques are not suitable either because PCF’s are not low-index contrast structures.

2. Description of the method

We consider the guidance problem, at a given frequency, in a triangular PCF which can have a nonlinear response in silica. That is, we solve the following nonlinear equation:

${\nabla }_t^2$ being the 2D-transverse Laplacian operator and k ₀ the vacuum wavenumber. We search for electric field solutions with well-defined constant polarization E(x⃗)=uϕ(x⃗). The linear refractive index profile function n ₀(x⃗) is one in the air-holes and equals n _silica in silica, whereas the nonlinear index profile function n ₂(x⃗) is different from zero only in silica ($n_{2\left(\text{silica}\right)}^2$≡3${\mathrm{\chi }}_{\left(\text{silica}\right)}^{\left(3\right)}$/2ε₀ cn _0(silica)). We approach this problem using a generalization of our modal method to describe linear propagation in PCF’s [9]. In our previous method, the linear differential equation for guided modes is transformed into an algebraic problem consisting in the diagonalization of a matrix representation of the differential operator. Its eigenvalues and eigenvectors provide the propagation constants and amplitudes of the guided modes, respectively. This method utilizes periodic boundary conditions. We implement them by putting the system into a finite two-dimensional rhomboid unit cell (Fig. 1(b)) and requiring that the fields fulfill periodic boundary conditions in the directions defined by the rhomboid sides. Periodic boundary conditions considerably simplify the algorithm since they permit the calculation of the matrix elements in an analytic manner [10]. This property converts this method in a precise and versatile tool. A similar simplification also occurs when the inhomogeneous nonlinear term is included. In the latter case, the equivalent algebraic problem becomes a nonlinear eigenvalue equation problem characterized by a second order tensor (the matrix representation of the linear part) and a fourth-order tensor (modal representation of the nonlinear index profile function $n_2^2$(x⃗)). Despite the complicated PCF geometry, it is also possible here to evaluate the fourth-order tensor analytically. This nonlinear eigenvalue equation is formally identical to that found for a soliton solution in a modal approach to dispersion-management systems [11]. By taking advantage of this fact, we use the same iterative projection technique reported in Ref.[11] to find self-consistent solutions of the modal representation of the nonlinear equation (1). Unlike other techniques, no further approximation beyond the modal expansion is needed to deal with the inhomogeneous nonlinearity.

 

Fig. 1. (a) Schematic representation of the transverse section of a PCF. (b) The box with dimensions D×D corresponds to the unit cell used to implement periodic boundary conditions. In our simulations we have chosen D=7Λ, Λ being the spatial periodicity, or pitch, of the photonic crystal cladding.

Download Full Size | PDF

3. Results

We take normalized field amplitudes (∫ ϕ*ϕ/A ₀=1), which is equivalent to define a power-dependent dimensionless nonlinear coupling γ=${\mathit{\text{Pn}}}_{2\left(\text{silica}\right)}^2$/A ₀, where P is the total power carried by the unnormalized field and A ₀ is a magnitude with dimensions of area characterizing the core size (here we choose A ₀=π(Λ/2)², Λ being the spatial periodicity—or pitch—of the PCF). Solutions of the nonlinear eigenvalue equation are found for different values of the PCF parameters and the nonlinear coupling. In all cases, the solution we find correspond to the eigenvector with highest β² of the nonlinear eigenvalue problem. The solution can be envisaged as the fundamental mode of the effective fiber generated by the combined effect of the PCF refractive index and the nonlinear index induced by the solution amplitude itself (n ²(x⃗)=$n_0^2$(x⃗)+$n_2^2$(x⃗)|ϕ_sol(x⃗)|²). Since the fundamental mode of the linear theory is also a localized solution, it is convenient to provide an efficient criterion to establish a distinction between localized nonlinear solutions and the fundamental mode of a linear PCF. We use as a measure of distinction the difference, or gap, (normalized to the given vacuum wavenumber k ₀) Δ≡(β_sol-β_fund)/k ₀, where β_sol and β_fund are the propagation constants of the solution of the nonlinear equation and the fundamental mode of the linear fiber, respectively. The generation of a gap implies that the nonlinear solution has a different shape than the fundamental linear mode.

Although Eq. (1) admits solutions for ordinary-scale PCF’s, their experimental feasibility is uncertain because of the high intensities involved when focusing light in such small cores (a typical core radius of a conventional PCF is approximately 1 µm). These high intensities can produce a silica breakdown when they overcome a critical value. In terms of our dimensionless nonlinear coupling, this critical intensity implies the existence of a maximum value of γ above which the fiber collapses (γ_max≈10^-3). For these reasons, it is more realistic to consider PCF structures that possess larger cores. One simple way of generating such structures is by magnifying the conventional PCF parameters, i.e., the hole radius a and the pitch Λ, (a→Ma, ΛMΛ). The core size area A ₀ is then scaled by a M ² factor and γ by 1/M ². In our simulation, we choose M=10 with respect to an ordinary PCF configuration characterized by Λ=2.3 µm and a ranging from 0.2 to 1.0 microns. That is, we simulate large-scale PCF’s with Λ=23 µm and a=2–10 µm. In this way, we guarantee that all solutions we find will lie below the breakdown intensity threshold.

Using these values, our simulation shows the existence of nonlinear guided solutions in all PCF configurations even for small nonlinear couplings. The shape of the soliton solutions change with γ in the form presented in Fig. 2(a)-(c). As the nonlinear coupling increases, the self-consistent solution is gradually narrower. In terms of the gap function Δ, this fact is reflected in a monotonously increasing behavior of the generated gap with γ, as shown in Fig. 3(a) (solid curves). However, the nonlinear coupling cannot be increased arbitrarily. There is a given value of γ (not to be confused with the maximum γ for silica breakdown), which can depend on a, at which the soliton acquires a critical shape and propagation constant (Fig. 2(d)). At the critical coupling, our calculations recover the rotationally invariant nonlinear solution obtained in an homogeneous medium (also known as Townes soliton) [1]. Above this critical value, our method finds infinitely narrowso lutions (within numerical precision) characterized by an infinite eigenvalue β² (Fig. 2(e)). The physical interpretation of these solutions is clear. As we increase the nonlinear coupling, the localized mode is more confined and it feels the photonic crystal cladding less and less. At the critical coupling, it stops seeing the cladding completely. Thus, its behavior corresponds to a nonlinear mode in an homogeneous medium. That is why we recover the rotationally invariant solution (no discrete symmetry is left) of an homogeneous medium. This argument also applies above the critical coupling, where we expect to recover the same physics as in an homogeneous medium. It is known that the Townes soliton in an homogeneous medium is unstable under power perturbations. In particular, if power is increased the solution collapses experimenting a process of filamentation induced by a self-focusing instability and generating an infinitely narrow solution with an infinite β² [12]. This is the same solution our method detects above the critical coupling (Fig. 2(e)).

It is interesting to notice that we have also simulated with our method the limit where all holes disappear (a→0), which corresponds to the homogeneous medium case. The differences and similarities with respect to the PCF case are enlightening. Below the critical coupling, the method does not detect any nontrivial solution and only finds the trivial zero solution of Eq. (1). This is also consistent with known results for an homogeneous medium, since no stationary nonlinear solutions are found below the critical power that corresponds to the Townes soliton [1, 12]. At and above the critical coupling, the Townes soliton profile and the infinitely narrow solution are recovered, respectively.

 

Fig. 2. Intensity distribution of different solutions of Eq. (1) for a PCF with pitch Λ=23µm, radius a=4µm, λ=1.55 µm and different nonlinear couplings: (a), linear mode (γ=0); (b) and (c), spatial soliton solutions for γ=0.0010 and γ=0.0015, respectively;(d), unstable nonlinear solution of the homogeneous medium (γ=γ_c=0.0017); and (e), self-focusing instability (γ>γ_c).

Download Full Size | PDF

 

Fig. 3. (a) Dependence of the gap function Δ on the nonlinear coupling γ for different hole sizes. As in Fig. 2, Λ=23µm and λ=1.55µm. (b) Diagram of existence of solutions for a nonlinear PCF. The shaded region is the nonlinear soliton phase. The other region corresponds to the homogeneous-instability phase. The inter-phase is given by the γ_c(a) curve. As before, Λ=23 µm and λ=1.55 µm.

Download Full Size | PDF

The soliton formation pattern depends on the interplay between the PCF geometry and the nonlinear coupling. It can be more easily visualized using the diagram of existence of solutions in Fig. 3(b). A systematic study of different soliton solutions of Eq. (1) for different values of the hole-radius a and the nonlinear coupling γ allows to establish two different regions, or phases, separated by their corresponding inter-phase curve. Since nonlinearities can dynamically induce a propagation constant gap with respect to the linear regime, the phase diagram in Fig. 3(b) is constructed according to the value of the gap function Δ. A finite non-zero gap characterizes the region where spatial soliton solutions with finite propagation constant exist. The region above the critical coupling corresponds to infinitely narrowso lutions characterized by an infinite gap. As mentioned previously, they are the result of a self-focusing instability acting on the Townes soliton of the homogeneous medium when the critical power is overcome. This fact justifies to call the white area in Fig. 3(b) as the self-focusing instability region. The curve that describes the inter-phase between the two regions is given by the function γ_c(a). Solutions lying right on the interface curve γ_c(a) constitute the family of Townes solitons of the homogeneous medium characterized by different beam diameters. PCF’s with larger holes radii a generate more localized solutions and, consequently, also give rise to Townes solitons with smaller beam diameters. Their propagation constant is thus larger and so the gap function is, as reflected in Fig. 3(a). Notice that our numerical simulation provides a nearly vertical inter-phase curve γ_c(a). This is consistent with known results for an homogeneous medium, where all solutions of the Townes soliton family are characterized by the same critical power, independently of the beam diameter [12, 13]. In terms of our dimensionless nonlinear coupling constant, γ=${\mathit{\text{Pn}}}_{2\left(\text{silica}\right)}^2$/A ₀, solutions lying right on the interface curve must share an identical critical γ_c(a). It is also interesting to point out that in the homogeneous case the diagram of existence would be reduced to the vertical inter-phase curve γ_c(a) exclusively.

 

Fig. 4. (2.55 MB) Evolution of the field amplitude in z for a large-scale PCF with Λ=23µm, a=8µm and λ=1.55µm. We show the transient from an initial Gaussian profile towards an asymptotic spatial soliton solution.

Download Full Size | PDF

Up to now, we have been only concerned with the existence of nonlinear solutions. We have found solutions of the nonlinear eigenvalue equation (1) but their stability has not been proven yet. We can only claim that nonlinear solutions lying on the inter-phase curve γ_c(a) are unstable since we recognized them as Townes solitons of an homogeneous medium and instability is a known feature of them. Stability is a dynamical issue which implies the study of the evolution equation associated to Eq. (1):

We have checked that nonlinear solutions corresponding to the shaded region in Fig. 3(b) are stable under small arbitrary perturbations. In order to do so, we simulate the evolution of the perturbed field using Eq. (2) to show that a nonlinear solution of Eq. (1) appears as its asymptotic state. Numerically, we use a modal beam propagation method to implement the evolution Eq. (2). We highlight here the main features of our results. Details of both the method and the stability analysis will be provided elsewhere. We have simulated the evolution of a spatial Gaussian profile launched into the fiber core for different PCF configurations. 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 Gaussian amplitude towards an asymptotically stationary profile corresponding to one of the nonlinear solutions of the diagram of existence in Fig. 3(b). A typical mechanism of dynamical stability is apparent. The initial profile re-shapes as it sheds energy in the form of outward dispersive waves until it reaches the shape of a stationary non-radiating solution. This transient is a fewcen timeters long and has the characteristic form of a damped oscillator. In order to visualize this mechanism more quantitatively we represent the evolution along the fiber axis z of the expectation value of the differential operator in Eq.(1)—let us denote it by L(ϕ)—with respect the field amplitude ϕ(x, y, z). This field verifies the 3D evolution equation (2), which in terms of the L(ϕ)-operator is expressed as L(ϕ)ϕ=-∂²ϕ/∂z ². Thus, the expectation value of the L(ϕ)-operator, 〈ϕ|L|ϕ〉, is z-dependent unless it represents an stationary solution of the nonlinear eigenvalue equation (1): ϕ(x, y, z)=ϕ_s(x, y)e ^iβz.

 

Fig. 5. Typical evolution behavior in z of 〈ϕ|L(ϕ)|ϕ〉 for a large-scale PCF with Λ=23µm, a=8µm and λ=1.55µm.

Download Full Size | PDF

In such a case, 〈ϕ|L|ϕ〉=β². In Fig. 5 we show a typical behavior of the evolution of 〈ϕ|L|ϕ〉 (z) where the stabilization mechanism becomes evident: the expectation value evolves (oscillating) until it asymptotically reaches a plateau indicating the existence of a spatial soliton solution. The asymptotic β² value corresponds to the spatial soliton solution of Eq. (1) for the same PCF parameters with a γ given by the asymptotic power. Notice that power is not conserved because of radiation and decreases during propagation, so that the asymptotic solution possesses a smaller power than any other intermediate solution. It is worth mentioning that we have also checked that these spatial solitons are stable under both small transverse displacements relative to the photonic crystal cladding and launching with small transverse momentum (slight offaxis illumination). Analysis of the evolution of the expectation value 〈ϕ|L|ϕ〉 (z) shows in both cases a similar asymptotic behavior as for the Gaussian profile in Fig. 5: damped oscillations tending towards a spatial soliton plateau.

4. Conclusions

We have demonstrated that large-scale PCF’s can support stable nonlinear localized solutions. They constitute a newclass of spatial soliton solutions characterized by a discrete symmetry and different from the so-called discrete spatial solitons. These new spatial solitons in PCF’s distinguish from discrete spatial solitons in the fact that they cannot be described by the discretized nonlinear Schrödinger equation. Despite both enjoy the sixth-fold symmetry of the photonic crystal, their physical properties are substantially different. We have shown that spatial solitons in PCF’s are closer to nonlinear solutions in a continuous medium. This feature is apparent in the fact that the unstable rotationally invariant solution of an homogeneous medium is recovered for specific large nonlinear couplings and also in the appearance of identical self-focusing instabilities.

On the other hand, our results point out the crucial role played by the photonic crystal cladding in the stabilization mechanism of these newso lutions. Besides the ordinary interplay between diffraction and nonlinearity, the photonic crystal cladding provides an extra and peculiar localization mechanism. In this way, even for small effective nonlinear couplings (or equivalently, for small input powers) stable nonlinear solutions can be generated. Without the existence of the additional confinement provided by the photonic crystal cladding, nonlinearities would not be able to build such solutions in this regime, as known in the homogeneous case. Remarkably enough, the photonic crystal structure does not only help to generate these solutions but it also gives rise to their stabilization.

Finally, we present a proposal that overcomes at first instance the problems derived from collapsing the fiber due to high intensities. Large-scale PCF’s prevent intensities from reaching values above the threshold silica breakdown.

We are thankful to Alex Gaeta 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 gratefully acknowledges Fundación Ramón Areces grant. D. Binosi gratefully acknowledges a grant, BFM2001-0262, from Ministerio de Educación y Cultura, Spain.