I. Introduction

The accelerating development and application of low-loss wideband optical fibers has been accompanied by extensive experimental and theoretical work to characterize their transmissive and dispersive properties.[1],[2] Most of the theoretical analyses have been based on modal expansions of the propagating field.[3] This type of description has the advantages of familiarity and ease of analytic manipulation, but it is difficult to apply quantitatively to multimode fibers with general refractive index profiles. A direct numerical solution of the Helmholtz equation, on the other hand, can provide a detailed and accurate description of the propagating field for a general class of practical fibers and a variety of realistic sources of illumination. A subsequent Fourier analysis of this numerically determined field with respect to axial distance can in turn provide such information relevant to a modal description as the modal eigenvalues or propagation constants and the relative strength of the individual modes present in the propagating field.

In this paper, we describe such a numerical method and give some detailed results applicable to the Corning 1151 fiber[4] with a parabolic graded refractive index profile. The numerical method is based on split operator and discrete Fourier transform techniques, which are known to be exceedingly accurate.[5] Discrete Fourier transform techniques offer an added advantage in providing an accurate dualistic description of the field: the configuration space solution provides detailed information about energy losses and confinement of the beam, while the Fourier transform of the solution gives detailed information about the beam’s angular properties.

Two representative types of illumination are considered here, although this by no means exhausts the possible forms of light coupling that can be treated by the method. In the first case, the fiber is assumed to be uniformly illuminated across its entire cross-sectional area by a beam with uniform intensity and phase. This will be referred to as the case of coherent illumination and corresponds to illumination by a collimated laser beam. In the second case, the initial beam is represented by a 2-D Fourier series in the transverse coordinates with coefficients of equal amplitude and random phase. This will be referred to as the case of incoherent illumination and could represent the illumination of the fiber by an incoherent light source such as an LED.

An obvious practical limitation of the numerical method is that computations can be carried out only over limited propagation distances. It turns out, however, that, for the axially uniform fiber considered in this paper, leaky modes are thoroughly attenuated, and steady-state propagation conditions are established within a fraction of a meter. This can be confirmed by taking the Fourier transform of the field along the fiber axis over a moving window in axial distance z and observing the decay of that portion of the spectrum that corresponds to unbound or leaky modes.

The solution of the Helmholtz equation in terms of operators is discussed in Sec. II, the implementation of this solution in terms of discrete Fourier transforms is treated in Sec. III, and numerical results obtained with the method are discussed in Secs. IV and V. Further discussion of the operator method of solution is provided in the Appendix.

II. Wave Equation and Solution Method

We begin with the usual assumption that the propagation of a single frequency component of light in a fiber waveguide can be described by the scalar Helmholtz equation

where E(ω,x,y,z) is the transverse component of the electric field, ω is the circular frequency of the light, and the refractive index n(ω,x,y) is assumed to depend only on the transverse coordinates x and y.

In this paper we shall be interested only in the propagation of light at a single frequency. A description of pulse dispersion, on the other hand, requires information about a number of frequency components propagating simultaneously. The latter problem will be addressed in a subsequent publication.

The solution to Eq. (1) at z = Δz may be written formally in terms of the field at z = 0 as

where ${\nabla }_{\perp }^2=\left({\partial }^2\right)/\left(\partial x^2\right)+\left({\partial }^2\right)/\left(\partial y^2\right)$. The square root in Eq. (2) can be written in the form

If n in the first right-hand member of Eq. (3) is replaced by the reference value n₀, where n₀ is typically taken to be the value of n in the fiber cladding, Eq. (3) becomes

with

The approximation in Eq. (4) is valid for sufficiently small variations in n(x,y), as shown in the Appendix, and should apply to a wide range of practical fibers. We wish to restrict the solution for E to a single wave propagating in the positive z direction. If the time dependence of E(ω,x,y,z) is exp(iωt), E can thus be expressed in the form

Substituting the expression in Eq. (6) into Eq. (2) and taking the negative sign gives

