1. Introduction

Photonic band gap (PBG) structures have opened new pathways for harnessing light. The PBG optoelectronic devices are designed to mold the properties of light by manipulating the dispersion, symmetry and amplitude of the field. Nonlinear optics in PBG’s has emerged as a way to make compact and efficient photonic devices. Many applications have been envisioned in one-dimensional systems consisting of multi-layer, dielectric stacks or patterned waveguides. Gap solitons was a term coined by Chen and Mills [1] to describe the hyperbolic-secant shaped stationary structures they observed in numerical simulations of periodic, nonlinear dielectrics. Their results were thoroughly analyzed by a number of researchers using envelope function formalisms; a summary can be found in the review of deSterke and Sipe [2]. Potential device applications include nonlinear optical limiters [3]; optical diodes [4]; a photonic band edge laser [5]; and a high-gain optical parametric amplifier for nonlinear frequency conversion [6].

In this paper we numerically solve the nonlinear, coupled equations for counter-propagating waves in a finite, periodic medium. An initial pulse is launched outside the medium and the ensuing interaction leads to reflections and pulse break up. The transverse and longitudinal modulation instability occurs at a threshold value of the field amplitude. Modulation instability phenomena occur in a variety of nonlinear optical systems. For instance, they are identified in a homogeneous optical media [7–9] and optical fibers [10–12]. Recently Lichinitzer et al. [12] examined the regions of spatio-temporal instabilities in a Kerr medium with a Bragg grating. Their linear stability analysis of an infinite medium predicts the appearance of transverse and longitudinal modulations of the field amplitude. We present a computational method, which is similar in spirit to the method introduced by Scalora and Crenshaw [13] to describe beam propagation with simultaneous forward- and backward-propagating waves.

2. Mode-coupled equations and computational approach

2.1 Mode-coupled Equations

A weakly periodic index profile is assumed in order to simplify the calculations for the electromagnetic modes of the structure. In this paper the nonlinear coefficient is assumed to be homogeneous throughout the medium. The field evolution is described by a set of coupled equations for the forward- and backward propagating waves [10]:

where δ=(k- k₀) is the detuning of the laser wavenumber in the medium, k, from half the wavenumber corresponding to the lattice periodicity, i.e. k₀=π/ʌ, where ʌ is the lattice constant ; κ=ω²Δε/2 k c², is the grating coupling coefficient; η=12πω²χ(3))/2kc² , is the nonlinear coefficient (related to the average Kerr coefficient of the medium), F is a diffraction parameter, and v is the group velocity . Our boundary conditions are applied at opposite ends of the sample:

The function S(x,y,t) is the pulse input to the sample. The linear equations are analytically solved by taking Fourier transforms of the time and transverse coordinates, using exp(-iωt) for the time-harmonic component. The solutions are

where Ω = δ + ω/v-q ²/F, and Δ² =Ω² -κ ². The solutions have the transversal and temporal coordinate dependence, albeit in Fourier transformed variables.

2.2 Numerical Computations for Coupled Mode Equations

Equations (1) with boundary conditions (2) are solved by a variation of the Slowly Varying Envelope in Time (SVEAT) method used for deep gratings [13]. Eqs. (1a) and (1b) are integrated by stepping forward in time; a split step operator form is used. The propagation is separated into a 2×2 diagonal matrix of a linear operator, whose components are

The second operator, V, is a matrix with diagonal components we denote by N

and off-diagonal components, which are denoted by K

Eqs. (1) are rewritten in the following vector form (defining U= (E_f, E_b)^T , where the superscript T denotes the transpose of the vector) is

The solution for the time increment t+Δt is expressed in the split-step operator form as:

The operator L is diagonalized by Fourier transforming the (x,z) variables. The split-step method therefore uses fast Fourier transforms to simplify the linear operator matrices. The central operator in Eq. (7) can be subsequently decomposed into a linear portion and a nonlinear contribution, namely

The decomposition of the equations may be done in several ways. Here Eqs. (7–8) provide a simply implemented algorithm. The operator L is solved by spectral methods; fast Fourier transforms (FFT) diagonalize the spatial derivatives. In the longitudinal direction we use 2048 points. We use only one transverse direction; in this paper the number of FFT points is 128.

Since this method steps the solution forward in time, the first boundary condition in Eq. (2) is modified to an initial condition by displacing the function from the boundary at z=0. The relation between the two conditions is F (x,z,0)=S(x,0, (z-z₀F)v), because the pulse propagates without dispersion in the medium outside the nonlinear medium. We also assume that the boundary is impedance matched, so that in the absence of the coupling κ no reflected wave is generated. The initial pulse has been displaced by z ₀ to a position outside the medium. The accuracy of the code was verified by comparing results with the analytic solutions of the linear equations and by verifying the conservation of energy at each time step.

