## Abstract

We present numerical modeling based on a combination of the Bidirectional Beam Propagation Method and Finite Element Method that completely describes the wavelength spectra of point by point femtosecond laser inscribed fiber Bragg gratings, showing excellent agreement with experiment. We have investigated the dependence of different spectral parameters such as insertion loss, all dominant cladding and ghost modes and their shape relative to the position of the fiber Bragg grating in the core of the fiber. Our model is validated by comparing model predictions with experimental data and allows for predictive modeling of the gratings. We expand our analysis to more complicated structures, where we introduce symmetry breaking; this highlights the importance of centered gratings and how maintaining symmetry contributes to the overall spectral quality of the inscribed Bragg gratings. Finally, the numerical modeling is applied to superstructure gratings and a comparison with experimental results reveals a capability for dealing with complex grating structures that can be designed with particular wavelength characteristics.

©2010 Optical Society of America

## 1. Introduction

The development of fiber Bragg grating (FBG) inscription using femtosecond lasers has resulted in two key approaches: the interferometric or phase mask method [1] and point-by-point (PbP) inscription [2–5]. The PbP approach is the focus if this paper as it provides a useful means to produce flexible grating structures, as variations in the grating length, pitch, and spectral response can be readily tailored according to the application. The small spot size (~1μm) permits the development of complicated structures in the fiber core, such as superimposed or sampled gratings [6–9]. Although the PbP method is an innovative and flexible way of modifying the filtering properties of optical fibers, the nature of the non-linear interaction of the light and material unavoidably results in refractive index changes that are confined to the laser beam focus volume. This local confinement of the index change can be a drawback, as a small lateral shift of the focal point can displace the inscribed grating from the fiber axis, resulting in symmetry breaking and changing the grating spectrum in strength and shape. The symmetry breaking or off-axis inscription may be intentional, but it is important to be able to predict how it will modify the FBG spectrum, and whether the grating will be acceptable for any proposed application. Hence, theoretical analysis or numerical modeling can become strong tools to support the inscription process. To date there are no significant studies of the cladding mode spectra arise from femtosecond laser inscribed gratings and this is particularly so for the PbP inscription technique.

Here we apply an advanced version of the bidirectional beam propagation method (BPM) that can be used to simulate wave-guiding structures with multiple reflecting interfaces [10] and combine it with the finite element method (FEM) thereby optimizing the solution convergence and accuracy for the propagating mode calculation while reducing the computation time required. Our analysis predicts the losses induced during inscription and gives prominence to the cladding modes. We confirm that any lateral shift of the grating leads to an exponential reduction in FBG reflectivity, with a concomitant modification of the cladding mode spectrum. Furthermore, as the reflectivity decays with the structural shift from the core center, a ghost mode appears which becomes stronger as the lateral shift increases. We note the strong convergence of the analysis with experimental results allowing us to accurately reproduce the physical parameters of the inscribed grating, such as the spatial extent and strength of the refractive index modulation, the grating length, its displacement from the core center and whether there is any device tilt.

## 2. Simulation method

The BPM is a very efficient, flexible and extensive calculation approach that can be readily applied to complicated structures. The bidirectional BPM algorithm in particular [10] is well suited to gratings in optical fibers as it considers coupled forward and backward traveling waves, readily accounting for multiple reflection phenomenon, including resonant effects as found in grating structures. Furthermore, the approach automatically includes the effects of both guided and radiating fields as well as mode coupling and conversion.

The well known Helmholtz equation under the scalar field assumption for monochromatic waves can be written as [10]

*k(x,y,z) = k*is the spatially dependent wavenumber, and the geometry of the problem is defined by the refractive index distribution

_{0}n(x,y,z)*n(x,y,z)*. To determine the launch field we use the slowly varying envelope approximation (with a slowly varying field

*u*proportional to the phase variation

*ϕ*) where the scalar field