where

To second order in Δz, Eq. (7) can be rewritten in the symmetrized split operator form

where the error arises from the noncommutation of ${\nabla }_{\perp }^2$ and χ(x,y).[6] The above expression is suitable for generating a numerical solution. Due to the unitarity of the operators in Eq. (9), the solution will be unconditionally stable.

The operation

is equivalent to solving the Helmholtz wave equation

with E(x,y,0) as an initial condition. Therefore, advancing the solution for ℰ(x,y,z) by repeated application of Eq. (9) is equivalent to propagating the beam through a periodic array of thin lenses (see Fig. 1). The first lens is located at z = Δz/2, and the remaining lenses are separated from one another by the distance Δz. Each lens imposes the phase front ϕ(x,y) = Δzχ(x,y) on the beam, and the propagation of the beam between lenses is governed by Eq. (10).

If ${\nabla }_{\perp }^2$ is neglected in comparison with k² in the denominator of Eqs. (7) and (9), one recovers the parabolic or Fresnel approximation. This approximation is valid for small beam divergences and has been used by many workers in a wide range of propagation studies.[5],[7] We have found it to be an excellent approximation to Eq. (9) for steady-state propagation of light in typical multimode optical fibers. In the early stages of propagation, however, plane waves with large angular deviations from the z axis may be present in the beam, and the parabolic approximation can break down. Under these conditions, the solution form in Eq. (9) should still give an accurate description of light propagation. Since a numerical solution is no more difficult to generate with Eq. (9) than it is with the parabolic approximation, Eq. (9) is to be preferred in applications to optical fibers.

III. Solution in Terms of Discrete Fourier Transforms

An accurate numerical representation of Eq. (9) can be obtained by expressing ℰ(x,y,z) as a 2-D Fourier series with a finite number of terms[5]:

where L is the length of the computational grid. Propagation of the beam through a distance Δz in a homogeneous medium of refractive index n₀ transforms ℰ_mn(0), according to Eq. (10), into

where κ_x and κ_y are the transverse wavenumbers

Equation (12) in conjunction with Eq. (11) provides an exact solution to Eq. (10) for an initial field of limited spectral bandwidth. In accordance with Eq. (9) the propagation step is followed by multiplication of ℰby the factor exp(−iΔzχ), whence ℰ(x,y,Δz) becomes

If the spectrum of ℰ′ remains finite and is bounded by that of ℰ, the Fourier coefficients of ℰ′ can be evaluated exactly in terms of the sampled values

where jΔx,lΔy are points on the computational grid, and there will be a one-to-one correspondence between the Fourier coefficients ${\mathcal{E}}_{nm}^{\prime }$ and the elements of the discrete Fourier transform[8]

The latter can be computed with the well-known FFT algorithm.

Following the next propagation step the numerical representation of Eq. (9) then remains exact. Thus, if the spectrum of ℰ(x,y,z) remains finite, it is possible to generate an exact numerical representation of Eq. (9). The spectral bandwidth of ℰ(x,y,z) is, in practice, never perfectly finite, but, for most optical fiber studies, it is possible to set up a configuration space computational grid with sufficient resolution to keep the spectral power on the boundaries of the corresponding wavenumber space grid extremely small. Spectral power on the mesh boundaries is normally monitored, making it possible to confirm the accuracy of a given calculation.

The angle between the direction of a representative plane wave with transverse wavevector (κ_x, κ_y) and the z axis is given by

The value of N in Eq. (16) will be determined by L and θ_max, the maximum value of θ for a beam propagating in the fiber core. The value of L will be of the order of the fiber diameter (including core and cladding), and ${\theta }_{\text{max}}\approx {\left(n_{\text{max}}^2-n_0^2\right)}^{1/2}$. Note that θ_max is the N.A. for a step index fiber. The minimum spatial bandwidth for ℰ required to accommodate the steady-state field is thus defined by the relations

IV. Field Properties for a Realistic Fiber

