# Analysis of Spatially Doped Fused Silica Fiber Optic by Means of a Hamiltonian Formulation of the Helmholtz Equation

Authors: Raúl E. Jiménez-Mejía; Rodrigo Acuna Herrera; Pedro Torres
Journal: Advances in Materials Science and Engineering
 (2018)

---

## Abstract

This paper discusses an alternative method for calculating modal parameters in optical fibers such as propagation constants, transverse distributions, and anisotropy, due to linear and nonlinear phenomena acting as perturbations caused by doped silica regions. This method is based on a Hamiltonian formulation of the Helmholtz equation and the stationary perturbation theory, which allows a full-vectorial description of the electric field components when linear anisotropic inhomogeneities and Kerr nonlinearity are included. Linear and nonlinear parameters can be found for each propagating mode, and its accuracy has been successfully tested when compared to numerical calculations from the vector finite element method, and the results are published in the literature. This method facilitates the calculation of the spatial-distributed perturbation effects on individual electric field components for each propagating mode.


---

## Body

## 1. Introduction

The characteristics of propagation of light in optical waveguides are governed by Maxwell’s equations. In the case of a monochromatic wave with time dependence ofexpjωt, these equations can be written as(1a)∇×∇×E=μ0ω2Dr→,(1b)∇×1ϵr→∇×Hr→=μ0ω2Hr→,(1c)∇⋅Dr→=0,(1d)∇⋅Hr→=0,where Er→ and Hr→ are the electric and magnetic field vectors and Dr→ represents the electric displacement vector given by Dr→=ϵ0Er→+Pr→, being Pr→ the polarization vector that includes linear and nonlinear responses given by: Pr→=PLr→+PNLr→ [1]. In uniform waveguides, with dependence expjωt±jβz attached to appropriated boundary conditions, solutions for radiated and guided modes in (1) can be found from an eigenproblem written in terms of the wavenumber β. These eigenfunctions are characterized by a specific set of parameters to describe the propagation characteristics as an unique entity such as spatial distribution for each field component, an effective refractive index, and the optical power distribution for each of the propagating modes [1–3]. Analytical solutions for the transverse electric field distributions found for some basic dielectric profiles ϵtr→⊥, such as rectangular slabs and cylindrical waveguides, can be used as basis solutions for finding analytical approaches of relatively more complex transverse dielectric distributions by using the perturbation theory [4–6]; this method allows spanning the perturbed scenario through a linear combination of the unperturbed solutions. In order to implement the perturbation theory to optical waveguides in the same fashion as in the Hamiltonian eigenproblem in quantum mechanics [4], a Hamiltonian formulation of Maxwell’s equations can be proposed [5–8]. In this formulation, the propagating parameters of a waveguide with perturbed transverse dielectric profile ϵ˜tr→⊥ are calculated from the unperturbed waveguide ϵtr→⊥, assuming that the perturbed dielectric profile can be considered as a small change from the unperturbed one as ϵ˜tr→⊥=ϵtr→⊥+Δϵtr→⊥. Such perturbation is included into the equations by means of a pair of very elaborated operators constructed to obtain an eigenproblem representation of (1) with a Hermitian Maxwell Hamiltonian [7, 8]. The first operator B^ acts as a longitudinal projector and the second one A^ includes the transverse characteristics of the waveguide [7]. By using this procedure, Maxwell’s equations can be written as a generalized eigenproblem, expressed in Dirac’s notation, can be handled by means of the standard perturbation theory, and can be written as(2a)A^Ψ=βB^Ψ,(2b)Ψ=EtHt,with an orthogonality condition between two eigenkets given by(3)βi−βj∬Eit×Hjt∗+Ejt∗×Hit⋅i^zdA=0,where βi and βj are nondegenerate propagation constants, and Eit, Hit, Ejt, and Hjt are their corresponding transverse mode distributions. The result in (3) coincides with the well-known orthogonality condition between two propagating modes in the coupled mode theory (CMT) [2, 9] and constitutes the basis function form for the perturbation expansions [7].Although the formulation in (2) allows the application of the stationary perturbation theory, one of the difficulties in its application is the assessment of the perturbation effects upon the field distribution individually, because the orthogonality condition cannot be applied in the usual sense between electric or magnetic field independently. This is mainly because the definition of the eigenkets in (2b) leads to the perturbation expansion coefficients to be calculated based in a condition of power-mode independence between two propagating modes, as stated in (3), instead of the effects upon the field individually. Evaluation of the perturbation effects over the electric and magnetic field in a separated way constitutes an advantage in the analysis of dielectric waveguides, particularly for the electric field which can be modified not only in its propagation characteristics but also in its transverse distributions, when immersed in dielectric profile perturbations. In the attempt to deal with (1a) and (1b) as independent eigenproblems, only (1b) satisfies the hermiticity condition for the resultant operator [8]. However, perturbative terms associated to linear and nonlinear contributions of the polarization vector cannot be included into the formulation [8]. On the other hand, (1a) permits the inclusion of the polarization vector effects, but resultant operator does not satisfy hermiticity [7, 8].In this paper, an alternative formulation for calculating the propagating parameters for a perturbed waveguide, as well as the distortions on the transverse electric field distribution is discussed. In the proposed formulation, the effects of perturbative terms upon the modal parameters can be calculated by using the electric field uniquely; thus, expansion coefficients found for the linear combination are relating the expected contribution of the perturbation over each modal electric field distribution. This fact offers a main advantage of the current formulation over conventional methods based on the CMT theory because instead of considering the whole transversal field in expansion coefficients, current formulation simplifies their calculation and facilitates the physical interpretation of the external perturbation over each electric field component that propagates throughout the fiber. As a consequence, the effects of the perturbation terms upon the waveguide modal parameters are more intuitive to be understood and calculated. This formulation can be obtained by operating over (1a) to yield the Helmholtz equation, and after an analogous Hamiltonian formulation of the Helmholtz equation (HFHE), the resultant Hamiltonian operator is reduced to be the Laplacian, and the rest of the terms can be included in the perturbation operator in a very natural manner. By means of this mathematical treatment, field discontinuities and full-vectorial perturbations, concerned to the polarization vector, can be directly included into the perturbative term. Therefore, it is possible to calculate full-vectorial perturbed waveguides scenarios from simple solutions to the Helmholtz equation for each electric field component in the unperturbed situation. A set of numerical experiments for fiber optics waveguides are discussed to show the accuracy of the HFHE when different types of perturbative terms such as inhomogeneities, anisotropies, and nonlinearities are included. Results obtained with HFHE shows an excellent agreement with FEM simulations, and the results are reported in literature.