*ϕ(x,y,z)*can be calculated along the propagation direction using [10],

To reduce the computation time and improve accuracy, we pre–calculate the mode supported by the structure using FEM. In this way the solution converges faster as the launch field is closer to the propagation mode. We note that FEM gives more accurate solutions for the supported modes and the effective refractive indices for each wavelength than the Beam Propagation Method [11]. Because of complexity of grating structures especially in the case of superimposed gratings and sample gratings, even small diversions in the initial parameters are crucial for the final solution and the calculated spectra. In order to optimize our final results we decided to combine those methods.

In order to handle simultaneous propagation along the negative z axis (i.e. 180°) the backward travelling waves are treated as a separate, though coupled part of the problem [10]. The guided wave propagation problem is divided into regions that are uniform along the propagation axis z, and the interfaces between these regions. At any point along the structure it is assumed that both forward and backward waves can exist, which are denoted by *u ^{+}(x, y, z)* and

*u*, respectively. In the uniform regions the forward and backward waves are decoupled, while the interfaces between these regions couple the forward and backward waves due to reflection. A transfer matrix approach is employed in which the individual matrices are differential operators. The transfer matrix problem, however, is formulated by assuming that both the forward and backward fields are known at the input of the structure, and an overall transfer matrix

^{−}(x, y, z)*M*describes the system as follows [10]:

The physical propagation problem requires two key pieces of information, the refractive index distribution, *n(x,y,z)*, and the input wave field, *u(x,y,z = 0)*. From these, the physics dictates the wave field throughout the rest of the domain, *u(x,y,z>0)*. Additional input is required in the form of numerical simulation parameters such as the finite computational domain, *{x ∈ (x _{min},x_{max})}*,

*{y ∈ (y*, and

_{min},y_{max})}*{z ∈ (z*, the transverse grid sizes,

_{min},z_{max})}*Δx*and

*Δy*, and the longitudinal step size,

*Δz*. Clearly the calculation accuracy is determined by the computational spatial resolution, the error in the effective refractive index of the supported mode and the step in the wavelength domain.

## 3. Simulation and experimental results – a comparison

As with any numerical method, the modeling has a minimum requirement for inputs and the accuracy of the analysis is very sensitive to them. The required inputs in our numerical analysis are the launch field, which determines the supported mode field distribution and the physical geometry of the fiber, characterized by refractive index distribution. Leaky modes have been taken into account for all our structures as their existence will add to the insertion losses (changes in the background transmission away from the mode resonances) because the measured spectra are normalized to the launch field. This is accounted for using the correlation method [10] to evaluate the effect of leaky modes, whereby an arbitrary field is launched into the structure through the BPM, and correlating the input field and the propagating fields, leads to the evaluation of the modal propagation constants and their fields. For our analysis the structure was set to be single mode fiber (SMF–28).

Figure 1
shows the pre-calculated dominant mode and its effective refractive index. Our analysis is very sensitive to initial parameters and that sensitivity is enhanced as the complexity of the structure increases, hence we are interested in the optimum calculation of the dominant core mode that is supported by the structure, and the effective refractive index for that mode. As a result of our initial study we conclude that FEM has improved accuracy on the effective refractive indices calculation compared to the BPM. Despite the differences remaining below 0.1% (n_{eff} = 1.444602 with FEM and n_{eff} = 1.444772 with BPM), their impact on the final solutions proves important. A comparison between the two computation methods shows a small deviation, with the FEM solutions closer to our experimental results. Finally, there is a significant reduction in the calculation time by using the FEM.