3. Results

Neglecting transverse coupling, the linear dispersion relation is a solution of Eqs. (3) for the form E_x = A_x exp(iKz), where the subscript is x = f, b. The result is a quadratic equation, which can be expressed it in the form

The detuning parameter takes the role of the fields’ angular frequency. The wavenumber is real as long as the magnitude of δ is larger than κ. A band gap separates the two branches. Tuning the laser close to one band or another determines the characteristics of the wave’s propagation. The transmission curve is determined from the coefficient of F(q,ω) on the right hand side of Eq.(3a) when z=L. In our numerical equations the space, time and field amplitude are scaled so that the following parameter values are used: κ= 1,v = 1 and η = 1. Only real, positive values of the nonlinear coefficient are discussed here and for simplicity we chose F=100. For a medium length L=6.1, the transmission curve is shown in Figure 1.

 

Fig. 1. Transmission curve and input pulse spectrum. The parameters have been scaled as discussed in the text. The frequency is shifted relative to the laser line, which is placed at the first transmission maximum on the high frequency side.

Download Full Size | PDF

The equations of motion are solved by stepping forward in time. The Gaussian-shaped pulse is initially displaced to a position outside the medium:

The initial displacement, z₀, of the pulse is chosen to be several times longer than its pulse width σ_z = 1/8π. This prevents any significant overlap between the medium and the leading edge of the wave. The transverse width was chosen as σ_x = π/64. The z-axis spans 20π and the x-axis span is 8π. In our computations F=100 so that the transverse coupling is weak and does not noticeably affect the pulse tuning at a selected frequency.

Figure 2 is an animation of a low-power pulse amplitude propagating through the sample (A= 0.2). The absolute amplitude of the field at each pixel is represented by a color. The vertical axis displays the transverse coordinate and the horizontal axis is the propagation direction. The evolution is stepped through a time sequence beginning with the pulse to the left of the nonlinear PBG medium. The time increment between frames is Δt=1. Three subfigures displayed are displayed in the animation. The top figure is the input Gaussian field. The initial pulse is displaced in the window and propagates to the right. The nonlinear medium begins at the far left of the top panel and is not shown. The middle figure is the continuation of the forward-propagating pulse as it passes into and through the medium. The boundaries of the medium are denoted by the green boxes positioned to the left side on the top and bottom of the figure. The bottom panel shows the propagation of the backward field in the medium. The green boxes on the right side of the figure denote the boundaries of the medium. The pulse width is shorter than the medium; it is reflected inside the medium as contributes to pulse break up along the longitudinal axis. This behavior can be expected by comparing the pulse spectrum with the transmission curve in Figure 1. The pulse is short enough to sample a large portion of the dispersive spectrum in the medium.

 

Fig. 2. Animation of the pulse evolution through the nonlinear medium with initial detuning δ = 1.12, i.e. the central frequency of the laser is tuned to the first transmission maximum on the high frequency branch. The initial amplitude is A=0.2. Top panel is the pulse amplitude launched outside the medium. The center panel is the forward-propagating pulse amplitude and the bottom panel is the backward propagating pulse amplitude. The numbers in brackets denote the maximum value of the intensity in that plot (196 KB gif animation).

Download Full Size | PDF

A similar animation is shown in Figure 3; however, now the input field amplitude has been increased to A = 0.6. In this case we do not plot the input field displayed in the top panel of Figure 2. The laser is still detuned from the center of the gap to the first transmission maximum on the high frequency side, δ = 1.12. As the field interacts with the medium the pulse begins to break up in both the transverse and longitudinal directions. There is also strong focussing of the pulse as it enters the medium.

 

Figure 3: The pulse amplitude launched in the forward direction is shown in the top panel. The bottom panel is the backward propagating pulse amplitude. See Figure 2 for a discussion of parameters. For this case A=0.6 (216 KB gif animation).

Download Full Size | PDF

 

Figure 4: A three-dimensional representation of the final frame in Fig. 3. The pulse amplitude in the forward direction is shown in the top panel. The bottom panel is the backward propagating pulse amplitude.

Download Full Size | PDF

The front end of the forward-propagating pulse is strongly altered by the passage through the medium, while the backward pulse has a relatively smooth leading edge. The tailing portion of both pulses has strong oscillations in both longitudinal and transverse portions. The pulse spread is increased over the small amplitude case. For higher intensities the pulse break up becomes stronger; the spatial frequencies observed in the transverse direction become higher and the pulse spread is stronger. Increasing diffraction effects by lowering F also leads to a larger spread of the pulse in passage through the medium and greater divergence of the nonlinearly modulated beam in free space.