## 2. Hamiltonian Formulation for the Full-Vectorial Helmholtz Equation

The polarization vectorPr→ is responsible to include linear and nonlinear response of dielectrics. In the case of isotropic materials, the linear term of the polarization vector can be expressed through a scalar relation with the electric field, that is, PLr→=ϵ0χEr→, where χ is the first-order electrical susceptibility of the material. However, when anisotropies are present, the polarization vector and the electric field vector must be related through a tensorial relationship given by PLr→=ϵ0χ^1Er→ [1, 10]. On the other hand, nonlinear response of the material implies to take into account a set of relatively more complex relations with the electric field [1]. In the majority of dielectrics, the nonlinear polarization vector can be written in terms of powers of the electric field as PNLr→=PE2r→,E3r→,…,χ^2,χ^3,…, where χ^j corresponds to the high-order susceptibility tensors. Apart from their natural resonances, this relation can be given by a sum of these terms: PNLr→=ϵ0χ^2Er→Er→+χ^3Er→Er→Er→+⋯ [1]. In the following, we show how to deal with polarization vector effects in the resultant eigenvalue problem when calculating modal parameters in optical waveguides.Consider az-directed propagating wave with propagation constant β and applying the vectorial identity ∇×∇×A=∇∇⋅A−∇2A, the eigenproblem in (1a) can be written as(4a)∇⊥2+k021+χ−β2+W^tEr→⊥=0,(4b)W^t=∇−∇ϵr→⊥ϵr→⊥+ω2Δχ^1c2+ω2c2χ^2:Er→⊥+χ^3Er→⊥:Er→⊥+⋯,where k0=ω/c is the wavenumber, c the speed of light in vacuum; β=neffk0 is the propagation constant through the waveguide with neff as the mode effective refractive index and χ is the linear, homogeneous, and isotropic susceptibility. The additional term W^t includes any inhomogeneity associated to the spatial distribution of the permittivity and polarization vector effects. Linear anisotropies are included through the term Δχ^1, and the nonlinearities are taken into account according to the electric field dependence [1]. This mathematical artifice permits to deal with the Helmholtz equation in the same way to the Hamiltonian eigenvalue problem in quantum mechanics [4], where the unperturbed scenario corresponds to the assumption of isotropic, linear, and lossless material, which reduces (4) to the well-known Helmholtz equation for each electric field component:(5a)∇⊥2+k02n2−β2Er→⊥=0,(5b)∇⊥2Er→⊥=ke2Er→⊥,where n=1+χ is the refractive index of the linear and isotropic medium and ke2=β2−k02n2 acts as the eigenvalue. It is worth noting that the full-vectorial characteristics of the electromagnetic field can still be considered in (5) as long as each electric field component is included in the formulation, any anisotropic effect that involves different electric field components can be addressed by means of the perturbative terms in the polarization vector as presented in (4b). Using Dirac’s notation, (5) can be written as(6)∇⊥2eαj=βj2−k02n2eαj,where eαi are the normalized kets that represent the modal spatial distribution for each electric field component α=x,y,z projected on the coordinate system. These normalized kets can be calculated from the solution to (5b) after imposing the corresponding boundary conditions associated to the waveguide characteristics, which leads to a set of functions Eir→⊥ that are able to propagate throughout the waveguide. For a step-index fiber optic, the set of solutions of each electric field component, Exir→, Eyir→, and Ezir→, describes the spatial distribution associated to the ith propagating mode, which are typically divided in four different families depending on the electric and magnetic field configuration, namely: TE, TM, HE, and EH modes [2, 3]. Once the solutions are found, normalized kets can be calculated by(7)r→eαi=Eαix,y∬Eαix,y2dxdy=Eαix,yAmj,where Amj is commonly called the modal area.As it is stated in the perturbation theory, eigenfunctions are able to form a basis that allows to span the space of solutions, if completeness and orthogonalityeαieαj=δij are satisfied [4], which are the main difficulties in applying the perturbation theory to (6). In the proposed formulation, completeness requirement is relaxed by means of a reduced basis analysis where both conditions can be achieved. It must be highlighted that even in the more elaborated formulation in (2), it is not possible to guarantee the existence of a complete orthogonal basis for the Maxwell’s equations [6, 7]. Therefore, in some scenarios, a reduced set of elements taken from the set of solutions could be enough to describe the perturbed scenario. In practice, few-mode fibers which allow applications with superior features over the standard single-mode fiber (SMF) in optical communications and sensors, are highly suitable for this type of analysis in as much as only some propagating modes can be excited [11].On the other hand, orthogonality condition must be tested between the set of solutions to (6). In the case of EH and HE families, it is easy to demonstrate orthogonality since their angular dependence leads to a null inner product independently of the radial functions [3]. However, when the orthogonality condition is tested between solutions of the same family, it will depend on Bessel’s functions associated to the radial coordinate [3]:(8)elieljr<a=πAliAlj∫0aJνuliarJνuljarrdr,elieljr>a=πAliAlj∫a∞KνwliarKνwljarrdr,where Al,Am are complex constants; a is the fiber core radius; uli, ulj, wli, and wlj are different roots of the transcendental equation for the EH or HE hybrid mode [3]. A simple inspection over the integral in (8) shows that the inner product between two general modes vanishes only when uli and ulj are the roots of the Bessel function Jνr with r∈0,a, which is not the case for the step-index fibers because uli and ulj are determined by the boundary conditions at the core-cladding interface and in general, they do not coincide with a root of the Bessel function [12]. This fact shows that hermiticity of the Laplacian operator in step-index fiber waveguides underlies on the boundary conditions, which supports the discussion presented in [5] about the nonhermiticity of the operator that results from (2a) after multiplying from the left by B^−1. Since the Laplacian operator in step-index fibers is non-Hermitian, it is not possible to use the typical procedures for finding their corresponding correction terms to the eigenvalue and eigenstate [4]. Indeed, it is required more sophisticated expressions in the formulation of perturbative terms in order to deal with non-Hermitian operators and propose an approached solution by using the whole set of nonorthogonal eigenstates [13]. Another possibility is to find a set of orthonormal functions through a systematic orthogonalization process such as Gramm–Schmidt [4], which allows to span each eigenstate eαi as a linear combination of a new orthonormal set of functions vαi that could be used as an orthogonal basis. This procedure works fine when degeneracy is present. However, in the case of nondegeneracy, this new-basis elements will not constitute eigenstates of the zeroth order Hamiltonian [13, 14].A more practical approach is to select some elements from the set of solutions in order to form a convenient orthogonal reduced basis that guarantees the diagonal matrix representation of the Laplacian operator [4, 8]. This assumption can make sense physically in fiber optics since propagating modes constitute a finite basis and under a small external perturbation, these modes will interact among each other instead of changing the mathematical nature of the solutions. The procedure for selecting solutions consists in testing the hermiticity of the operator between the functions of the set of solutions. By solving (6), a finite number of functions eαj with different eigenvalues is found, kαj2=βj2−k02n2. We can write two different eigenvalue equations: ∇⊥2eαj=kαi2eαj, and ∇⊥2eαj=kαj2eαj. Assuming that each ket eαj has associated a corresponding bra eαi, Hermitian condition of the Laplacian operator can be tested through(9)eαj∇⊥2eαi=eαi∇⊥2eαj∗,kαi2eαjeαi=kαj2∗eαjeαi.Since the eigenvalues obtained in (6) are real valued, conjugation makes no effects in the test. When there is not degeneracy between propagating modes, kαi2≠kαj2, (9) is satisfied only by these spatial distributions that are strictly orthogonal each other. Whether the inner product between two nondegenerate solutions is different from zero, that is, eαieαj≠0, hermiticity test of the Laplacian operator upon these solutions will fail, and they are not suitable functions to be used in the reduced basis. An important case to take into account appears when eαieαj=0, (9) is satisfied but eigenvalues could be degenerated, kαi2=kαj2, which implies to use a different formulation of perturbation theory [4]. Although this is not found in step-index fibers because eigenvalues are different each other, this could happen under the weakly-guiding approximation [3], where degeneracy can be overcome grouping degenerated modes into a new set of propagating modes, known as linearly polarized (LP) modes, which can also be analyzed by using the same formulation. Certainly, the accuracy of the result will depend on the number of modes that are included into the basis as well as the perturbation nature. It will be shown that perturbations involving spatial inhomogeneities, as well as some anisotropies and nonlinearities can be analyzed through this simplification leading to accurate results in the propagation parameters, as long as the optical fiber propagates few modes only. However, other type of perturbations that impose strong changes over the number of propagating modes must be addressed by different approaches.Once the basis of eigenkets,eαj is constructed by orthonormal and nondegenerated eigenfunctions, the first-order stationary perturbation theory permits to calculate the correction terms for the propagating parameters through the expressions [4](10a)β˜2=β02+eα0W^teα0,(10b)e˜α=eα0+∑j=1,j≠0∞eαjW^teα0β02−βj2eαj,where β˜ and e˜α are the propagation constant and mode field distribution after perturbation, respectively. From (10a), we can see that(11)β˜=β01+Δβ2β02≈β0+Δβ22β0,where Δβ2=eα0W^teα0. Here, it was assumed that Δβ2≪β02.It is worth noting that perturbation operatorW^t can include any inhomogeneity, linear anisotropy, and nonlinear effects, depending on which terms are included in (4b).


