Optimisation of microstructured waveguides in z-cut LiNbO3 crystals

We present a practical approach to the numerical optimisation of the guiding properties of buried microstructured waveguides, which can be fabricated in a z-cut lithium niobate (LiNbO3) crystal by the method of direct femtosecond laser inscription. We demonstrate the possibility to extend the spectral range of low-loss operation of the waveguide into the mid-infrared region beyond 3{\mu}m.


Introduction
Tight focusing of intense radiation of femtosecond (fs) laser inside transparent dielectrics leads to an alteration of the material properties, e.g., the refractive index (RI), around the focal region only. This enables unprecedented flexibility in the design of waveguide (WG) geometries [1] and the fabrication of various integrated optics components and light-guiding circuits [2]. In crystals, by writing multiple tracks with a slightly reduced RI around the unmodified volume of material it is possible to produce a depressed-index cladding [3,4] with the central volume serving as the core of a WG. Tunable, narrow-and broadband light sources and detectors in the mid-infrared (IR) wavelength range (2)(3)(4)(5) µm) are the focus of intensive research because of a wide range of enabling applications [5], including ones in remote sensing, spectroscopy, biomedical sensing, imaging, and others [6][7][8][9]. After the advent of quasi-phase matching by poling of ferroelectric crystals [10], lithium niobate (LiNbO 3 ) has become one of the most commonly used materials in devices such as acousto-optical filters, frequency converters and optical parametric generators. Microstructured WGs with low losses, high damage threshold, and controlled dispersion in a broad spectral range are highly demanded for both traditional applications and applications in recently emerged fields, such as integrated quantum optics and mid-IR range frequency comb generation. Therefore, it is of interest to establish appropriate design principles and develop suitable fabrication technologies. Recently, depressed-cladding WGs have been fabricated by fs laser direct writing into bulk optical glasses (isotropic materials) [11,12]. Note that the numerical model used in these works builds upon a scalar diffraction approximation and, thus, cannot be used to describe vector propagation as required by birefringent hosts such as LiNbO 3 crystals.
In a previous work [13] we have demonstrated numerically that the guiding properties of mi-crostructured WGs in z-cut LiNbO 3 crystals can be controlled by the WG geometry. Obviously, with virtually an infinite number of structural parameters the design of such structures is not as trivial. Our study ranged over the parameter space: track size (D), period or pitch (a), number of cladding layers (N clad ), and RI contrast (δ n), that is accessible experimentally. In particular, the number of track rings revealed to play a major role in the control of the confinement losses, i.e., the losses due to the finite transverse extent of the confining structure. We observed that for the typical induced RI contrast δ n = −0.01, increasing the number of equidistant cladding layers of identical tracks from two to seven in a microstructured WG with an hexagonal symmetry can result in a reduction of the losses in both ordinary (O) and extraordinary (E) polarisations by more than three orders of magnitude at the telecommunication wavelengths. Clearly, the range of possible structural parameters that can be investigated for further optimisation is far from being fully explored. We have also shown [13] that WG designs with track diameters that differ from one ring to another [14,15] can extend the spectral region of low-loss operation to longer wavelengths. However, it emerged from our study that a systematic procedure is required to find optimum laws for the variation of WG parameters such as the track size among different rings of tracks. Another motivation for the present work is that the RI contrast and the track size are not independent parameters -if the track size varies the RI contrast changes too, thus a meaningful optimisation should also account for this dependence. Additionally, the direct fs laser inscription always generates losses both within and around the tracks [16][17][18], whose effect should be properly accounted for. All the fore-mentioned effects become more significant at the long wavelengths, where we observed an unusual variation of the confinement losses in the O and E polarisations [13]. In order to confirm the accuracy of our numerical model we simulated known microstructured optical fibre (MOF) structures [15] using the same material dispersion relations. Quantitative agreement for both the dispersion profile and confinement losses was obtained within a broad wavelength range.
In this paper, we present a practical approach to the numerical optimisation of the guiding properties of depressed-cladding, buried WGs that can be formed in a z-cut LiNbO 3 crystal by fs laser writing. The approach accounts for both a suitable variation of the track size among different track layers, the relationship between track size and induced RI contrast, and the intrinsic losses due to fs laser inscription. For our modeling, we used the COMSOL simulation software based on the finite element method (FEM), which enables an accurate description of both anisotropic materials and materials with absorption (indicated by the imaginary part of the RI). To truncate computational grids for simulating Maxwell's equations and, thus, minimize the effect of boundary reflections, we used a rather wide perfectly matched layer absorber surrounding the cross-section of the WG structure [13]. However, this limits the number of modes that can be followed simultaneously over a wide wavelength range. In order to explore the behaviour of a larger number of modes, we compared our method with the plane-wave method (PWM), which has already been used for the analysis of anisotropic waveguides [19] and MOFs [20][21][22]. The PWM can be computationally much faster than the FEM, and it is also well suited for the study of mode interactions as well dispersion characteristics changes under external influences. However, in its current formulation, the PWM cannot provide an accurate description of leaky WGs, for which the imaginary parts of the effective mode indices are large.