The solution method described in Sees. II and III has been applied to the Corning 1151 fiber, which has a 125-μm outer diameter and a 62.5-μm core diameter (see Fig. 2). The refractive index of the fiber as a function of radius r is described by[9]

where a is the core radius, n₀ = 1.5, and Δ = 0.008.

Calculations were performed for two forms of ℰ(x,y,0). In the coherent illumination case, ℰ was taken to be constant in amplitude and phase over the entire fiber cross section. In the incoherent illumination case, the coefficients in the Fourier series in Eq. (11) were selected with random phases and equal amplitudes. For computational efficiency the field was assumed to be symmetric with respect to reflections about either the x or y axes. This made possible a cosine series solution for Eq. (11).

The wavelength of the light was taken to be 1 μm, and the upper right quadrant of the fiber cross section was represented by a 64 × 64 computational grid with Δx = Δy = 0.98 μm, creating a maximum possible spatial bandwidth of Δk_max = 2π/Δx = 6.41 × 10⁴ cm⁻¹. In the case of the incoherent beam, the initial spectrum was filtered, so that only half of this bandwidth was used. The axial space increment Δz was taken to be 10 μm. In all cases the return by reflection of radiation that has reached the outer fiber boundary was prevented by placing a strong absorber on the outer circumference of the fiber.

Calculations were carried out for propagation distances of 18 cm and 20 cm for the coherent and incoherent illumination cases, respectively. In the incoherent illumination case, a steady state appeared to be established within a few centimeters, whereas a distance of the order of 15 cm was required to establish a comparable steady state for the coherently illuminated case, due to significant initial excitation of leaky but almost guided modes.

Figures 3–8 apply to the case of coherent illumination. Figure 3 shows the axial intensity of light as a function of axial distance over two separate 1-cm path segments; the first extending from 0 to 1 cm, the second from 17 cm to 18 cm. An almost periodic focusing and defocusing pattern is set up at once with the foci separated by a constant distance of 0.78 mm. In Fig. 3(b) the pattern has reached steady-state behavior, but a certain amount of amplitude modulation persists. Both the small scale and the large scale features of the amplitude modulation are displayed in Fig. 4, which is plotted for the first 15 cm of propagation. This behavior should be compared with that of a Gaussian beam launched in a quadratic lenslike medium. As is well known,[10] the Gaussian beam retains its shape and focuses with a spatial period that depends solely on the curvature of the refractive index. For the fiber under consideration, this period would be 0.776 mm, which is close to the period exhibited in Figs. 3(a) and 3(b). It should be emphasized, however, that perfect periodicity in the focusing pattern can be expected only for a beam that is initially Gaussian and that propagates in an infinite square law refractive medium. The intermediate peaks and amplitude modulation in Figs. 3 and 4 are thus due to the non-Gaussian initial beam shape and the finite fiber core diameter.[11]

Spatial (r,z) contours representing the radii of circles containing a specific percentage of the local beam power as a function of z are exhibited in Figs. 5(a) and 5(b) for the path segments corresponding to Figs. 3(a) and 3(b). In the order of increasing radius these contours correspond to 20%, 40%, 60%, and 80%. The contours for the 0–1-cm segment clearly indicate the presence of significant power in the cladding. The same contour set, however, is well contained within the core in Fig. 5(b).

Figure 6 shows spectral (|κ|,z) contours representing the radii of circles in transverse wavenumber space that contain the same percentages of the beam power. The path segments are the same as for Figs. 3 and 5. The spectral contours may also be interpreted as angular contours by means of Eq. (17). A corresponding scale of θ in degrees is provided on the right-hand vertical axes in Fig. 6.