## 3. Numerical Experiments

This section presents a set of numerical experiments to show the accuracy of the HFHE method described above, in the calculation of modal propagation parameters under the presence of external perturbations. Inhomogeneity and anisotropies of refractive index profiles can be designed to enhance the intermodal interaction, which can exploit the induced distortion on the mode field profile in the design, for instance, of spatial-multiplexing processes and all-optical switching. By means of the present formulation, the external perturbations can be conveniently engineered, such as each component of the electric field in the propagating mode is affected selectively. These perturbations can be achieved by introducing selectively external materials to the fused silica matrix by different doping techniques [15]. Optical poling has been also one of the possible techniques to induce anisotropies and spatially modulated inhomogeneities to oxide glasses [16]. The application of the current formulation is very suitable in few-mode fibers where the amount of supported propagating modes constitutes an orthogonal basis for each of the electric field components. However, its application on single-mode fibers with arbitrary polarization is also possible since the electric field components can be conveniently selected to constitute an orthogonal basis. Different types of analysis were performed in order to illustrate how the perturbations can be included into the analysis, and the accuracies are also discussed.


### 3.1. Linear Inhomogeneities

This section considers a step-index few-mode fiber with the following parameters:ncore=1.46, nclad=1.3, and radius ρ=1.5μm; the wavelength was set to be λ=1.6μm in order to allow only one propagating mode for each family. In this case, the set of solutions for the unperturbed Helmholtz equation is composed by orthogonal functions with propagating modes TE01,TM01,HE11,EH11,HE21. Perturbation consisted in the inclusion of spatial inhomogeneities for the linear refractive index in both the core and the cladding region, such as the perturbation strength, Δn=2nr,θΔχr,θ, was imposed from the unperturbed case Δn=0, to refractive index changes about Δn≈10−2. This term can be included in (4b), where only the first two terms are considered. The presence of inhomogeneities in the spatial distribution of the permittivity makes to appear a polarization charge density at the interfaces between the inhomogeneous regions. It is worth noting that the first term of the operator in (4b) represents the polarization charge that can be found by ρp=−ϵ0∇ϵ/ϵ⋅E. In the case under study, the azimuthal change of susceptibility Δχθ will impose a surface charge density at each azimuthal interface θ=θ0′. Using cylindrical coordinates, polarization charge density at each azimuthal interface can be written by(12)ρpr′,θ′,z′=ϵ0r′ϵ1−ϵ2ϵ2E1θjr′,θ0′,z′δθ−θ0′,where subscript 1,2 defines both regions at the interface in which the azimuthal component of the electric field is directed from 1 to 2. Subscript j relates the mode function for the corresponding propagating mode. As it is stated in (12), the magnitude of the perturbation depends on the relative change of permittivity due to the inhomogeneity. For the refractive index contrasts under consideration, induced polarization charge at the interfaces can be neglected; thus, the perturbation operator can be given by a simple expression W^t=μ0ϵ0ω2Δχ1x,y, which is a spatial-dependent perturbation.In order to compare the accuracy of the proposed approach in the prediction of the effects due to the inhomogeneitiesΔχx,y, the results obtained with the HFHE formulation were compared with those obtained from the vector FEM approach. Comparisons were performed for both the effective refractive index and the distortion in the electric field distribution that undergoes each mode. It is worth noting that the simulations with the FEM are performed assuming that fiber does not suffer any perturbation, that is, inhomogeneities are included as initial conditions of the problem, so it is not rigorously speaking an induced distortion by a perturbation, but an initial distribution of refractive index. This is an important advantage of the HFHE method because it can describe the transition from an initial spatial distribution of the guided mode into a distorted one due to the presence of an external perturbation, which can be used as an strategy for modal division multiplexing [17].Figure1 shows the spatial dependence of the perturbative terms. Once the perturbation is included, both the effective refractive index and the change in the field distribution for each propagating mode are calculated as a function of the perturbation magnitude. Figure 2 presents the dependence of the effective refractive index for each propagating mode as a function of the perturbation strength. From this figure, we can see that our results are in good agreement with the vector FEM solutions. We found a maximum absolute error about 0.1% between them at the strongest perturbation. Sources of the mismatch can come from neglecting the induced polarization charges due to the imposed inhomogeneities in the polarization vector. An additional comparison was performed on the calculation of the Ex field profiles of the guided modes due to the presence of the inhomogeneity. In [7] correction terms were applied only to propagation constants. However, by using HFHE formulation and through the perturbation method described in Section 2, the spatial-distribution correction can be also well estimated after perturbation as presented in Figure 3.Figure 1
Spatial-dependent linear inhomogeneity and orthogonal propagating modes: (a) Unperturbed fiber; (b) perturbed fiber with azimuthal changesΔχθ of periodicity π/2.
(a)
(b)Figure 2
Effective refractive index of the perturbed optical fiber in Figure1(b) for the four propagation modes versus Δn obtained by our HFHE formulation (lines) and calculated by the vector FEM (dots).Figure 3
Mode field distributions (Ex component) of the perturbed optical fiber in Figure 1(b). (a) HE11, (b) TE01, (c) TM01. When these field distributions are compared with the unperturbed ones in Figure 1(c), the effect of the perturbation is to rotate the mode distribution in a clockwise direction.
(a)
(b)
(c)It must be noted that transverse modes are most affected because they do not present symmetry with respect to the spatial inhomogeneity, increasing the magnitude of the coupling coefficient between modes in (10b). When these field distributions are compared with the unperturbed ones, the effect of the perturbation is to rotate the mode in a clockwise direction in which the rotation angle depends on the perturbation magnitude, and it is directly related with the coupling coefficient between TE and TM modes. In practice, this effect can be exploited by doping the fiber in specific regions where inhomogeneities can be controlled externally allowing mode conversion processes.