Numerical model for FEM calculations
For this study, we modeled microstructured WGs with a hexagonal structure. The depressed cladding was formed by up to seven rings of cylindrical tracks whose centers were arranged hexagonally, as shown in Fig. 1. These tracks can be written in the bulk of a z-cut LiNbO 3 crystal by direct fs laser irradiation using a transverse inscription geometry. We found experimentally (the results will be published elsewhere) that the track diameter D (in [µm]) and the induced RI contrast δ n depend on the laser pulse energy E (in [nJ]) as follows: where E th = 36.45 nJ is the energy threshold of the 100× micro-objective used for inscription (numerical aperture NA = 1.25). These relations were found to be valid for a laser repetition rate of 11 MHz and the optimum (sample translation) inscription speed of 10 or 20 mm/s (determined by a trade-off between inscription depth and laser pulse energy), up to a maximum available laser pulse energy of approximately 75 nJ. It follows that δ n can be at most −0.015 or −0.02. Note that at laser pulse energies approaching the maximum value the smoothness of the inscribed tracks may be compromised due to a strong self-action of the circular laser beam [23]. In our inscription experiments in LiNbO 3 (which will be described in detail in a future work), RI contrast reconstruction was performed by quantitative phase microscopy (QPM) at the optical microscope. The QPM reconstruction method relies on processing "z-stack" images of each individual track by using both commercial and in-house software. This procedure was calibrated by measuring standard optical fibers with known specifications. Further, the experimentally observed intensity distributions of various WG modes were found to be in good quantitative agreement with the results of numerical modeling. As mentioned in the introduction, a natural strategy to extend the spectral range of low-loss operation of the WG structure, is to allow the track diameter to differ from one ring to another with the exterior rings that have large tracks. The rate of growth of the track size from the innermost to the outermost ring can be parameterized with a single parameter p > 0, so that the track diameter D n in the n-th ring is: where D max and D min are the respective maximum and minimum diameters. Examples of how the growth rate parameter p changes the track diameter in a seven-ring WG structure and the cross-sections of structures for different values of p are given in Fig. 2.
Maxwell's equations were solved to find out the complex effective RIs n eff o,e of the modes of the structure for the O and E polarisation states. The confinement losses L (in decibels per centimeter) were computed from the imaginary part of effective RI I (n eff ) through the formula L = 40πI (n eff ) × 10 4 /[λ ln (10)], where λ is given in micrometers. Details of the numerical computation can be found in [13].

Eigenmodes of anisotropic microstructured WGs by PWM
Maxwell's equations for monochromatic optical waves of angular frequency ω propagating along the x axis of the coordinate system in a transparent, uniaxial crystal with the y optical axis and of permittivityε yield the following wave equation for the magnetic field vector H: where K = ω/c = 2π/λ is the free-space wave number, and Here, n o,e are the O and E RIs of the (intact domain of) crystal, and the distance-and timedependence of the electromagnetic field: exp(iβ x − iωt) is assumed, where β = Kn eff is the propagation constant of the mode propagating in the structure. For a crystal hosting a depressedcladding WG with the track RI contrast δ n < 0, we can use the approximation introduced in [24] and express the RIs as: To apply PW decomposition, the WG structure should be put in a box and tiled periodically along the y-and z-axes. The box should not be too large in order not to include a large area of unmodified material around the structure and, at the same time, should not be too small in order to avoid overlapping of modes in adjacent computational cells. Thus, the size of the box has to be a parameter for the calculations.
The second equation in Eq. (3) can be used to exclude H x (y, z) from the first equation and, thus, obtain: Equations (4) with unknown complex amplitudes H y,z m,n , and G y m,n = 2πm/a, G z m,n = 2πn/a. Substituting these expansions into Eq. (4) and equalising the coefficients of the same exponential factors one can obtain the following linear equations for the modal amplitudes H y m,n , H z m,n : where From Eq. (6), one can determine β 2 and H y,z m,n . The transverse components of the electric field E y,z and the longitudinal component of the Poynting vector S x will be then given by The solution of Eq. (6) involves the calculation of Fourier harmonics of the mode profiles. For arbitrary shapes of the tracks of the WG structure, numerical methods have to be employed to calculate the Fourier coefficients, including, when necessary, digital image processing of real structure cross-sections [25]. In the special case of cylindrical tracks (circular cross-sections), however, analytical expressions for the Fourier coefficients can be found. Let X denote any of η o,e , 1/η o,e , with X(y, z) = X 1 if (y, z) ∈ tracks and X(y, z) = X 2 if (y, z) ∈ otherwise. For cylindrical tracks of radii r i , one can obtain: , m, n = 0 (8) Here S is the area of the surrounding box, J 1 is the Bessel function of the first kind, and y i , z i are the coordinates of the center of the ith track.