Although the contour patterns in Fig. 6(b) are more regular than those exhibited in Fig. 6(a), the former do not represent significantly different beam divergence patterns, because the beam was collimated to begin with. The maximum excursion in θ for the 80% power contour in Fig. 6(b) is about 6°. This corresponds to a N.A. of 0.16. The manufacturer’s quoted value for the 90% power N.A. is 0.16 ± 0.02.[4] The simple formula[12] $\mathrm{N}.\mathrm{A}.={\left(n_1^2-n_2^2\right)}^{1/2}$ for a step fiber also gives a N.A. value of 0.16. Thus, Fig. 6 is in rough agreement with simple theory, but it shows a great deal more detail.

The fractional power gained or lost by the fiber core in the coherent illumination case is shown in Fig. 7. Initially the power in the core is seen to double as energy streams into the core from the cladding. However, this additional energy slowly leaks back out of the core over many bounces of the wavefront, and, by the time the steady state is established, the core power is only slightly above its initial value. A calculation of power in the cladding (not shown) reveals that by 20 cm the value fluctuates about an average of 0.4%.

The uncertainty product ΔrΔκ over a 1-cm distance is exhibited in Fig. 8 as a function of axial distance z, where again Δr and Δκ are the radii containing the fractions 0.2, 0.4, 0.6, 0.8 of the local beam power. The most obvious feature of Fig. 8 is the fluctuating behavior of the uncertainty product. A Gaussian shaped beam would have a constant uncertainty product expressible in terms of the power fraction f as ΔrΔκ = −ln(1 − f). For f = 0.8 this expression would give a value of 1.61. For the same value of f, Fig. 8 shows values of the uncertainty product that exceed this value by at least an order of magnitude. (The ratio of the uncertainty product to the corresponding value for a Gaussian beam gives the beam quality or the number of times diffraction limited.)

Figures 9–12 refer to the case of incoherent illumination. Figure 9 shows the fractional power radius as a function of z over the first centimeter of propagation. This should be compared with Fig. 5(a) corresponding to coherent illumination. Clearly a steady state is approached in a much shorter distance in the incoherent illumination case due to the presence initially of a large proportion of the total beam energy in large angle waves, which are quickly lost from the core.

The peak excursions in radius are similar to those exhibited in Fig. 3, but the minimum radii attained for the incoherent illumination case are larger than those attained for the incoherent illumination case. Clearly, the incoherent beam does not focus as sharply as the coherent beam, even though the two beams are similarly confined.

The spectral or angular contours for the incoherent illumination case are shown for the interval from 0 to 1 cm in Fig. 10, which should be compared with Fig. 5. The maximum angular excursions in the two cases are close, although initially the incoherent beam contains substantially wider angles. The incoherent beam also clearly shows evidence of having reached a steady state within a fraction of the centimeter in agreement with Fig. 9. The angular spread of the source, it will be noted, is roughly a factor of 2 greater than the N.A. of the fiber. The corresponding uncertainty product is displayed in Fig. 11. The excursions corresponding to a given power fraction are higher in Fig. 11 than they are in Fig. 7.

Figure 12 shows the fractional change in the core power for the incoherent illumination case. Initially the core gains power, but, after approximately 0.5 mm of propagation, it begins to lose power rapidly. The core power reaches a steady state after approximately 0.5 cm at which point it contains less than 40% of its original energy.

V. Axial Spectrum Behavior

Fourier transformation of the field ℰ(x′,y′,z) with respect to z for a fixed transverse position (x′,y′) will reveal the normal mode eigenvalues or propagation constants β_n, where ℰ(x,y,z) is expressible in terms of the modal eigenfunctions u_n(x,y) as

By transforming over a moving window in z, one can observe changes in the modal content of the field as a function of propagation distance z. The amplitudes of the guided modes or the modes that correspond to bound states in quantum mechanics will remain constant, but the amplitudes of leaky modes corresponding to the continuum in quantum mechanics will decay. The establishment of a propagational steady state can be determined when the amplitudes of the leaky modes have been observed to decay sufficiently. Figures 13–15 show spectra of the axial field ℰ(0,0,z) for both the coherent and incoherent illumination cases calculated using discrete Fourier transforms for windows centered at different axial positions. To avoid spectral aliasing problems that arise from transforming a record of finite length, the sampled values have been multiplied by the Hanning truncation function.[13]