The optimized data pre-calculated by FEM was set as the launch field or the forward propagated field at the entrance to the fiber (*u ^{+}_{in}*) in the matrix problem (Eq. (3), whereas the BPM was used to solve the problem; this takes advantage of the strong points of each calculation method. Figure 2
shows a typical comparison between calculated and measured spectra in a fiber containing a Bragg grating. A first order grating was inscribed in the center of the core of a standard telecom fiber. For the inscription a femtosecond laser system (HighQ Laser Femtoregen IC335) emitting 300 fs pulses at 1035 nm and operating at a repetition rate of 1 kHz was used while the light was focused with a Mitutoyo, long working distance microscope objective ( × 50, NA 0.42). Key parameters of the grating such as the period and the size of the spots have been measured using an optical microscope. Those parameters are required for our analysis. The size of the core was set to be 8.3 μm and the period was 536 nm. The width of each spot was 1 μm and the length of the grating set to be 2 cm. As we can see from Fig. 2, the analysis predicts the exact number of the cladding modes with an average wavelength determination error of 0.05 nm, or 0.005%. Furthermore, the small average deviation between calculated and measured insertion loss (4.6 × 10

^{−3}dB) is indicative of the reliability of the predictions. Small deviations of the analysis compared with the experiment relate to the differences of the spot shape associated with the index change. Our model assumes that the spots are homogenous but in the fiber the refractive index change has an inhomogeneous distribution, because of the nonlinear laser pulse interaction with the material. The modulated refractive index amplitude of the model was set as 4 × 10

^{−4}to be in agreement with the measured value obtained from phase contrast microscopy and the application of the inverse Abel transform. Stronger refractive indices increase the reflectivity of the grating but at the some time give rise to far stronger cladding modes in the spectrum. For sensing applications where the quality of the grating is the dominant factor instead of the reflectivity, cladding mode suppression is very important. In the experimental work all gratings were recorded just above the non-linear threshold value so that we could make distinct optical microscopy measurements in order to unambiguously determine the refractive index value used in the modeling. We carefully controlled the femtosecond laser energy such that we have a near zero mean index change and a net positive modulated index change. This meant that the gratings were not particularly strong, but were of greatest relevance to sensing applications, with good spectral quality.

Is well known that femtosecond laser inscribed Bragg gratings are often accompanied by strong cladding mode resonances thereby reducing the spectral quality and limiting their use for some applications. For example, it would be beneficial for a multiplexed array of optical fiber sensors to have weaker cladding mode contributions to aid wavelength shift recovery. This is achievable by inscribing the grating structure in the central position of the fiber core with a uniform shaped index distribution. The coupling to cladding modes is a result of the lateral shift of the grating in the fiber core (see Fig. 5 , section 4) and the non–uniform shape of the grating spots. For the grating shown in Fig. 2 the cladding mode strength of the experimentally inscribed Bragg grating is quite small, because the grating is very close to the fiber axis (improved alignment techniques) and also because of the fine control of the inscription pulse energy ensuring that the shape of the spots are well defined.

For further investigation of the quality of the simulation we have inscribed gratings with different characteristics, such as longer and non-centered gratings with different spot sizes. Figure 3 presents the spectrum from a grating inscribed within the core, just over 1 μm from the core center. Because of the higher pulse energy used for the inscription (compared to the case of Fig. 2), the size of the spots was larger (~1.5 μm). The greater lateral shift and the bigger spot size explain the strong cladding mode resonances that appear on the shorter wavelength side of the Bragg resonance. Despite the increased complexity in the refractive index distribution of the structure, our model still predicts the transmission spectrum with high accuracy. This shows that the main features of the transmission spectrum can be modeled even for asymmetric structures; of course there is a trade-off between the complexity of the structure and the computation time.

Because of the very good agreement between our analysis and experimental results, we expanded our study to structures with even lower symmetry. Modeling structures with a low degree of symmetry is more challenging and we anticipate that simulation results may diverge more from experimental data. We have simulated superimposed gratings in a single mode fiber and have compared them with experimental results. The inset in Fig. 4 (a) shows the index profile of the superimposed gratings. For this study two similar first order Bragg gratings were placed in the core of the fiber; one is located on the fiber’s axis while the other is placed laterally 3μm from the first grating. The period of the gratings is 0.532 μm and the spot size set to be 1 μm for both of the gratings. Such gratings can only be inscribed with femtosecond lasers because of the non linear interaction of light with matter. Figure 4 (a) shows the simulated spectrum predicting 12 cladding modes between the two grating peaks, which are in agreement with experimental data shown in Fig. 4 (b). The relative strength of the peaks, resonance wavelengths and the shape of the spectra are also in good agreement.

Despite the fact that the structures are now more complicated and less symmetric, deviation between the modeling and experimental results is on the same scale as discussed earlier where there was only one grating inscribed in the core. However, the mesh in FEM must be denser in order to extract sufficiently accurate solutions. Also, we observed that non–uniform meshing close to areas with small index contrast improves the accuracy of the solution. This happens because the accuracy of the underlying simulation algorithms can heavily depend on the correct modeling of field variations that can sometimes be more pronounced near interfaces than in bulk regions. For more complicated structures, computation cells on index interfaces were set to be five time smaller than in the bulk region. Finally, grid cells were set to be aligned in a symmetric way with respect to interfaces between material regions. These results reveal the potential of the analysis and its use for PbP FBG inscription optimization, which is key to the fabrication process.

## 4. Analysis of more complicated structures

Given the previous success we continued with investigations of more complicated and less symmetric structures. Figure 5 shows the calculated transmission spectra for gratings inscribed in different positions of the core. We observed that the lateral shift in the grating position breaks the structural symmetry and a strong “ghost mode” appears close to the Bragg resonance. This is also observed in Fig. 3 where a grating was inscribed that was offset from the fiber axis by approximately one micron.

We believe that the “ghost mode” is a spatial mode supported by the new geometry. With increasing offset the reflectivity of the grating drops whilst the ghost mode increases. The asymmetry perturbs the lower order cladding modes that are closer to the core, rather more than the higher order modes. We also noted that gratings inscribed on the core–cladding interface (labeled boundary in the figure) display a stronger ghost mode with a simultaneous decay in the Bragg resonance reflectivity. The cylindrical fiber symmetry no longer exists with the lateral shift of the grating resulting in more cladding modes and an increase in loss. Finally, the reflectivity reduction of the grating is due to the smaller overlap integral of field and index perturbation (the grating overlaps with the wings of the field profile as it is shifted from the fiber axis).

To reveal the role of the symmetry of the structure in the transmission spectrum a special setup was created. Initially, a grating was placed on fiber’s axis. In order break the structure’s symmetry a second grating was placed in the cladding close to core–cladding interface (see Fig. 6 (a) ). In that way the light that escapes form the core interacts with different refractive index distributions in the opposite directions. Finally, in order to restore the symmetry of the structure a third similar grating was placed in the opposite side of the core to the second grating (see Fig. 6 (b)).

We calculated the transmission spectra for the three different structures focusing on the appearance of the ghost mode and its relationship to asymmetric grating inscription and Fig. 7 presents our results. We observe that asymmetry in the refractive index distribution increases the ghost mode strength in general.

According to our hypothesis, a recovery of the symmetry should erase characteristics introduced in the transmission spectrum by the second grating. The green line spectrum in Fig. 7, which represents the transmission spectrum after the symmetry reconstruction, has similar characteristics with the transmission spectrum of the single, centered grating. The ghost mode disappears while the reflectivity of the grating increases. Moreover, the small red shift that appears when more gratings have been inserted to the structure can be explained by the fact that the additional gratings increase the effective refractive index (the gratings have a positive refractive index change). Finally Fig. 7 shows an increase in the FBG resonance strength, because gratings in the cladding interact with the wings of the field profile adding to the total reflectivity of the structure.

To expand our analysis we calculated the reflectivity of a grating according to the lateral shift from the center of the core. In Fig. 8 we observe that for displacements greater than 1 μm the gradient of the slope increases rapidly leading to lower reflectivity even for small displacements. This is in agreement with experimental observations. Empirical mathematical analysis shows that reflectivity reduction follows a Gaussian function according to the displacement of the grating and is given by

*x*is the displacement of the grating from the center of the core. A Gaussian fit on the calculated values lead to Eq. (4), which is in agreement with our experimental observations noticing that grating’s reflectivity is strongly correlated with the position in the fiber core. Here we have to say that Eq. (4) is valid for a grating offset 0 >

*x*> 5μm to avoid negative reflectivity values.

All previous analyses are based on the assumption that the grating is parallel to the core axis; however, in practice there is always some degree of grating tilt (Fig. 9 (a) ). Figure 9 (b) shows that as the angle increases the strength of the grating drops, ghost modes arise, and for small angles there is no spectral wavelength shift. This tilt also introduces asymmetry but with a lower impact than when the entire grating was laterally displaced by one micron. Transmission spectrum calculations have been made for gratings with different relative rotations from the fiber’s axis. The start of the grating is located at the center of the core and the calculations have been made for four different angles starting with the grating parallel to fiber’s axis and tilted at 1°, 2° and 4°. On careful inspection we note that there is an interesting observation regarding the cladding mode strength. As shown at point A, cladding modes from less axially tilted gratings are stronger than those from parallel gratings, whereas at point B the opposite is true. This provides useful information on the grating position and tilt. Furthermore, phase mismatching appears even for very small rotations and that is translated to a drop in grating reflectivity. This fact reveals the importance of symmetry to the quality of the spectrum.

## 5. Superstructure fiber Bragg gratings

Superstructure fiber gratings (SFGs) are a particular kind of grating with important optical properties for applications in telecommunications and sensing [12,13]. The SFG is a composite grating structure formed when the Bragg grating is exposed to a slowly varying refractive index (amplitude) modulation along the grating length thereby inducing counter-directional mode coupling, resulting in a comb filter response viewed in reflection.

Light penetration along the structure is crucial, so the Bragg grating reflectivity must remain below of 25% to fall within the Fourier design limit. For weak Bragg gratings the reflection spectrum is linearly correlated to the Fourier transform of the refractive index modulation. Here we are going to compare our numerical analysis results with experimental measurements that have been made for typical PbP SFG. In our analysis the Bragg gratings have a period of 539nm, a spatial separation of 600μm for a 50:50 duty cycle, the refractive index modulation is 8 × 10^{−4} and the superstructure’s total length is 2cm. Figure 10 (a)
shows a refractive index profile of the superstructure. The modeling solves Maxwell equations and the reflected signal is calculated.

As we observe from Fig. 10 (b) the experimental measurements and modeling show very strong agreement with regard to the resonances wavelength, reflectivity and peak separation of the resonances. Furthermore, side–band resonances around the central peak and the prediction of peaks with far lower power levels confirm that our modeling is very accurate over a wide range of wavelengths. Despite the agreement between experimental and simulated results there is a difference in the FWHM between spectra, nonetheless this is the first accurate modeling of superstructure gratings via this method, and further refinement is under way. Finally, when modeling these gratings we set the superstructure to be highly symmetric (near the center of fiber core), and we note the excellent agreement with the experimental results, which indicates that the inscribed superstructure is the very close to the fiber axis.

## 6. Conclusion

We have applied both FEM and bidirectional BPM to completely describe the transmission and reflection spectra of PbP femtosecond laser inscribed FBGs, predicting insertion loss and all dominant cladding and ghost modes that are strongly dependant on the FBG position in the core. We have proven that symmetry breaking is responsible for phase mismatching and the transmission spectrum can be affected even by a minute tilt between the grating and fiber axis. Moreover, superstructure gratings have been successfully modeled. Further investigation of superstructure characteristics will follow in the future. This analysis and approach can be used to improve and optimize the inscription of PbP fiber Bragg gratings that are simple or complex in nature.

## Acknowledgement

The authors acknowledge the financial support of the European Commission via project PHOSFOS (PHOtonic Skins For Optical Sensing), and the comments of the reviewers that have improved the quality of the paper.

## References and links

**1. **S. J. Mihailov, C. W. Smelser, P. Lu, R. B. Walker, D. Grobnic, H. Ding, G. Henderson, and J. Unruh, “Fiber Bragg gratings made with a phase mask and 800-nm femtosecond radiation,” Opt. Lett. **28**(12), 995–997 (2003). [CrossRef] [PubMed]

**2. **A. Martinez, M. Dubov, I. Khrushchev, and I. Bennion, “Direct writing of fibre Bragg gratings by femtosecond laser,” Electron. Lett. **40**(19), 1170–1172 (2004). [CrossRef]

**3. **G. D. Marshall, M. Ams, and M. J. Withford, “Point by point femtosecond laser inscription of fibre and waveguide Bragg gratings for photonics device fabrication,” Proc. 2nd Pacific International Conference on Application of Lasers and Optics, 360–362 (2006).

**4. **Y. Lai, K. Zhou, K. Sugden, and I. Bennion, “Point-by-point inscription of first-order fiber Bragg grating for C band applications,” Opt. Express **15**(26), 18318–18325 (2007). [CrossRef] [PubMed]

**5. **T. Geernaert, K. Kalli, C. Koutsides, M. Komodromos, T. Nasilowski, W. Urbanczyk, J. Wojcik, F. Berghmans, and H. Thienpont, “Point-by-point fiber Bragg grating inscription in free-standing step-index and photonic crystal fibers using near-IR femtosecond laser,” Opt. Lett. **35**(10), 1647–1649 (2010). [CrossRef] [PubMed]

**6. **A. Arigiris, M. Konstantaki, A. Ikiades, D. Chronis, P. Florias, K. Kallimani, and G. Pagiatakis, “Fabrication of high-reflectivity superimposed multiple-fiber Bragg gratings with unequal wavelength spacing,” Opt. Lett. **27**(15), 1306 (2002). [CrossRef]

**7. **M. Harumoto, M. Shigehara, and H. Suganuma, “A novel superimposed sampled long-period fiber grating,” in *Bragg Gratings, Photosensitivity, and Poling in Glass Waveguides*, 2001 OSA Technical Digest Series (Optical Society of America, 2001), paper BThC16.

**8. **B. J. Eggleton, P. A. Krug, L. Poladian, and F. Ouellette, “Long periodic superstructure Bragg gratings in optical fibres,” Electron. Lett. **30**(19), 1620–1622 (1994). [CrossRef]

**9. **N. G. R. Broderick and C. M. de Sterke, “Theory of grating superstructures,” Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics **55**(3), 3634–3646 (1997). [CrossRef]

**10. **H. Rao, R. Scarmozzino, and R. M. Osgood, “A bidirectional beam propagation method for multiple dielectric interfaces,” IEEE Photon. Technol. Lett. **11**(7), 830–832 (1999). [CrossRef]

**11. **L. Vincetti, A. Cucinotta, S. Selleri, and M. Zoboli, “Three-dimensional finite-element beam propagation method: assessments and developments,” J. Opt. Soc. Am. A **17**(6), 1124–1131 (2000). [CrossRef]

**12. **J. Chow, G. Town, B. Eggleton, M. Ibsen, K. Sugden, and I. Bennion, “Multiwavelength generation in an erbium-doped fibre laser using in-fibre comb filters,” IEEE Photon. Technol. Lett. **8**(1), 60–62 (1996). [CrossRef]

**13. **F. Ouellette, P. A. Krug, T. Stephens, G. Dhosi, and B. Eggleton, “Broadband and WDM dispersion compensation using chirped sampled fibre Bragg gratings,” Electron. Lett. **31**(11), 899–901 (1995). [CrossRef]