### 3.2. Linear Anisotropies

Linear anisotropies can be also included in the HFHE formulation, and these can be considered in the second terms of (4b). In this case, the perturbation term is an imposed anisotropy that additionally has a spatial variation in the refractive index. The same spatial dependence shown in Figure 1(b) is used for the anisotropic perturbation, but in this case, the perturbation term is defined as a diagonal tensor given by(13)W^t=μ0ϵ0ω2Δχ1x,y000−Δχ1x,y000Δχ1x,y.This situation differs greatly from the previous case because we must now include the effect of the perturbation on each component of the electric field in the propagating mode exploiting the full-vector characteristics of the formulation. A negative perturbation is included for they-component in order to induce a linear uniaxial birefringence in the waveguide. Figure 4 presents the comparison between the full-vectorial calculation by the FEM and the HFHE for those propagating modes that are most affected.Figure 4
Effective refractive index of the optical fiber with uniaxial birefringence having the same spatial dependence of Figure1(b) versus Δn obtained by the HFHE formulation (lines) and calculated by the vector FEM (dots). Solid lines correspond to Ex component and dashed lines to Ey component.
(a)
(b)As can be seen from Figure4, the results obtained from the two schemes are in good agreement. In this instance, a maximum mismatch about 0.6% was found between them. This case is of great importance in propagation analysis because it allows the calculation of individual effects over the electric field components of the propagating modes under an external perturbation. A significant difference between both predictions occurs when calculating effective refractive index for the EH11 mode. As can be observed in Figure 4, the FEM method does not report birefringence for this mode, and the predicted effective refractive index remains on the calculated values for the x-component in the HFHE method. This discrepancy could be related with the corresponding mode-profile for the hybrid EH11 mode and the spatial distribution of the imposed inhomogeneity. FEM solves the vectorial eigenvalue problem independently of the coordinate axis orientation and unperturbed mode distributions, whereas the HFHE method starts from an analytical solution expressed in well-defined coordinate axes and calculates the effects of the perturbation upon the unperturbed mode distributions leading to correction terms for each electric field component as shown in Figure 4.


### 3.3. Kerr Nonlinearity