where Z represents the length of the window. In all cases the widths of the spectral peaks are determined solely by the sample length.

In general the magnitudes of propagation constants for the guided modes will be bounded according to 0 < |β_n| < kΔ. To insure that the axial spectrum is accurately represented it is thus necessary for the axial sampling distance Δz to satisfy Δz < π(kΔ)⁻¹. This condition is satisfied in the present example by a factor of 5.

Figure 13 shows the spectrum obtained for the coherent illumination case by transforming over the first 2.56 cm of propagation path. The propagation constants β_n, expressed relative to the value k, can be read from the positions of the spectral peaks. The peaks to the left of 0 represent the guided or trapped modes of the fiber, while spectral components to the right of 0 represent leaky modes. (In analogy with quantum mechanics the value β = 0 corresponds to the top of the potential well.) The positions of the guided mode eigenvalues[14] correspond closely to the eigenvalues of the Schrödinger equation for the 2-D harmonic oscillator potential:

The propagation constants defined by the spectral peaks in Fig. 13 are compared with values calculated from Eq. (22) in Table I. Because the axis is a center of symmetry, only the even parity modes are exhibited.

The presence of a sharp peak just to the right of 0 in Fig. 13 is due to the near coincidence of one of the eigenvalues in Eq. (22) and the top of the harmonic oscillator potential well. This peak is analogous to the virtual level that exists for the deuteron and that leads to the low lying resonance in the scattering of neutrons by protons.[15] Evidence that the above peak corresponds to a virtual state (leaky mode) is seen in Fig. 14, which shows the axial spectrum for a window extending from z = 16.13 cm to z = 18.58 cm. The peaks representing the trapped modes have the same relative height as in Fig. 12, but the peak near β = 0 has nearly vanished. The indicated exponential decay length for the peak near β = 0 is ~4 cm. The overwhelming preponderance of power in the trapped modes shown is evidence that, for practical purposes, the field is in a steady state.[16]

Axial spectra for the case of incoherent illumination, calculated for windows extending from 0.45 cm to 1.05 cm and from 0.96 cm to 5.87 cm, are shown in Figs. 15(a) and 15(b), respectively. The excitation of the virtual state mode is evident for this case as well, but the fraction of the total power that it contains is substantially less than it was in the coherent illumination case. Consequently, a steady state is well established within a propagation distance of about 5 cm. As would be expected, the propagation constants contained in the spectra of Figs. 15(a) and 15(b) agree with those of Fig. 14 to within the resolution of the spectral peaks (2.5 cm⁻¹). Note the random contribution of the various guided modes to the axial spectrum.

VI. Summary and Conclusions

We have described an accurate and general numerical method for solving the Helmholtz equation for the electric field in an optical fiber. The method has been applied to determine the detailed properties of the field for a practical graded-index fiber. The principal conclusion to be drawn from these results is that weakly guided modes decay within a short propagation distance, allowing steady-state propagation conditions to be set up within centimeters.

Appendix: Further Discussion of the Solution Method for the Helmholtz Equation

We wish to find a solution to the Helmholtz equation

(A1)$$\frac{{\partial }^{2}E}{\partial z^2}=-{\nabla }_{\perp }^{2}E-\frac{{\omega }^2}{c^2}{n}^2\left(x,y\right)E,$$

that has the form

(A2)$$E\left(x,y,z\right)=\text{exp}\left(\pm iAz\right)E\left(x,y,0\right),$$

where the application of the operator exp(±iAz) is defined by a Taylor series expansion

(A3)$$\text{exp}\left(\pm iAz\right)=1\pm iAz-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.A^2{z}^2+\dots .$$

[Equation (A3) may also be viewed as a perturbation series representation of the solution to Eq. (A1).] Substitution of the expression in Eq. (A2) into Eq. (A1) gives the following relation that must be satisfied by A:

