Numerical Analysis of Mode Coupling in Multimode Graded Index Optical Fibers with Bending

In this paper, we present a new and more realistic theoretical framework for lightwave propagation in a multimode graded index (GRIN) optical fiber when the fundamental mode is selectively excited into the fiber with constant radius bending, causing coupling between the various modes of the fiber. First, a wave equation is formulated to represent the light behavior in the GRIN fiber and solved numerically by an eigenvalue method using a difference equation representation of the differential equation, resulting in the mode amplitudes. Next, the local normal mode fields at a succession of infinitesimal corner bends are matched to calculate the bending-induced mode coupling. Finally, the power in the fundamental mode vs the distance along the propagation direction is calculated, assuming that this is the only mode that is excited initially. For large bend diameters, the output data indicates that, when the fundamental mode is excited, the light remains in a set of low-order modes. However, for small radius bends (less than 1 cm), the oscillations become irregular and the power is not completely located in the fundamental mode when z > 0. While this is consistent with experimental observations, it contradicts the predictions of the previous oversimplified model.


Introduction
The research in this paper was motivated by an interest in using multimode graded index optical fibers as the transmission medium in intrusion-alarmed communication systems. (1,2) In such systems, data is transmitted using the lowest order mode, while the intrusion-monitor signal is simultaneously transmitted using high-order modes. An attempted intrusion results in greater attenuation of the monitor signal, thereby setting off an alarm. However, for such a system to be effective, it is essential that a perturbation (e.g., fiber bend) introduced by an intruder removes the light from the high-order modes more effectively than that from the low-order modes. Thus, in this context, it is meaningful to investigate the bend-induced mode couplings and lightwave behavior in graded index multimode fibers undergoing bending.
Optical propagation in graded index fibers subject to mode-coupling perturbations or mode-selective attenuation has already been studied for decades owing to its implications for communication systems using these fibers. Furthermore, not only means of communication links but also bending properties have been used for optical fiber sensor applications such as pressure, strain, vibration, and displacement sensors.
As a result, several studies of graded index fiber modes and bending have been reported with various objectives. (3)(4)(5)(6)(7)(8)(9)(10)(11) Mode couplings in GRIN fibers have been studied. (5,6) The effects of bending and its sensor application for graded index fibers in which a few low-order modes are excited have been discussed. (8) Studies of the bending effects in optical fibers have also been reported. (9)(10)(11) In addition, mode coupling at bends in optical fibers supporting one or only a few guided modes has been analyzed by considering the local normal modes for the corresponding straight waveguide. (11) From a previous theoretical study, (2) it is known that for sufficiently large bend radii, little coupling from the fundamental mode to high-order modes is expected. For small bend radii, the fundamental mode power has been reported to oscillate as a function of the propagation distance in the bent fiber, where the period of oscillation is equal to the beat length of the graded index fiber modes. However, this model is oversimplified on the basis of the assumption that the quadratic variation in refractive index extends to infinity. Accordingly, the objective of this research is to develop a more realistic model of mode coupling in a graded index optical fiber experiencing bending around a constantdiameter mandrel. (2) In particular, the calculations are related to bend-induced mode coupling in a graded index fiber with fundamental mode excitation.

Formulation
The objective is to determine the effect of the "finite core" on bend-induced mode coupling in a graded index fiber. This deals with an unresolved issue from the previous model, (2) which treats the bend-induced mode coupling problem for an index variation, which extends to infinity ("infinite core case"), where n 1 is the refractive index of the core and a is the core radius. It is already known that for the "finite core" case with strong intermode coupling (small radius bend), losses occur in the low-order modes, and these losses need to be quantified. To solve the "finite core" issue, solutions to the wave equation are needed for the mode fields of a graded index fiber with an index variation, where a is the core radius and b is the fiber radius. The mode fields u lm (r,θ) are the solution to where k 0 = 2π/λ, and r and θ are the cylindrical coordinates. The solutions to the equation can be written as where m is an arbitrary integer and l is a positive integer, such that the number of "zeroes" in f lm (r) equals l for any m. The equation for f then becomes This f is only a function of r; thus, it can be determined by the given boundary conditions. Inside the core, the solutions have the form f (r) = P(r)e −α 2 r 2 /2 , where P(r) is a polynomial in r and α is defined by α = (2πNA/λa) 1/2 , where NA is the numerical aperture of the fiber and is given by NA = (2n 1 Δn) 1/2 . For example, if NA equals 0.2, λ = 1.3 μm, and the core radius equals 25 μm, then α = 0.196 /μm. Outside the core, the solutions can be expressed in terms of modified Bessel functions. The boundary conditions are that f and f′ are continuous at r = a, and f(b) = 0. As Bessel functions are difficult to deal with, this equation is solved as an eigenvalue problem by numerical analysis.