In silica, optical fibers’ second-order susceptibility termχ^2 is null due to the material symmetry [1]. However, χ^3 can cause several nonlinear effects, such as four-waves mixing phenomena (FWM), third harmonic generation (THG), self-phase modulation (SPM), and cross-phase modulation (XPM). In practice, FWM and THG effects require phase-matching conditions that are typically of great difficulty to achieve in optical fibers [18]; therefore, Kerr effect is typically included through the refractive index dependence on the electric field intensity, which is related with the SPM. This refractive index change is given by [18](14)Δnω,E=38nℜχxxxx3E2=n¯2E2,where ℜχxxxx3 is the real part of the third-order susceptibility under the assumption of a constant state of polarization and n¯2 is known as the nonlinear index coefficient. In optical fiber analyses, it is very common to report the nonlinear Kerr parameter or simply the nonlinear refractive index, n2, which is calculated by n2=2n¯2/ϵ0nc [18]. To calculate the induced change in susceptibility due to the Kerr-type nonlinearity, it is known that Δχ=2nΔn, thus(15)Δχ=34ℜχxxxx3E2.It has been shown in literature that scalar approach is not always accurate enough to correctly determine the effective refractive index associated with the nonlinear effect [19]. Therefore, a vectorial description of the Kerr nonlinearity must be considered for an accurate description of its associated nonlinear effects. By means of a similar procedure, used in Section 3.2, a full-vectorial relation can be included in the calculation of the effective nonlinear refractive index through a perturbation term given by(16)W^t=μ0ϵ0ω2Δχx,y,where Δχx,y is a diagonal tensor that will affect each of the electric field components based on the material symmetry properties. For silica fibers, the full-vectorial effect of the polarization vector can be arranged through a second-rank diagonal tensor by the following equations [20, 21]:(17)Δχx,yx,y=34ℜχxxxx3Ex2+Ey2+13Ez2,Δχx,yz=34ℜχxxxx313Ex2+Ey2+Ez2.These full-vectorial perturbation terms can be directly included into the third term in the formulation of (4b), and using similar procedure as in Section 3.2. 2, both effective propagation constant and mode distortion correction can be calculated through (10b) when SPM nonlinearity is induced due to the electric field intensity.In order to test the validity of the proposed approach, the effective refractive index by using the HFHE method is calculated for a nonlinear waveguide with the same parameters reported in [20], where an iterative method was used for finding the effective propagation constant based on the input power. Optical fiber parameters are defined as follows: radius a=0.5μm, wavelength λ=1.55μm, linear refractive indexes nco=1.45, ncl=nair=1.00, and the nonlinear Kerr coefficient is given by n2=3.2×10−20m2/W.Figure5 shows the calculation of the effective refractive index for a single-mode fiber as a function of the optical power. The results are compared to those obtained with the iterative solution reported in [20]. Effective refractive index is presented for the Ex field component of the fundamental mode HE11. As can be seen from Figure 5, a very good agreement is found for a wide input power range. In addition, the maximum absolute error is about 0.5% at the highest optical excitation power.Figure 5
Effective refractive indexneff for a silica rod as a function of the input optical power. The dashed line is the neff given by the HFHE formulation, whereas the dots represent the results calculated by a full-vectorial iterative approach (FVI) [20].This result shows the validity of the proposed method for dealing with vectorial perturbation terms in single-mode fibers. It is worth noting that for single-mode fibers, hermiticity analysis of (9) is not required because there is only one propagating mode. Birefringence analysis are typically performed under the assumption of two orthogonal polarized-modes, but in that particular case, the analysis is more suitable to be performed through the formulation in [6].


### 3.4. Nonlinear Parameter Calculation

Nonlinear parameterγ is defined in the analysis of optical pulse propagation when solving the nonlinear propagation equation [10]; particularly for fiber optics, γ is related mainly with the SPM effect. It is worth noting that [20, 22, 23] have proposed different expressions for extending the calculation of γ in order to take into account all full-vectorial Kerr-type perturbations. Indeed, different definitions of γ are still under discussion and research [24]. Typical definition of the nonlinear parameter in single-mode fiber relates the propagation constant change due to the nonlinear effects and the optical power that carries the transverse electromagnetic field [18]. However, when multimode fibers are under analysis mode, interactions due to nonlinearities must be considered [22]. These propositions have been derived based on the power-orthogonality relation discussed in (3), and though they can include full-vectorial properties of the waveguide such as anisotropies and losses, these definitions remain on the perturbation expansion associated to power independence, and as it was discussed in Section 3.1, external perturbations in multimode fibers can cause also spatial-mode distortions when perturbative terms enable the interaction between propagating modes as stated in (10). A simple extension of the definition for γ in [18] can be achieved based on the HFHE method, in which the distortions on spatial distribution can be included in the calculation of the optical power carried by the perturbed mode and consequently in the nonlinear parameter evaluation. With the aim to take into account these possible mode distortions for each j-th propagating mode, the nonlinear parameter γj can be written by(18a)γj=2Δβ∬ℜS˜⋅i^zdA,(18b)γj=eαjW^teαjcϵ0βjAmje˜αjnx,ye˜αj,where S˜=E˜×H˜∗ is the complex Poynting’s vector for the perturbed propagating mode, with E˜ and H˜ the perturbed fields calculated from (10). Magnetic field components can be determined directly from the perturbed electric field by using the Maxwell’s equations (3). For single-mode fibers, and considering that the usual nonlinearities are only in the core region, Poynting’s vector S^ can be simplified under the assumption of dominant electric field components as: S^z=E˜x2/η, being η the intrinsic impedance η=1/ncϵ0. This simplification allows to write the optical power for the resultant field in terms of perturbed kets, leading to the expression in (18b). This latter is similar to that proposed in [22] and presents the same advantages in the analysis of the nonlinearity separated into parts, namely, contributions of linear and nonlinear regions. Besides, (18b) allows associating the intermode interaction effects through the inclusion of perturbed kets in the denominator. From (18b), it is also possible to reproduce the expression proposed in [18] for the nonlinear parameter. This can be achieved under the following considerations: (1) No modal interactions are considered and only the x-component of the electric field is propagating, i.e. e˜xj≈e˜xj; (2) Kerr effect is induced due to x-polarized field, thus the perturbation operator can be written as W^t=ω2/c22nx,yn¯2x,yExx,y2; (3) nonlinear effects occur only in the core region; (4) the propagation constant for single-mode fiber can be approached by β≈nω/c.Under these assumptions (18) for the nonlinear parameter calculation leads to:(19)γ=ωc2n¯2ϵ0n∬Exx,y4dxdy∬Exx,y2dxdy2,γ=ωn2cAeff,Aeff=∬Exx,y2dxdy2∬Exx,y4dxdy,where Aeff is the standard definition of the effective area [18].In order to show the accuracy of the derived expression in (18) for calculating nonlinear parameter when full-vectorial Kerr nonlinearity is considered, a set of step-index rod configurations with a high-index contrast were considered. Rod configurations were the same as those analyzed in [22] by using the expression for γ in [23]. Figure 6 displays the comparison between γ obtained with (18) and from those reported in [22]. As it can be seen in Figure 6, the calculation through (18) allows to reproduce the nonlinear parameter with relatively good accuracy. Maximum differences were found when higher refractive contrasts were considered, which could be due to the overestimation of the optical power transported in the core region of the rods with high refractive contrasts, which affects the calculation of the intensity-dependent perturbation. Experimental tests reported in [25] show that effectively high nonlinearities are expected when full-vectorial Kerr effect description is included.Figure 6
Nonlinear parameterγ at λ=800nm as a function of the core diameter obtained by the HFHE formulation (continuous line) and results given in [22, 23] (dots). (a) Silica rod: nc=1.45, ncl=1.0, n2=2.6×10−20m2/W. (b) Bismuth rod: nc=2.05, ncl=1.0, n2=3.2×10−19m2/W. (c) Silicon rod: nc=3.45, ncl=1.45, n2=4.5×10−18m2/W.
(a)
(b)
(c)