(A4)$$A^2={\nabla }_{\perp }^2+\frac{{\omega }^2}{c^2}{n}^2={\nabla }_{\perp }^2+k^2{\left(1+\frac{\delta n}{n_0}\right)}^2,$$

where k = n₀ω/c and n = n₀ + δn(x,y).

We can use Eq. (A4) to test the accuracy of the following approximate forms for the solution operator A:

(A5a)$$A_1=\frac{{\nabla }_{\perp }^2}{2k}+k\left(1+\frac{\delta n}{n_0}\right),$$

(A5b)$$A_2={\left({\nabla }_{\perp }^2+k^2\right)}^{1/2}+k\delta n/n_0.$$

The operator A₁ when used with Eq. (A2) gives a solution to the Fresnel or parabolic approximation of the Helmholtz equation. The operator A₂ is equivalent to the expression in Eq. (4). We shall refer to the solution obtained by substitution of A₁ into Eq. (A2) as the extended Fresnel or parabolic approximation. Squaring A₁ gives

(A6)$$A_1^2=\left(1+\frac{\delta n}{2n_0}\right){\nabla }_{\perp }^2+{\nabla }_{\perp }^2\frac{\delta n}{2n_0}+\frac{{\nabla }_{\perp }^4}{4k^2}+k^2{\left(1+\frac{\delta n}{n_0}\right)}^2.$$

We can judge the accuracy of a solution obtained with A₁ by assuming that δn = δn₀ = const and that

(A7)$$E\left(x,y,0\right)=E_0\text{exp}\left[i\left({\kappa }_{\perp }·{\mathit{r}}_{\perp }\right)\right],$$

where r_⊥ = (x,y) and E₀ is a constant. Applying A² [the expression in Eq. (A4)] is equivalent to writing

(A8)$$A^2={-\kappa }_{\perp }^2+k^2{\left(1+\frac{\delta n^0}{n_0}\right)}^2,$$

while applying $A_1^2$ to Eq. (A7) is equivalent to writing

(A9)$$A_1^2=-\left[1+\frac{\delta n^0}{n_0}-{\left(\frac{{\kappa }_{\perp }}{2k}\right)}^2\right]{\kappa }_{\perp }^2+k^2{\left(1+\frac{\delta n^0}{n_0}\right)}^2.$$

Comparing Eq. (A9) with Eq. (A8), we see that Eq. (A9) is a good approximation to Eq. (A8) provided

(A10)$$\begin{array}{cc}\frac{\delta n^0}{n_0}\ll 1,& \frac{{\kappa }_{\perp }}{k}\ll 1.\end{array}$$

Thus we conclude that the Fresnel approximation is a good one for small refractive index changes and waves whose directions do not deviate much from the z axis. These conditions are, of course, well known.

We now square A₂ to obtain

(A11)$$\begin{array}{ll}{A}_2^2=\hfill & {\nabla }_{\perp }^2+k^2+k{\left({\nabla }_{\perp }^2+k^2\right)}^{1/2}\frac{\delta n}{n_0}\hfill \\ \hfill & +k\frac{\delta n}{n_0}{\left({\nabla }_{\perp }^2+k^2\right)}^{1/2}+k^2{\left(\frac{\delta n}{n_0}\right)}^2.\hfill \end{array}$$

If use of the expansion

(A12)$${\left({\nabla }^2+k^2\right)}^{1/2}=k\left(1+\frac{{\nabla }_{\perp }^2}{2k^2}-\frac{1}{8}\frac{{\nabla }_{\perp }^4}{k^4}+\dots \right)$$

is made in Eq. (A11), one obtains

(A13)$$\begin{array}{ll}{A}_2^2=\hfill & \left(1+\frac{\delta n}{2n_0}\right){\nabla }_{\perp }^2+{\nabla }_{\perp }^2\frac{\delta n}{2n_0}-\frac{1}{8}\frac{{\nabla }_{\perp }^4}{k^2}\frac{\delta n}{n_0}-\frac{1}{8}\frac{\delta n}{n_0}\frac{{\nabla }_{\perp }^4}{k^2}\hfill \\ \hfill & +k^2{\left(1+\frac{\delta n}{n_0}\right)}^2+\dots .\hfill \end{array}$$