Computational analysis
As the mode amplitude is obtained by solving differential eq. (6) induced from the wave equation, the best way to solve the equation numerically is treating it as an eigenvalue problem using a difference equation representation of the differential equation. (12) This avoids having to deal with the complexity of Bessel functions. Thus, to convert eq. (6) into a difference equation form, the continuous terms need to be replaced with discrete terms. Hence, the first and second derivative terms are converted into discrete terms as ( , and r j = jΔr. By collecting the terms and replacing the expressions in eqs. (7) and (8) with related terms, eq. (6) can be rewritten as where Q = 2 − n j ·(Δr) 2 , Q p = 1 + 0.5·Δr/r j , and Q m = 0.5·Δr/r j − 1. From eq. (9), f j+1 can be calculated by knowing f j−1 and f j . Therefore, f 0 and f 1 are determined to get the process started in the m = 0 and m ≠ 0 cases, respectively. When m = 0, assuming f (r) = e −α 2 r 2 as a trial solution for a small r, then f 0 = 1 and . Likewise, f 0 = 0 and f 1 = (∆r) m 2 when m ≠ 0. The following procedure is then executed to solve the difference equation.
Here, when performing the above procedure, Δr = 0.1 micron and ε = 10 −6 or less. The program is iterated until the new trial Δ(β 2 ) value makes the f N value converge near 0 for each m value and all possible numbers of zero crossings.

Coupled Mode Analysis
Next, the coupling between modes caused by a constant-diameter bend is considered. The normalized field distributions u(r,θ) for various modes can be stored and integrated to determine the coupling constants k ij . Once these are known, the equation representing the amplitude variation of each mode as a function of the fiber distance, such as eq. (18) from ref. 2, can be solved numerically to determine the power in each mode as a function of z. Of particular interest is when A 00 = 1; A ij (0) = 0 for (i,j) ≠ (0,0), in which case, all the power is in the lowest order mode at z = 0. To determine the power distribution between the fiber modes as a function of distance in a uniform bend, the coupled mode analysis of a circular bend in a fiber is performed.
First, a single corner bend is used, where the fiber axis changes direction through the angle φ. Figure 1 shows the local coordinate axis for the corner bend. For the light approaching the bend, it is assumed that the electric mode field distribution is u lm (r,θ) = f lm (r)e imθ , whereas after the bend, the electric mode field distribution is where A l'm' s are constants.
A bend in an optical fiber can be considered as a continuous succession of infinitesimal corner bends. Thus, at a corner bend under a field matching condition for a small φ (φ<<1), and since z = xsinφ, the following can be written as where β c = 2πn 1 /λ. For a small φ (φ<<1), e iβ c xφ ≈ 1 + i β c x φ . When multiplying both sides of eq. (11) by u * l m , this produces ], eq. (12) can be rewritten as which equals 0 unless m' = m ± 1. While the above results are for a corner bend, they can also be transformed into expressions that are applicable to a continuous bend. (11) Thus, a continuous bend is modeled as a succession of corner bends of Δφ separated by Δz in the propagation direction, and then Δz → 0 as Δφ/Δz remains constant. It can be easily shown that Δφ/Δz = 1/R, where R is the radius of the curvature of the bend. This then allows where Now, eq. (14) can be numerically solved to determine the mode coupling that results from a continuous bend. The key parameters are the coupling constants, given by eq. (16), and the phase mismatch factor Δ (l'm')(lm) given by eq. (15), representing the difference in propagation constant between the modes of the graded index fiber. Numerical values for the constants k (l'm')(lm) and Δ (l'm')(lm) can then be calculated for the given fiber parameters. In this paper, it is assumed that n 1 = 1.46, λ = 1.3 μm, and α = 0.196 /μm. From eq. (16), it can be seen that the coupling constants are inversely proportional to the constant bend radius R. Thus, a decreasing bend radius produces a stronger coupling.
A program is written to solve eq. (14) numerically in order to determine the power distribution between the modes as a function of the distance z for a uniform, continuous bend. It is assumed that all the power is in one mode at z = 0. In addition, the coupling perturbation is proportional to x = rcosθ = r(e iθ + e −iθ )/2. This is a one-dimensional problem as coupling only occurs between a mode with m = m 0 and modes with m = m 0 ± 1, in which the field varies in the direction perpendicular to the bend axis. For wellguided modes (low mode index), the "finite core" mode solutions u(r,θ) should approach the infinite core solutions in ref. 2. Therefore, this is a good way to confirm the validity of the numerical solutions to eqs. (6) and (14). Calculation of the coupling constants, designated k (l'm')(lm) in eq. (16), and the propagation constant differences Δ (l'm')(lm) in eq. (15), followed thereafter.