## 3.1. Linear Inhomogeneities

This section considers a step-index few-mode fiber with the following parameters:ncore=1.46, nclad=1.3, and radius ρ=1.5μm; the wavelength was set to be λ=1.6μm in order to allow only one propagating mode for each family. In this case, the set of solutions for the unperturbed Helmholtz equation is composed by orthogonal functions with propagating modes TE01,TM01,HE11,EH11,HE21. Perturbation consisted in the inclusion of spatial inhomogeneities for the linear refractive index in both the core and the cladding region, such as the perturbation strength, Δn=2nr,θΔχr,θ, was imposed from the unperturbed case Δn=0, to refractive index changes about Δn≈10−2. This term can be included in (4b), where only the first two terms are considered. The presence of inhomogeneities in the spatial distribution of the permittivity makes to appear a polarization charge density at the interfaces between the inhomogeneous regions. It is worth noting that the first term of the operator in (4b) represents the polarization charge that can be found by ρp=−ϵ0∇ϵ/ϵ⋅E. In the case under study, the azimuthal change of susceptibility Δχθ will impose a surface charge density at each azimuthal interface θ=θ0′. Using cylindrical coordinates, polarization charge density at each azimuthal interface can be written by(12)ρpr′,θ′,z′=ϵ0r′ϵ1−ϵ2ϵ2E1θjr′,θ0′,z′δθ−θ0′,where subscript 1,2 defines both regions at the interface in which the azimuthal component of the electric field is directed from 1 to 2. Subscript j relates the mode function for the corresponding propagating mode. As it is stated in (12), the magnitude of the perturbation depends on the relative change of permittivity due to the inhomogeneity. For the refractive index contrasts under consideration, induced polarization charge at the interfaces can be neglected; thus, the perturbation operator can be given by a simple expression W^t=μ0ϵ0ω2Δχ1x,y, which is a spatial-dependent perturbation.In order to compare the accuracy of the proposed approach in the prediction of the effects due to the inhomogeneitiesΔχx,y, the results obtained with the HFHE formulation were compared with those obtained from the vector FEM approach. Comparisons were performed for both the effective refractive index and the distortion in the electric field distribution that undergoes each mode. It is worth noting that the simulations with the FEM are performed assuming that fiber does not suffer any perturbation, that is, inhomogeneities are included as initial conditions of the problem, so it is not rigorously speaking an induced distortion by a perturbation, but an initial distribution of refractive index. This is an important advantage of the HFHE method because it can describe the transition from an initial spatial distribution of the guided mode into a distorted one due to the presence of an external perturbation, which can be used as an strategy for modal division multiplexing [17].Figure1 shows the spatial dependence of the perturbative terms. Once the perturbation is included, both the effective refractive index and the change in the field distribution for each propagating mode are calculated as a function of the perturbation magnitude. Figure 2 presents the dependence of the effective refractive index for each propagating mode as a function of the perturbation strength. From this figure, we can see that our results are in good agreement with the vector FEM solutions. We found a maximum absolute error about 0.1% between them at the strongest perturbation. Sources of the mismatch can come from neglecting the induced polarization charges due to the imposed inhomogeneities in the polarization vector. An additional comparison was performed on the calculation of the Ex field profiles of the guided modes due to the presence of the inhomogeneity. In [7] correction terms were applied only to propagation constants. However, by using HFHE formulation and through the perturbation method described in Section 2, the spatial-distribution correction can be also well estimated after perturbation as presented in Figure 3.Figure 1
Spatial-dependent linear inhomogeneity and orthogonal propagating modes: (a) Unperturbed fiber; (b) perturbed fiber with azimuthal changesΔχθ of periodicity π/2.
(a)
(b)Figure 2
Effective refractive index of the perturbed optical fiber in Figure1(b) for the four propagation modes versus Δn obtained by our HFHE formulation (lines) and calculated by the vector FEM (dots).Figure 3
Mode field distributions (Ex component) of the perturbed optical fiber in Figure 1(b). (a) HE11, (b) TE01, (c) TM01. When these field distributions are compared with the unperturbed ones in Figure 1(c), the effect of the perturbation is to rotate the mode distribution in a clockwise direction.
(a)
(b)
(c)It must be noted that transverse modes are most affected because they do not present symmetry with respect to the spatial inhomogeneity, increasing the magnitude of the coupling coefficient between modes in (10b). When these field distributions are compared with the unperturbed ones, the effect of the perturbation is to rotate the mode in a clockwise direction in which the rotation angle depends on the perturbation magnitude, and it is directly related with the coupling coefficient between TE and TM modes. In practice, this effect can be exploited by doping the fiber in specific regions where inhomogeneities can be controlled externally allowing mode conversion processes.