Without any loss of generality we can omit from further discussion any higher order terms in ${\nabla }_{\perp }^2$.

For δn = δn⁰ application of Eq. (A13) to the expression in Eq. (A7) is equivalent to writing

(A14)$$A_2^2=-\left[1+\frac{\delta n^0}{n_0}+\frac{1}{4}\frac{\delta n^0}{n_0}{\left(\frac{{\kappa }_{\perp }}{k}\right)}^2\right]{\kappa }_{\perp }^2+k^2{\left(1+\frac{\delta n^0}{n_0}\right)}^2.$$

Equation (A14) is clearly a good approximation to Eq. (A8) if δn⁰/n₀ ≪ 1. However, unlike the case of Eq. (A9), there is no restriction on the size of κ_⊥ except that it be less than k.

In conclusion, the extended Fresnel approximation permits operator splitting, which is crucial to the use of spectral methods of solution, but it is not restricted to waves or beams that propagate at small angles to the z axis.

The authors acknowledge the valuable assistance of P. F. Thompson in the development of the postprocessor program used to digest the large volume of data generated by these calculations.

This work was performed under the auspices of the U.S. Department of Energy by the Lawrence Livermore Laboratory under contract W-7405-ENG-48.

Figures and Table

 

Fig. 1 Algorithm for solving the Helmholtz equation replaces fiber by a system of lenses. In between lenses field satisfies the Helmholtz equation for a homogeneous medium. Algorithm treats large angle waves accurately and is accurate for small refractive index gradients.

Download Full Size | PDF

 

Fig. 2 Geometry of calculation showing Corning 1151 fiber. At the edge of the cladding, a strong absorber is placed to prevent waves from reentering cladding or reflections from computational mesh boundary.

Download Full Size | PDF

 

Fig. 3 On-axis intensity as a function of axial distance for the case of coherent illumination. Shown are intensity plots over two separate 1-cm path segments.

Download Full Size | PDF

 

Fig. 4 On-axis intensity as a function of axial distance for the coherent illumination case over first 15 cm of propagation.

Download Full Size | PDF

 

Fig. 5 Fractional power contours over two propagation ranges for coherent illumination case. Plotted as functions of axial distance are radii of circles containing 20%, 40%, 60%, and 80% of total beam power at the given z position.

Download Full Size | PDF

 

Fig. 6 Fractional power contours in transverse wavenumber space. Plotted are spectral radii containing specific fractions of spectral power as functions of propagation distance z. Fractions are the same as for Fig. 5. Note that the spectral radii are expressed both in terms of wavenumber and angle in degrees.

Download Full Size | PDF

 

Fig. 7 Percent core power change as a function of propagation distance for the case of coherent illumination.

Download Full Size | PDF

 

Fig. 8 Uncertainty product as a function of position for coherent illumination and power fractions 0.2, 0.4, 0.6, and 0.8. Curves with increasing radius indicate increasing power fraction. The uncertainty product for a Gaussian beam would be constant for a fixed power fraction.

Download Full Size | PDF

 

Fig. 9 Fractional power contours for the incoherent illumination case. Power fractions are 0.2, 0.4, 0.6, and 0.8.

Download Full Size | PDF

 

Fig. 10 Fractional power contours in transverse wavenumber space for the incoherent illumination case. The same information is plotted as in Fig. 6.

Download Full Size | PDF

 

Fig. 11 Uncertainty product ΔrΔκ over first cm of propagation. Power fractions are 0.2, 0.4, 0.6, and 0.8.

Download Full Size | PDF

 

Fig. 12 Percent change in core power as a function of propagation distance for incoherent illumination case.

Download Full Size | PDF

 