The break up of the pulses is illustrated in the three-dimensional plots in Fig. 4. The top frame is the forward-propagating pulse amplitude is shown for the last frame in Fig. 3. The extent of the pulse’s transverse spreading is more obvious in this figure as the smooth Gaussian envelope break up into multiple peaks. The longitudinal instability is obvious in the lower frame showing the backward-propagating field amplitude. The backward-propagating pulse in Fig. 3 displays a rib structure that is a characteristic of the longitudinal instability.

The modulation instability is readily apparent in our simulations for A=0.4. But a closer look at the transverse spatial spectra of the forward propagating field, defined by

shows deviations from the initial transverse Gaussian. The spectrum H(q) is plotted for three amplitudes in Figure 5. The pulses have propagated through the medium. The spectrum for A=0.2 in Figure 5 is Gaussian. However, the spectrum for A=0.3 is already deformed from a Gaussian shape; the wings of the initial Gaussian profile acquire a broad shoulder. At threshold the transverse coordinate has relatively low spatial frequencies. Above threshold the spatial frequencies become higher; we find the maximum develops in the nonlinear regime around q=2.2 and moves towards larger values.

 

Fig. 5: The spatial spectrum, H(q), of the pulse after traversing the PBG, as defined in Eq. (11) and the peak is normalized to unity. The three curves correspond to the amplitudes A=0.2 (solid blue line), 0.3 (dash-dotted red line), and 0.4 (dashed black line).

Download Full Size | PDF

The pulse break up depends on the laser tuning. The normal dispersion regime occurs for tuning frequencies below the band gap, so its ability to shape pulses is changed [12]. The animation in Figure 6 shows pulse evolution for the case where the laser detuning from the gap center is δ = -1.12. A large field amplitude is chosen for this example A= 0.7. The pulse does undergo a pulse focussing effect in the medium due to the dominance of the nonlinearity. However, the pulse break up does not occur. For our chosen parameters the absence of modulations instability is expected based on the linear stability analysis of Lichtinitzer et al. [12]. The gain curves are peaked at zero wavenumbers in both longitudinal and transverse directions. The one distinctive difference between the upper and lower branches is the opposite signs of their dispersion. The lower branch dispersion is normal, while the upper branch is anomalous. In the normal dispersion regime simulations run to amplitudes of unity have not resulted in pulse break up.

 

Fig. 6. Top panel is the pulse amplitude in the forward-propagating direction and the bottom panel is the backward propagating pulse amplitude. The parameters are the same as in discussed earlier, except δ = -1.12 and A= 0.7 (156 KB gif animation).

Download Full Size | PDF

4. Conclusions

We have presented a computational algorithm to numerically solve simultaneous forward- and backward-propagating, time-dependent fields in a nonlinear media. In a periodic medium the modulation instability is manifested in the transverse and longitudinal directions. The modulation instabilities lead to a rapid spread of the pulse in the transverse direction and temporal break up of the pulse at threshold intensity. Our threshold modulation instability amplitude is around A=0.3. This compares favorably with the value A=0.14 we deduce from the paper of Lichinitzer et al. [12] by using the 1/e widths of the initial Gaussian profiles in Fourier transformed variables, i.e. transverse wavenumber, k_t=16/π ², and longitudinal wavenumber, k_p = 0.4. This comparison is only qualitative though, since the stability analysis is made for an infinite medium. In a finite medium we chose the central frequency at the first transmission maximum, because the field is locally enhanced. Our pulses are short and the spectral width of the input pulse samples a spread of dispersion values. Also the modulation instability gain needs a sufficient path length to grow a local perturbation and our sample was kept short. A stability analysis for the finite system is required for a quantitative comparison with our simulations, but the problem may not be manageable by analytic techniques.

The instability in the transverse direction can be applied to designing optical limiters. On the high frequency branch the low intensity transmission is about 0.4 when tuned to the first transmission maximum; it decreases to 0.19 for A=0.7. On the high frequency branch and tuned to the first transmission maximum the transmission coefficient increases with the input intensity from about 0.4 at low intensity to about 0.6 for A=0.6; however, the diffraction reduces the on-axis energy, so that by using an apertured detector the off-axis portion is eliminated from the system. To put values on the threshold intensity for the modulation instability consider a sample made from two dielectrics with dielectric contrast, 0.01 and a nonlinear optical coefficient n₂= 10^-10 cm²/W, the corresponding threshold intensity in finite systems is about 3 MW/cm². At threshold the nonlinear change of the index is about 0.0003.