## 3.2. Linear Anisotropies

Linear anisotropies can be also included in the HFHE formulation, and these can be considered in the second terms of (4b). In this case, the perturbation term is an imposed anisotropy that additionally has a spatial variation in the refractive index. The same spatial dependence shown in Figure 1(b) is used for the anisotropic perturbation, but in this case, the perturbation term is defined as a diagonal tensor given by(13)W^t=μ0ϵ0ω2Δχ1x,y000−Δχ1x,y000Δχ1x,y.This situation differs greatly from the previous case because we must now include the effect of the perturbation on each component of the electric field in the propagating mode exploiting the full-vector characteristics of the formulation. A negative perturbation is included for they-component in order to induce a linear uniaxial birefringence in the waveguide. Figure 4 presents the comparison between the full-vectorial calculation by the FEM and the HFHE for those propagating modes that are most affected.Figure 4
Effective refractive index of the optical fiber with uniaxial birefringence having the same spatial dependence of Figure1(b) versus Δn obtained by the HFHE formulation (lines) and calculated by the vector FEM (dots). Solid lines correspond to Ex component and dashed lines to Ey component.
(a)
(b)As can be seen from Figure4, the results obtained from the two schemes are in good agreement. In this instance, a maximum mismatch about 0.6% was found between them. This case is of great importance in propagation analysis because it allows the calculation of individual effects over the electric field components of the propagating modes under an external perturbation. A significant difference between both predictions occurs when calculating effective refractive index for the EH11 mode. As can be observed in Figure 4, the FEM method does not report birefringence for this mode, and the predicted effective refractive index remains on the calculated values for the x-component in the HFHE method. This discrepancy could be related with the corresponding mode-profile for the hybrid EH11 mode and the spatial distribution of the imposed inhomogeneity. FEM solves the vectorial eigenvalue problem independently of the coordinate axis orientation and unperturbed mode distributions, whereas the HFHE method starts from an analytical solution expressed in well-defined coordinate axes and calculates the effects of the perturbation upon the unperturbed mode distributions leading to correction terms for each electric field component as shown in Figure 4.


## 3.3. Kerr Nonlinearity

In silica, optical fibers’ second-order susceptibility termχ^2 is null due to the material symmetry [1]. However, χ^3 can cause several nonlinear effects, such as four-waves mixing phenomena (FWM), third harmonic generation (THG), self-phase modulation (SPM), and cross-phase modulation (XPM). In practice, FWM and THG effects require phase-matching conditions that are typically of great difficulty to achieve in optical fibers [18]; therefore, Kerr effect is typically included through the refractive index dependence on the electric field intensity, which is related with the SPM. This refractive index change is given by [18](14)Δnω,E=38nℜχxxxx3E2=n¯2E2,where ℜχxxxx3 is the real part of the third-order susceptibility under the assumption of a constant state of polarization and n¯2 is known as the nonlinear index coefficient. In optical fiber analyses, it is very common to report the nonlinear Kerr parameter or simply the nonlinear refractive index, n2, which is calculated by n2=2n¯2/ϵ0nc [18]. To calculate the induced change in susceptibility due to the Kerr-type nonlinearity, it is known that Δχ=2nΔn, thus(15)Δχ=34ℜχxxxx3E2.It has been shown in literature that scalar approach is not always accurate enough to correctly determine the effective refractive index associated with the nonlinear effect [19]. Therefore, a vectorial description of the Kerr nonlinearity must be considered for an accurate description of its associated nonlinear effects. By means of a similar procedure, used in Section 3.2, a full-vectorial relation can be included in the calculation of the effective nonlinear refractive index through a perturbation term given by(16)W^t=μ0ϵ0ω2Δχx,y,where Δχx,y is a diagonal tensor that will affect each of the electric field components based on the material symmetry properties. For silica fibers, the full-vectorial effect of the polarization vector can be arranged through a second-rank diagonal tensor by the following equations [20, 21]:(17)Δχx,yx,y=34ℜχxxxx3Ex2+Ey2+13Ez2,Δχx,yz=34ℜχxxxx313Ex2+Ey2+Ez2.These full-vectorial perturbation terms can be directly included into the third term in the formulation of (4b), and using similar procedure as in Section 3.2. 2, both effective propagation constant and mode distortion correction can be calculated through (10b) when SPM nonlinearity is induced due to the electric field intensity.In order to test the validity of the proposed approach, the effective refractive index by using the HFHE method is calculated for a nonlinear waveguide with the same parameters reported in [20], where an iterative method was used for finding the effective propagation constant based on the input power. Optical fiber parameters are defined as follows: radius a=0.5μm, wavelength λ=1.55μm, linear refractive indexes nco=1.45, ncl=nair=1.00, and the nonlinear Kerr coefficient is given by n2=3.2×10−20m2/W.Figure5 shows the calculation of the effective refractive index for a single-mode fiber as a function of the optical power. The results are compared to those obtained with the iterative solution reported in [20]. Effective refractive index is presented for the Ex field component of the fundamental mode HE11. As can be seen from Figure 5, a very good agreement is found for a wide input power range. In addition, the maximum absolute error is about 0.5% at the highest optical excitation power.Figure 5
Effective refractive indexneff for a silica rod as a function of the input optical power. The dashed line is the neff given by the HFHE formulation, whereas the dots represent the results calculated by a full-vectorial iterative approach (FVI) [20].This result shows the validity of the proposed method for dealing with vectorial perturbation terms in single-mode fibers. It is worth noting that for single-mode fibers, hermiticity analysis of (9) is not required because there is only one propagating mode. Birefringence analysis are typically performed under the assumption of two orthogonal polarized-modes, but in that particular case, the analysis is more suitable to be performed through the formulation in [6].