Numerical results by PWM and comparison with FEM
For practical applications, microstructured WG structures built of low-RI-contrast tracks should feature low confinement losses and, thus, need to have a fairly large number of cladding layers [13]. For these calculations, we considered seven-ring structures with the pitch a = 2µm, the track radius r = 0.8µm and the track RI contrast δ n = −0.01. The Fourier coefficients and matrices U o,e , P o,e , etc. were calculated for values of the indices m, n up to ±35, and 578 PWs were used. These PWs included mostly modes propagating along the x-axis. The eigenvalue problem was solved numerically using the "cg.f" program from the "EISPACK" package (NetLib). Figure 3 shows the computed modes of the structure. We can note an overlap between the O and E polarisations at the wavelength 1.5 µm. This overlap will lead to additional losses in the E mode due to a perturbation induced in the WG. In Fig. 4, the effective RI profiles for the fundamental mode are compared to those obtained using the FEM. The PWM and FEM results are in excellent agreement for the O polarisation, whereas a small deviation within the wavelength range 0.5 µm to 1.5 µm can be observed for the E polarisation. This is likely due to a leakage of the E mode with lower RI into the higher-RI O polarisation [13], which is not accounted for by PWM calculations. The effective mode index values at different wavelengths as obtained from PWM and FEM calculations are given in Table 1.