Results
The above mode coupling analysis was performed numerically to determine the mode power distribution between modes for a uniform continuous bend. All the power was assumed to be initially in the fundamental mode and considered as a one-dimensional problem, since coupling only occurs between modes in which the mode field varies in the direction perpendicular to the bend axis, e.g., between modes in the x-direction for bending around a mandrel with its axis in the y-direction. It was also assumed that the parabolic index variation extended to the radius a, as in a real fiber, rather than to infinity, as in the previous model.
As a prerequisite to presenting the results of the coupled mode calculations, Table 1 presents the number of guided modes with β > n 1 k 0 − dbmax. This table is based on the eigenvalues found in § 2, using the following conversion: dbmax = 0.0712*db2. A small value for dbmax indicates that only the most strongly guided modes are included. All the core modes are included if dbmax = 0.0662. When enumerating the number of modes, the solutions for m = 1, 2, 3, … are counted twice to account for the modes with negative m-values (i.e., twofold degeneracy).
The relative optical power in the fundamental mode as a function of the propagation distance for various bend radii when dbmax = 0.01 is shown in Figs. 2-5. In addition, the figure showing the previous results of Asawa and Taylor cited from Taylor is shown in Fig. 6, which helps in comparing the results shown in Figs. 2-5. Figure 7 shows the number of core modes vs the bend radius R when the relative powers contained in those modes in the first peak were 99, 90, and 50% of the incident power.

Discussion
The difference between the model presented in this paper and the previous model is that the previous model assumes that the quadratic index variation extends to infinity, (2) which means that as β becomes smaller, the mode solutions of both f(r) and β with equally spaced βs become less accurate representations of the true mode solution. With the previous model, periodic power variations are predicted for both large and small bend radii, as in Fig. 6, which presents the calculated power in the fundamental mode of a graded index fiber as a function of the fiber length and long bends with radii of 1 and 0.5 cm. Thus, when using the actual fiber modes in the previous model, small bend radii result in excitation of the modes in a regime where the β spacing is not uniform; in addition, the spatial periodicity in the mode amplitudes can no longer be predicted. This is particularly evident when comparing the results in Figs. 2-4 in the present study with the previous results in Fig. 6. In Fig. 2, with a large bend radius of 3.53 cm, the power oscillates in and out of the fundamental mode with a spatial period of 1.15 mm. This agrees with the predictions of the previous Asawa and Taylor model, as only loworder modes are excited and the "infinite quadratic refractive index variation" is a good approximation for those modes. The results in Fig. 3 with a 0.9413 cm bend radius still agree fairly well with the previous model, yet only 99% of the power returns to the fundamental mode vs 100% predicted using the previous model. Note also that the spatial period of the oscillation is decreased to 0.95 mm. However, for the smaller bend radii of 0.5883 and 0.4243 cm in Figs. 4 and 5, respectively, the behavior is radically different from the previous model, as the power oscillates irregularly owing to the fact that high-order modes with irregular β spacings are excited with the bending of the fiber.

Conclusions
While the predictions of the previous "oversimplified" model are assumed to be valid for large radius bends, this study investigated both the bend radius regime when the previous model breaks down and the mode coupling behavior in this regime. The analysis relied upon numerical (difference equation) solutions of a partial differential equation to (1) obtain mode solutions for a graded index fiber, taking a finite extension of the core and cladding into account, and (2) obtain solutions to the coupled mode equations vs the distance of propagation in a constant-radius bend to obtain the amplitudes of all the fiber modes. The results presented here confirm the predictions of the previous model for bend radii above 1 cm. In this case, when the fundamental mode is initially excited, the power oscillates back and forth between the fundamental mode and high-order modes with a constant spatial period. However, for smaller bend radii, the power tends to oscillate irregularly and never completely returns to the initially excited mode. This behavior is attributed to the excitation of the high-order modes for which the propagation constants are not regularly spaced. Furthermore, in contrast to the previous model, the proposed model provides an easier approach that avoids the need for complicated Bessel functions or quantum mechanics.