## 3.4. Nonlinear Parameter Calculation

Nonlinear parameterγ is defined in the analysis of optical pulse propagation when solving the nonlinear propagation equation [10]; particularly for fiber optics, γ is related mainly with the SPM effect. It is worth noting that [20, 22, 23] have proposed different expressions for extending the calculation of γ in order to take into account all full-vectorial Kerr-type perturbations. Indeed, different definitions of γ are still under discussion and research [24]. Typical definition of the nonlinear parameter in single-mode fiber relates the propagation constant change due to the nonlinear effects and the optical power that carries the transverse electromagnetic field [18]. However, when multimode fibers are under analysis mode, interactions due to nonlinearities must be considered [22]. These propositions have been derived based on the power-orthogonality relation discussed in (3), and though they can include full-vectorial properties of the waveguide such as anisotropies and losses, these definitions remain on the perturbation expansion associated to power independence, and as it was discussed in Section 3.1, external perturbations in multimode fibers can cause also spatial-mode distortions when perturbative terms enable the interaction between propagating modes as stated in (10). A simple extension of the definition for γ in [18] can be achieved based on the HFHE method, in which the distortions on spatial distribution can be included in the calculation of the optical power carried by the perturbed mode and consequently in the nonlinear parameter evaluation. With the aim to take into account these possible mode distortions for each j-th propagating mode, the nonlinear parameter γj can be written by(18a)γj=2Δβ∬ℜS˜⋅i^zdA,(18b)γj=eαjW^teαjcϵ0βjAmje˜αjnx,ye˜αj,where S˜=E˜×H˜∗ is the complex Poynting’s vector for the perturbed propagating mode, with E˜ and H˜ the perturbed fields calculated from (10). Magnetic field components can be determined directly from the perturbed electric field by using the Maxwell’s equations (3). For single-mode fibers, and considering that the usual nonlinearities are only in the core region, Poynting’s vector S^ can be simplified under the assumption of dominant electric field components as: S^z=E˜x2/η, being η the intrinsic impedance η=1/ncϵ0. This simplification allows to write the optical power for the resultant field in terms of perturbed kets, leading to the expression in (18b). This latter is similar to that proposed in [22] and presents the same advantages in the analysis of the nonlinearity separated into parts, namely, contributions of linear and nonlinear regions. Besides, (18b) allows associating the intermode interaction effects through the inclusion of perturbed kets in the denominator. From (18b), it is also possible to reproduce the expression proposed in [18] for the nonlinear parameter. This can be achieved under the following considerations: (1) No modal interactions are considered and only the x-component of the electric field is propagating, i.e. e˜xj≈e˜xj; (2) Kerr effect is induced due to x-polarized field, thus the perturbation operator can be written as W^t=ω2/c22nx,yn¯2x,yExx,y2; (3) nonlinear effects occur only in the core region; (4) the propagation constant for single-mode fiber can be approached by β≈nω/c.Under these assumptions (18) for the nonlinear parameter calculation leads to:(19)γ=ωc2n¯2ϵ0n∬Exx,y4dxdy∬Exx,y2dxdy2,γ=ωn2cAeff,Aeff=∬Exx,y2dxdy2∬Exx,y4dxdy,where Aeff is the standard definition of the effective area [18].In order to show the accuracy of the derived expression in (18) for calculating nonlinear parameter when full-vectorial Kerr nonlinearity is considered, a set of step-index rod configurations with a high-index contrast were considered. Rod configurations were the same as those analyzed in [22] by using the expression for γ in [23]. Figure 6 displays the comparison between γ obtained with (18) and from those reported in [22]. As it can be seen in Figure 6, the calculation through (18) allows to reproduce the nonlinear parameter with relatively good accuracy. Maximum differences were found when higher refractive contrasts were considered, which could be due to the overestimation of the optical power transported in the core region of the rods with high refractive contrasts, which affects the calculation of the intensity-dependent perturbation. Experimental tests reported in [25] show that effectively high nonlinearities are expected when full-vectorial Kerr effect description is included.Figure 6
Nonlinear parameterγ at λ=800nm as a function of the core diameter obtained by the HFHE formulation (continuous line) and results given in [22, 23] (dots). (a) Silica rod: nc=1.45, ncl=1.0, n2=2.6×10−20m2/W. (b) Bismuth rod: nc=2.05, ncl=1.0, n2=3.2×10−19m2/W. (c) Silicon rod: nc=3.45, ncl=1.45, n2=4.5×10−18m2/W.
(a)
(b)
(c)


## 4. Conclusions

This paper presented a discussion for including anisotropic inhomogeneities and nonlinearities in fiber-optic propagating modes by using the stationary perturbation theory for each component of the electric field. This method takes advantage of the Hamiltonian formulation of the Helmholtz equation (HFHE) to include the perturbation of the polarization vector for each electric field component independently. As it was discussed above, this formulation requires the assumption of an orthogonal basis for each component of the propagating modes. Therefore, its validity is closely related to the number of modes supported by the optical fiber under analysis. The results obtained with the HFHE were successfully contrasted with numerical results from the vector finite element method and the previous literature. As it was shown, for a few-mode fiber with an inhomogeneous perturbation with refractive index changes aroundΔn=0.01, only four modes were required to constitute the orthogonal basis to implement the formulation. The approach predicted the perturbed propagation constant and their spatial distribution with a maximum error about 0.1%, and for the same fiber with induced anisotropy, the error was under 0.6%. For a single-mode fiber with a full-vectorial fundamental mode, it was shown that the presented approach reproduces propagation constant changes due to strong Kerr-type nonlinearities with a maximum error of 0.5%. It was also shown that the nonlinear parameter can be written in terms of the perturbed fields, and through current formulation, it leads to the standard expression for single-mode fibers in the scalar limit.



---

*Source: 5806947-2018-05-14.xml*