Optimisation of WG structural parameters by FEM simulations
To minimise the confinement losses in the WG structure at the long wavelengths and, thus, extend the spectral range where the loss figures for the modes are acceptably low, we performed FEM simulations. This is because the PWM in its current formulation does not allow to calculate the imaginary parts of the effective mode indices. To find the leakage losses with the PWM, one should compute the localization coefficients for the modes or, alternatively, include an absorbing medium at the boundary of the computational supercell. The first approach would require the "calibration" of the losses for given localization coefficients, while the second approach is similar to using a PML absorber as in the FEM. We assume here 1 dB/cm to be an acceptable loss level for technological applications. A simple idea for addressing the numerical optimisation problem described in this section came from the observation that the confinement losses become monotonic functions of wavelength at sufficiently long wavelengths. Thus, one can vary different WG structural parameters at a fixed wavelength of interest, and only when the best loss figures are obtained, perform the full wavelength scan. This makes the optimisation procedure much less time consuming and, thus, practically feasible. The FEM results presented here refer to the fundamental mode of the structure, which was selected using the criterium of 'minimum effective mode area' during the wavelength scanning [13]. The RI contrast induced by fs inscription is the most important parameter for mid-IR operation (and successful optimisation) of the WG. Figure 5 illustrates the dependence of the confinement loss on the track diameter at the wavelength λ = 1.55 µm for a seven-ring, uniform (p = 0) structure with different values of δ n. It is seen that for δ n = −0.01 and δ n = −0.02 there is a "plateau" of low losses  The influence of a varying track diameter among the different cladding layers on the loss properties of the WG at the wavelengths 1.55 µm and 3 µm is illustrated by Fig. 6. Our simulations showed that the most interesting range of growth rate parameter values is 0 < p ≤ 1. We can see from Fig. 6 that as p decreases from 1, the confinement loss also decreases up to much lower values than that of a uniform structure (corresponding to p = 0 when D n = D max ∀n). At λ = 1.55 µm, the loss profile is flat over the p range p → 0 (yielding D 1 = D min and D n ≈ D max , n > 1) to p = 0.5 for the maximum track diameter D max = 2.4 µm. For a uniform structure D max = 2.2 µm yields lower confinement loss than D max = 2.4 µm. On the other hand, the loss for a structure with p = 1 is higher than that for the uniform structure at D max = 2.2 µm. Differently, at λ = 3 µm the decrease of the confinement loss with decreasing p values is approximately linear for both D max = 2.2 µm and D max = 2.4 µm. The smallest loss value L = 0.5 dB/cm is obtained for p very close to 0 and D max = 2.4 µm. The possibility of achieving such low loss figures at λ = 3 µm makes the WG suitable for mid-IR applications. Note that the transparency region of the WG can be further extended to longer wavelengths by increasing D max and properly fitting the p parameter to the RI contrasts being used. Figure 7 shows the variation of the confinement loss across the wavelength range 0.3 µm to 3 µm for WG structures with various growth rate parameters. One can see that while for a uniform structure the losses in both O and E polarisations are below 1 dB/cm for wavelengths up to approximately 1.8 µm, the spectral range where the losses are acceptably low is extended up to λ = 3 µm for a structure with p = 0.01. The optimisation of the growth rate parameter results in a reduction of the losses in both polarisations by two orders of magnitude at λ = 3 µm. One may also notice that while the confinement losses of a WG with p = 0.5 are lower than those of a WG with p = 1 for wavelengths below 2.6 µm, a reversal of trend happens at the wavelengths above 2.6 µm. Further, Fig. 6 highlights the distinctly different behaviours of the confinement losses in the O and E polarisations at the low wavelengths. The confinement loss for the E wave stops decreasing with decreasing wavelength below some wavelength which is  The results presented so far were obtained by assuming that the RI change induced in the material by direct fs laser inscription is a real value. In fact, as mentioned in the introduction, the fs irradiation always induces material absorption, which needs to be accounted for in the WG design. To this end, we computed the confinement loss in a WG with p = 0 and p = 0.5 at λ = 1.55 µm, and with p = 0.01 at λ = 3 µm for a range of fs laser-induced loss values. The results are shown in Fig. 8, which reveals that the higher is the confinement loss of the WG, the lower is the WG sensitivity to the imaginary part of the induced RI contrast. Indeed, induced losses of up to 1 dB/cm do not affect the confinement loss at λ = 3 µm, whereas the effect of induced loss is more important at λ = 1.55 µm, where the WG exhibits lower confinement loss. Fs laser−induced loss (dB/cm) Confinement loss (dB/cm) p=0.01, λ=3µm p=0.5 , λ=1.55µm p=0, λ=1.55µm We note that the imaginary part of the modified RI could be measured by using, for example, the Born scattering interferometry method recently proposed in [16].
Finally, we included in our design procedure also the relationship between induced RI contrast δ n and track size D. Indeed, as we mentioned before, the dependence of δ n and D on the laser pulse energy makes such parameters correlated if the sample translation speed is fixed (Eq. (1)). Note that it is possible to experimentally trim these parameters to the desired values by tuning both the laser pulse energy and the sample translation speed, as both produce albeit connected but not identical changes to δ n and D. Simulation results are presented in Fig. 9, which shows the variation of the confinement losses in the O and E polarisations as a function of wavelength for a seven-ring WG structure with a growth rate parameter of p = 0.01 and where the RI contrast was changed from one ring of tracks to another following the change in the track size. It is seen that the spectral region where the losses in both polarisations are below 1 dB/cm extends up to 3.5 µm for this optimised structure.

Conclusions
We have presented a numerical approach to the optimisation of the guiding properties of depressed-index cladding WGs that can be fabricated in a z-cut LiNbO 3 crystal by direct fs laser inscription. The approach accounted for both a suitable variation of the track size among different cladding layers, the relationship between track size and induced RI contrast, and the losses induced on the tracks by fs irradiation. We have shown that the spectral region where the confinement losses in both O and E polarisations are acceptably low (below 1 dB/cm) can extend up to a wavelength of 3.5 µm for optimised, hexagonal WG structures with seven rings of tracks. This makes such structures suitable for mid-IR applications. We note that the possibility of further extending the low-loss operation region of these WGs depends on the ability to experimentally achieve higher RI contrasts than those that are feasible by use of current fs micro-fabrication technology, as well as to increase the track sizes.