Fig. 13 Axial spectrum of coherently excited field over first 2.56 cm of propagation path.

Download Full Size | PDF

 

Fig. 14 Axial spectrum of coherently excited field over window extending from 16.13 cm to 18.58 cm.

Download Full Size | PDF

 

Fig. 15 Axial spectrum of incoherently excited field (a) over window extending from 0.45 cm to 1.06 cm, and (b) over window extending from 5.26 cm to 5.87 cm.

Download Full Size | PDF

Table I. Comparison of Modal Eigenvalues for an Infinite Square Law Refractive Index Medium and an Optical Fiber with a Parabolic Graded-Index Profile ^a

View Table

References

1. For a comprehensive survey of research in this field see, for example, Optical Fiber Technology, D. Gloge, Ed. (IEEE Press, New York, 1976).

2. See also J. A. Arnaud, Beam and Fiber Optics (Academic, New York, 1976).

3. Examples of modal theory can be found in D. Gloge and E. A. J. Marcatili, Bell Syst. Tech. J. 52, 1563 (1973); D. Gloge, IEEE Trans. Microwave Theory Tech. MIT-23, 106 (1975); [CrossRef]  D. Marcuse, Theory of Dielectrical Optical Waveguides (Academic, New York, 1964); and R. Olshansky and D. B. Keck, Appl. Opt. 15, 483 (1976) [CrossRef]   [PubMed]  .

4. EVA Buffered Corguide Fibers Product Bulletin No. 2 (Telecommunication Products Dept., Corning Glass Works, Corning, N.Y. 14830, 1 May 1976).

5. See, for example, J. A. Fleck Jr., J. R. Morris, and M. D. Feit, Appl. Phys. 10, 129 (1976) [CrossRef]  .

6. Operator splitting results in a separation of the propagation and phase updating parts of the calculation. Since the operators inside the bracket in Eq. (7) do not commute, splitting must involve an approximation that holds for limited propagation distances Δz. By symmetrizing the splitting, the commutation error is further reduced. For a more complete discussion of this question see the Appendix of Ref. [5].

7. For a comprehensive review of the application of the Fresnel approximation to problems in nonlinear optics see, for example, J. H. Marburger, Prog. Quantum Electron. 4, 35 (1975) [CrossRef]  .

8. This follows from the sampling theorem; see, for example, E. O. Brigham, The Fast Fourier Transform (Prentice-Hall, Engle-wood Cliffs, N.J., 1974), pp. 99–102.

9. It is more usual to write n as n = n₀[1 − Δ(r/a)²] rather than in the form of Eq. (19). However, in the present development, the peak refractive index plays no role whereas the cladding index does.

10. See, for example, P. K. Tien, J. P. Gordon, and J. R. Whinnery, Proc. IEEE 53, 129 (1965) [CrossRef]  .

11. It has been shown that beam aberration and an erratic focusing pattern are general consequences of propagation in a nonquadratic lenslike medium. M. D. Feit, J. A. Fleck Jr., and J. R. Morris, J. Appl. Phys. 48, 3301 (1977) [CrossRef]  .

12. S. E. Miller, E. A. J. Marcatili, and T. Li, Proc. IEEE 61, 1703 (1973) [CrossRef]  .

13. E. O. Brigham, Proc. IEEE 61, 141 (1973).

14. We have expressed the modal z dependences in the form exp(−ikz + iβ_nz) in order to emphasize the similarity between the waveguide problem and the quantum mechanical problem of a particle in a potential well.

15. J. M. Blatt and V. F. Weisskoff, Theoretical Nuclear Physics (Wiley, New York, 1952), p. 64.

16. The short decay length for the leaky mode radiation may seem surprising in the light of previous discussion of the subject. See, for example, A. W. Snyder, Appl. Phys. 4, 273 (1974); and [CrossRef]  A. W. Snyder and D. J. Mitchell, J. Opt. Soc. Am. 64, 599 (1974) [CrossRef]  .