Modeling of transient modal instability in fiber amplifiers

Benjamin G. Ward

doi:10.1364/OE.21.012053

1. Introduction

Recent evidence suggests that dynamic thermal modal instability (TMI) affects the beam quality of large mode area high average power fiber amplifiers [1–5]. The resulting degradation severely limits the usefulness of these devices. Efforts to improve fiber amplifier designs to mitigate this degradation benefit from numerical models that reveal the effects of amplifier parameters on performance. Several previous investigations have isolated the cause to refractive index gratings caused by spatial-temporal temperature variations that arise from quantum defect heating through the thermo-optic effect [6–10]. A key aspect of the physical description of this phenomenon is the time behavior of the temperature profile within the fiber. Under the assumption that the temperature at any point within the fiber is a periodic function of time, a steady-state solution to the optical wave equation arises that describes the process of stimulated thermal Rayleigh scattering (STRS). Numerical models have been presented that capture the essential features of STRS in fiber amplifiers [9–12]. Observations of periodic modal fluctuations in the output beams of amplifiers operating near the TMI threshold are consistent with the STRS picture of TMI [3,7]. This regime of TMI relies on the amplifier optical and temperature fields achieving a dynamic equilibrium which may be unstable with respect to external perturbation [9].

Although significant understanding is provided by the steady-state model, several factors motivate further study of the transient behavior of TMI. Amplifiers operated with a rapid turn-on may not reach steady-state equilibrium prior to the onset of TMI thus the threshold calculated using the periodic method may be higher than values that characterize realistic operating conditions. One definition of the TMI threshold has been given in terms of the buildup time required to observe the effect [5]. Also, observed TMI is often characterized by chaotic behavior [3,4,7] not able to be described in the context of a periodic model. The goal of the work presented here is to present a computational model that can enable a greater understanding of the chaotic and transient regimes of TMI. This model incorporates a 3D beam propagation method (BPM), including laser gain, for the spatial evolution of the signal field as well as a 3D time-dependent thermal model. This model is also capable of describing the steady-state periodic case.

This paper begins with a description of the BPM which is of a new type, followed by a description of the thermal solver, also of a new type. Example results are then presented and discussed.

2. Hybrid finite element-harmonic beam propagation method

The problem at hand requires the beam propagation calculation to be performed over the entire length of the fiber thousands of times, once for each time step. Each propagation step must occur over a fraction of the inter-modal beat length of the fiber leading to a requirement to evaluate upwards of one million BPM steps. Thus the computational speed of the BPM is at a premium. Recently, a new type of BPM was presented, the azimuthal harmonic expansion beam propagation method [13], which enables much faster computation than either the split-step fast Fourier transform BPM (FFT-BPM) or finite difference BPM (FD-BPM). This method is based on a hybrid discretization in cylindrical coordinates that combines a finite-difference discretization in the radial coordinate with a harmonic Fourier expansion in the azimuthal angle. However, a 1D finite element approach in the radial direction further aids treatment of fiber designs with discontinuous refractive index profiles, such as photonic crystal fibers, as well as enabling variable grid spacing to conserve the total number of unknowns in a natural fashion. A simplification of such a method has been previously applied to the study of stimulated Brillouin scattering in fibers [14].

The starting point is the scalar optical wave equation Eq. (1) valid for weakly-guiding fibers

\nabla^{2} E (\vec{r}) + n^{2} (\vec{r}) k_{0}^{2} E (\vec{r}) = 0,

where

k_{0} = ω / c

is the free space wave-vector and

n (\vec{r})

is the refractive index distribution. Invoking the slowly-varying envelope and paraxial approximations yields Eq. (2)

[- 2 i β \frac{\partial}{\partial z} + \nabla_{t}^{2} + n^{2} (r, φ, z) k_{0}^{2} - β^{2}] E (r, φ, z) = 0

where

E (r, φ, z)

now refers to the slowly-varying envelope,

β

is the wave-vector describing rapid oscillations in the propagation direction, and

\nabla_{t}^{2}

is the transverse Laplacian operator. Setting

\partial E / \partial z = 0, \partial n / \partial z = 0

yields equation Eq. (3) for the stationary propagating modes of the fiber

[\nabla_{t}^{2} + n^{2} (r, φ) k_{0}^{2} - β^{2}] E (r, φ) = 0.

To obtain a finite element solution of this equation the variational form Eq. (4) is introduced

S = \frac{1}{2} \iint_{Ω} [- {| \nabla_{t} E (r, φ) |}^{2} + n^{2} (r, φ) k_{0}^{2} {| E (r, φ) |}^{2} - β^{2} {| E (r, φ) |}^{2}] d A

where

Ω

represents the fiber cross-section. Expanding as shown in Eq. (5)

E (r, φ) = \sum_{q = - Q}^{Q} E_{q} (r) \exp [i q φ]

where

Q

is the truncation order of the series leads to Eq. (6)

\nabla_{t} E (r, φ) = \frac{\partial E}{\partial r} \hat{r} + \frac{1}{r} \frac{\partial E}{\partial φ} \hat{φ} = \sum_{q = - Q}^{Q} [\frac{\partial E_{q}}{\partial r} \hat{r} + \frac{i q}{r} E_{q} \hat{φ}] \exp [i q φ]

where Eq. (7) and Eq. (8) define the refractive index in Eq. (4),

n^{2} (r, φ) = n_{0}^{2} + \sum_{q = - Q}^{Q} V_{q} (r) \exp [i q φ]

V_{q} (r) \equiv \frac{1}{2 π} \int_{0}^{2 π} δ n^{2} (r, φ) \exp [- i q φ] d φ

and

n_{0}

is a suitable background refractive index. With these definitions the variational form is given by Eq. (9)

S = - \frac{1}{2} \sum_{q = - Q}^{Q} \iint_{Ω} [{| \frac{\partial E_{q}}{\partial r} |}^{2} + (\frac{q^{2}}{r^{2}} - k_{0}^{2} n_{0}^{2} + β^{2}) {| E_{q} |}^{2} - k_{0}^{2} \sum_{q^{'} = - Q}^{Q} V_{q - q^{'}} E_{q}^{*} E_{q^{'}}] d A .

To accomplish the radial integral it is convenient to divide the integration domain into 1D finite elements represented by

N - 1

line segments with end points

r_{k} : k \in 1, N

. The integrals over the fiber cross section are then approximated as shown in Eq. (10) and Eq. (11)

\iint_{Ω} f (r) \exp [i q φ] d A \approx π δ (q) \sum_{k = 1}^{N - 1} f (\frac{r_{k} + r_{k + 1}}{2}) (r_{k + 1}^{2} - r_{k}^{2})

{\iint_{Ω} [f^{'} (r)]}^{2} \exp [i q φ] d A \approx π δ (q) {\sum_{k = 1}^{N - 1} (\frac{f (r_{k + 1}) - f (r_{k})}{r_{k + 1} - r_{k}})}^{2} (r_{k + 1}^{2} - r_{k}^{2}) .

The variational form is then conveniently expressed using matrix notation as shown in Eq. (12)

S = - \frac{1}{2} E^{*} \cdot [K - β^{2} M] \cdot E

the minimization of which leads to the generalized eigenproblem given in Eq. (13)

[K - β^{2} M] \cdot E = 0

where

E

is a vector of dimension

(2 Q + 1) N

and

K

and

M

are square matrices with the same dimension. If now the envelope is allowed to vary slowly along the length of the fiber such that

\partial E / \partial z \neq 0

, then the longitudinal wave vector value

k

is freely chosen leading to the matrix form of the paraxial wave equation as shown in Eq. (14)

[- 2 i k M \frac{d}{d z} + K - k^{2} M] \cdot E = 0.

Differencing

E

with respect to

z

and taking one explicit and one implicit step each of distance

Δ z / 2

results in the usual Crank-Nicholson update rule for propagation given in Eq. (15)

E_{m + 1} = {(M + \frac{i Δ z}{4 k} (K_{m + 1} - k^{2} M))}^{- 1} (M - \frac{i Δ z}{4 k} (K_{m} - k^{2} M)) E_{m}

where the cumulative propagation distance is

z = (m - 1) Δ z : m \in 1, M

. In this scheme, the waveguide refractive index profile, the gain within the core, and the thermally-induced refractive index profile are incorporated as shown in Eq. (16)

δ n^{2} (r, φ, z) = δ n_{wg}^{2} (r, φ) + i n_{0} g (r, φ, z) + 2 n_{0} \frac{d n}{d T} Δ T (r, φ, z)

where

δ n_{wg}^{2}

is the waveguiding index profile,

g

is the optical gain,

d n / d T

is the thermo-optic coefficient, and

Δ T

is the local temperature rise. To a good approximation, the gain is instantaneously determined by the local pump and signal intensities while the local temperature rise depends on the past history of heat deposited in the neighborhood of each point within the fiber according to the time-dependent heat equation.

The matrix $K$ changes at each step due to its dependence on the signal intensity profile. To keep the gain contribution space-centered during the propagation step, the field is updated first using the initial matrix $K_{m}$ for both the explicit and implicit steps. Then these field values are used to calculate $K_{m + 1}$ which is then used to repeat the implicit step over the second half of the propagation interval. Furthermore, the average refractive index change can increase significantly over the length of the fiber necessitating resetting the reference refractive index between steps. This is done automatically in the computer implementation.

3. Gain and heat deposition

The following discussion pertains to active media for which the atomic energy level populations are determined by two-level rate equations where instantaneous decay from the lower level to the ground state is assumed. Under these conditions the local fraction of ions in the upper state is given by Eq. (17)

f (I_{s}, I_{p}) = (\frac{I_{s} σ_{a s} τ}{h ν_{s}} + \frac{I_{p} σ_{a p} τ}{h ν_{p}}) / (1 + \frac{I_{s} (σ_{a s} + σ_{e s}) τ}{h ν_{s}} + \frac{I_{p} (σ_{a p} + σ_{e p}) τ}{h ν_{p}})

where

I_{s, p}

are the signal and pump intensities,

h ν_{s, p}

are the energies of each signal and pump photon,

σ

are the absorption and emission cross-sections at the pump and signal wavelengths, and

τ

is the upper-state lifetimes. The gain is then given by Eq. (18)

g (r, φ, z) = N_{ion} [f (r, φ, z) σ_{e s} - (1 - f (r, φ, z) σ_{a s})]

where

N_{ion}

is the doping density. Similarly, the pump gain is given by Eq. (19)

g_{p} (r, φ, z) = N_{ion} [f (r, φ, z) σ_{e p} - (1 - f (r, φ, z) σ_{a p})] .

For cladding-pumped amplifiers, the paraxial approximation is generally not valid for the pump light guided in the cladding due to its large numerical aperture. Therefore, constant pump intensity throughout the cladding is assumed to be maintained through the process of mode mixing so that the evolution of the pump intensity is given by Eq. (20)

\frac{d}{d z} I_{p} (z) = \pm \frac{1}{A_{clad}} \iint_{Ω} g_{p} (r, φ, z) I_{p} (r, φ, z) d A

where

A_{clad}

is the area of the pump cladding and the

\pm

represents counter versus co-pumping arrangements. The heat load due to the quantum defect neglecting a small contribution from spontaneous emission is given by Eq. (21)

Q (r, φ, z) = (\frac{ν_{p}}{ν_{s}} - 1) g (r, φ, z) I_{s} (r, φ, z) .

The normalization of the optical fields is chosen so that

I_{s} = {| E |}^{2}

. Starting from some known temperature distribution, the temperature evolves in time at each spatial point according to the time-dependent heat equation given by Eq. (22)

\frac{\partial}{\partial t} T (r, φ, z, t) = \frac{K}{ρ C} \nabla^{2} T (r, φ, z, t) + \frac{1}{ρ C} Q (r, φ, z, t)

where

K

is the thermal conductivity,

ρ

is the mass density, and

C

is the heat capacity within the fiber. This equation is amenable to a solution approach similar to that used for the optical propagation equation. Both

T

and

Q

are approximated in

r

and

φ

using the 1D finite element scheme in the radial direction and the azimuthal harmonic expansion. The longitudinal component of the Laplacian is evaluated using finite differences in the temperature values at the

z

points defining the optical beam propagation steps as shown in Eq. (23)

{\frac{\partial^{2}}{\partial z^{2}} T (r, φ, z) |}_{m} \approx \frac{1}{Δ z^{2}} [T_{m + 1} (r, φ) - 2 T_{m} (r, φ) + T_{m - 1} (r, φ)] .

To a good approximation, the thermal properties of the fiber including conductivity and heat capacity are constant throughout the cross-section. Significantly, this causes the different azimuthal harmonic terms to de-couple yielding

(2 Q + 1)

independent Crank-Nicholson time step update equations for the temperature distribution as shown in Eq. (24)

T_{q, l + 1} = {[M_{TH} + \frac{Δ t}{2} \frac{K}{ρ C} K_{q,TH}]}^{- 1} \cdot [(M_{TH} - \frac{Δ t}{2} \frac{K}{ρ C} K_{q,TH}) \cdot T_{q, l} + \frac{Δ t}{2 ρ C} M_{TH} \cdot (Q_{q, l + 1} + Q_{q, l})]

where

l

is the time step index,

Δ t

is the time step,

q

is the harmonic component,

T_{q, l}

are vectors of length

M \cdot N

and

M_{TH}, K_{q, TH}

are sparse matrices of this same dimension. These sparse matrices are assembled by inserting blocks of values that when multiplied by the

T_{q, l}

approximate the integrals over the fiber cross-section as described by Eq. (25) and Eq. (26)

M_{TH} = \sum_{k = 1}^{N - 1} \sum_{m = 1}^{M} Δ A_{k} [v_{k, m} \cdot v_{k, m}^{T}]

\begin{matrix} K_{q,TH} = - \sum_{k = 1}^{N - 1} \sum_{m = 1}^{M} Δ A_{k} [\frac{1}{Δ r_{k}^{2}} u_{k, m} \cdot u_{k, m}^{T} + \frac{4 q^{2}}{{(r_{k} + r_{k + 1})}^{2}} v_{k, m} \cdot v_{k, m}^{T} \\ + \frac{1}{Δ z^{2}} (v_{k, m + 1} \cdot v_{k, m + 1}^{T} - 2 v_{k, m} \cdot v_{k, m}^{T} + v_{k, m - 1} \cdot v_{k, m - 1}^{T})] \end{matrix}

where

Δ A_{k}

is the area of the annulus with inner and outer radii of

r_{k}

and

r_{k + 1}

,

v_{k, m}

are column vectors that have two non-zero entries each with value

1 / 2

at locations

m \cdot N + k

and

m \cdot N + k + 1

, and

u_{k, m}

are column vectors that have two non-zero entries with values −1 and 1 at locations

m \cdot N + k

and

m \cdot N + k + 1

. The choice to order the vectors

T_{q, l}

first in

k

is arbitrary and leads to contiguous blocks non-zero entries of dimension

N

in the matrices.

To keep the time step update centered in time, the harmonic components of the heat load given by Eq. (27)

Q_{q, l z m} = \frac{1}{2 π} \int_{0}^{2 π} Q (r_{k}, φ, z_{m}, t_{l}) \exp [- i q φ] d φ

must be evaluated at the beginning and end of each time step. Therefore the updated temperature is calculated first using the initial heat distribution. This temperature distribution is then used to update the optical field leading to the updated heat distribution. This is then used to repeat the temperature update, this time incorporating the heat distributions at the beginning and end of the time step. To start the simulation an initial temperature distribution, input mode field, and input pump power must be specified.

The input mode field is taken to be a superposition of guided modes with pre-specified amplitudes. Additionally, a time-dependent phase shift may be applied to one or more of the input modes to incorporate a frequency shift between modes as shown in Eq. (28)

E (r, φ, 0, t) = \sum_{n = 1}^{N} \sqrt{P_{n}} E_{n} (r, φ) \exp [i f_{n} (t)]

where

P_{n}

and

E_{n} (r, φ)

are the power and mode field distribution of guided mode

n

. This has been shown to induce inter-modal coupling through the process of stimulated thermal Rayleigh scattering [9–12].

4. Computational implementation

The model described above is a 3 + 1D fiber amplifier model with a hybrid discretization scheme in cylindrical coordinates employing finite elements in the radial direction, a harmonic expansion in the azimuthal direction, and a finite difference grid in the direction of optical propagation. For a fiber length of 1.63 meters with a beat length of 22.3 millimeters, 50 samples per beat length lead to a longitudinal grid with 3650 points. The finite element approach in the radial dimension permits variable radial point spacing thus reducing the number of radial points to about 100 that span the entire fiber cladding radius, resulting in a matrix order of approximately 3.7 × 10⁵ for each of the $(2 Q + 1)$ thermal matrices. Since the region involved in optical propagation is only at the center of the fiber, a subset of the thermal radial points is used in the optical propagation equations. This region is typically 2~3 times the core size. Figure 1 depicts the computational steps required to implement this method. To reduce the number of time steps required for the amplifier to reach thermal equilibrium, the initial temperature distribution is set to the equilibrium for the case of seeding with the fundamental mode only.

Fig. 1 Sequence of steps required to carry out the time dependent amplifier simulation.

Download Full Size | PDF

This approach can be used for any waveguide profile for which the harmonic expansion can be calculated, however, those with stronger azimuthal variations require higher truncation orders. It is convenient to use separate truncation orders $Q$ for the optical and thermal problems. For modeling photonic crystal fibers with small capillaries truncation orders of $Q_{opt} \sim 10 - 20$ are appropriate while for step index fibers, $Q_{opt}$ can be chosen to include the highest azimuthal order propagating modes supported by the fiber. For typical large mode area fibers $Q_{opt} \sim 2 - 4$ suffices. Likewise, the averaging effect resulting from heat diffusion means that $Q_{therm} \sim 2 - 4$ for all types of fibers. Nevertheless, the orders may be increased until the desired level of convergence is achieved. The simulation time is most sensitive to $Q_{opt}$ because this affects the speed of the beam propagation portion of the calculation which is the most time consuming. The simulation time of a large pitch photonic crystal fiber was observed to be about five times greater than that of a step index fiber.

Within this approach, the primary computational tasks are sparse matrix multiplication and sparse linear system solution. Various software packages are available for these tasks that take advantage of modern multi-processor high performance computing architectures. The implementation reported here relies on parallelization using the message passing interface (MPI) [15]. Within MPI, each collective task on the processor grid is mediated through an MPI communicator. To solve all harmonic components of the thermal problem simultaneously, the processor grid is divided into sub-communicators, one for each harmonic component. Therefore, as long as processors are available, the thermal harmonic order may be increased with negligible increase in the time required to accomplish the thermal update.

During the beam propagation portion of the calculation, the temperature distribution remains stationary in memory scattered across the processing grid and the required temperature values are broadcast to all processes at each spatial propagation step. The hybrid discretization scheme greatly reduces the number of floating point operations required for each update step compared to other discretization techniques. Therefore the optical solution is obtained faster by the group of processors on each sub-communicator performing its own optical update than by spreading the optical update matrix operations across the entire processing grid due to the communications overhead required. As the thermal order is increased, additional communications time is required to retrieve the additional temperature information for the optical update step leading to a gradual increase in the total time to solution.

5. Example results

As a first example, it makes sense to compare the results of this model to a prior model that employed coupled mode theory for the optical propagation [7]. The double-clad Ytterbium-doped amplifier parameters are given in Table 1. Due to the calculation speedup, a step index fiber approximately equivalent to the photonic crystal fiber discussed before was simulated.

Table 1. Simulated Fiber Parameters

View Table

Furthermore, there was no material, scattering, or bend loss assumed. The outer boundary of the fiber cladding was maintained at a fixed reference temperature thus assuming perfect conductive cooling. Also, the pump and signal linewidths were assumed negligible. A counter-pumped configuration was assumed for this simulation in which the total pump absorption throughout the fiber was approximated to be constant so that the pump intensity at the seeded end could be set as a boundary condition.

Figure 2 and the accompanying movie show the time dependence of the optical intensity within the fiber core and the output intensity distribution over time. The simulation was started with 9.5 Watts of seed power launched into the LP₀₁ mode and 0.5 Watts in the LP_11e mode and no launched pump power. The pump was then linearly increased over the first 10ms of the simulation. In agreement with the coupled mode theory results [7] the instability was confined to the latter portion of the amplifier and was more severe near the output end. Furthermore, as time progressed the stable region grew in extent as the thermally induced grating reached stable equilibrium as shown in Fig. 2. This equilibrium region was possible due to the fact that the two modes were launched at the same optical frequency and so no moving grating existed at the seed end [6]. The presence of an equilibrium region is also in agreement with the coupled mode theory results [7].

Fig. 2 Plot of the optical intensity as a function of length at the center of the fiber core and at points offset from the center of the fiber core by half the core radius after a simulation time of 14 ms (a) and 20 ms (b) for an output power of 287 Watts. The accompanying movie (Media 1) shows the evolution of these intensity probes over a time period of 40 ms. The insets shows the optical field at the output. The markers in the insets show the locations of the probes.

Download Full Size | PDF

While there are some similarities between the appearances of Fig. 2 and Fig. 2 of [7] different quantities are being plotted. Figure 2 shows the local intensity probe at three points while the prior published figure shows the modal content. The full optical field that is required to calculate the modal content was not stored at every longitudinal position for every time step due to the large file size that would be required. Nevertheless, rapid spatial oscillations characteristic of non-adiabatic power transfer [8] were evident in both. The beam propagation model described here also captures the effect of thermal lensing on the mode field area of the fundamental mode of the fiber. The mode field area of the cold fiber was 2750 µm² and at the output end of the heated fiber, the area decreased to 2100 µm² at an output power of 287 Watts.

Figure 3 shows the time evolution of the output modal decomposition. The modal content was calculated with respect to the cold fiber modes as shown in Fig. 4. The overall mode coupling behavior exhibited three distinct regimes. The first was the onset regime where the pump power was still being ramped up from 0 to 10 ms. Although significant coupling into the LP_11e mode immediately begins to occur in this regime as the pump level is increased, the time variation is still relatively slow. The second regime was characterized by very rapid mode coupling behavior that included additional modes and did not exhibit any periodicity.

Fig. 3 Modal content of the amplifier output over time as well as total power from the integrated optical output intensity and the sum of the first 8 modal powers for amplifier output powers of 241 (a), 270 (b), and 287 (c) Watts. Plot (d) is a zoomed in view of (c) to show the details of the rapid fluctuations between 10 and 22 ms. The modal content was calculated with respect to the cold fiber modes. The pump was ramped on linearly in time over the first 10 ms in each case.

Download Full Size | PDF

Fig. 4 Plots of the intensity profiles of the first 8 guided modes of the cold fiber incorporated in the simulated amplifier: LP₀₁(a), LP_11o(b), LP_11e(c), LP_21o(d), LP_21e(e), LP₀₂(f), LP_31o(g), LP_31e(h).

Download Full Size | PDF

All modes of the same azimuthal and radial order are perfectly degenerate in this implementation; therefore the symmetry axes of the different pairs of mode are oriented differently so the symmetry axis of the LP_21e is not the same as that of the LP_11e mode as seen in Fig. 4. For times near the end of the simulation the power was primarily in the two launched modes although the LP₀₂ content increased due to thermal lensing. The onset of this regime coincided with the stable region of the fiber reaching its maximum extent as can be seen from the intensity profiles shown in Media 1. The gradual decrease in the severity of the mode fluctuations in time may be attributable to the decrease in effective mode area caused by increased thermal lensing as the fiber continues to heat up. As the pump power was increased, the modal fluctuations increased in depth.

Since the thermal boundary condition was symmetric in this simulation, it is not surprising that the output intensity distribution mostly retained mirror symmetry about the axis of the launched LP_11e mode. The observed small temporary asymmetries were attributed to the fact that the sparse linear solvers used to accomplish the beam propagation were of the iterative variety so that each propagation step and thermal update introduced some error into the optical fields and temperature profile depending on the chosen convergence tolerance. Simulations performed with larger tolerances exhibited more severe asymmetry and the asymmetry became more noticeable at higher powers.

The total output power was found to fluctuate slightly. These slight fluctuations may be consistent with reported observations of constant total output power within experimental precision [4]. This is attributable to the fact that the spatial fluctuations of intensity within the fiber inhibit efficient extraction of the available gain thus causing the uniform pump absorption approximation to break down. This observation suggests an additional mechanism for feedback that can drive the instability. It has been shown that amplitude fluctuations at the seeded end of amplifiers can cause a steady-state amplifier solution to become unstable [9]. If fluctuations in the spatial intensity cause varying amounts of pump to be absorbed in a counter-pumped configuration, this would effectively modulate the seed level present near the input end of the amplifier thus providing the feedback mechanism modulating the amplitude causing instability. Also, it is important to note that the total power contained in the first 8 modes was very close to the total optical power at the output indicating that this is an appropriate number of modes to use to characterize the output of this particular amplifier.

It is also apparent that the incorporation of the beam propagation model significantly lowers the instability threshold. The maximum pump power of 340 Watts in the simulations presented here was below the previously observed approximate instability threshold of 1060 Watts [7] for the conductively-cooled case and yet instability is clearly present. This suggests that models of the transient regime based on coupled mode theory should include all guided modes in order to accurately capture instability behavior. This also agrees with the observation that larger cores that support more transverse modes exhibit lower instability threshold powers. Experimentally-observed instability thresholds for this type of amplifier are yet lower so clearly the incorporation of the beam propagation method improves prospects for agreement of the model with experimentally observed instability thresholds.

Figure 5 shows the frequency spectra of the transverse modes at the fiber output for output powers of 241, 270, and 287 Watts. No sharp frequency peaks are visible indicating a degree of randomness in the output. Also, the high frequency tail was more prominent in the lower-order modes, which is not immediately evident just by examining the time series shown in Fig. 3. The temperature was recorded in time at the same sampling points as the optical intensity throughout the fiber and is presented in Fig. 6. As expected, the temperature is static in the region of the fiber not experiencing instability but is irregular where the optical intensity fluctuates. Due to the length of the thermal diffusion time, the temperature cannot keep up with the rapidly fluctuating intensity. Rather, it exhibits slowly varying behavior that captures a windowed time average of the quantum defect heating profile. Both the center temperature and the amplitude of the temperature fluctuations are irregular and not uniform across the core as revealed by the difference between the profiles at the core center and off center. Observing the time evolution of the temperature profile from the beginning of the simulation shown in Media 2 reveals that the temperature grating strengthens while the interference pattern is stationary. As soon as the interference pattern begins to move rapidly, the intensity of the grating decreases as the minima and maxima are averaged out.

Fig. 5 Frequency spectra of the first 8 modes of the simulated fiber calculated over the simulation period of 100 ms for output powers of 241 (a), 270 (b), and 287 (c) Watts.

Download Full Size | PDF

Fig. 6 Temperature distribution as a function of length at the fiber core and at points offset from the center of the fiber core by half the core radius at a time of 14 ms or the amplifier operating at 287 Watts of output power. The accompanying movie (Media 2) shows the evolution of these temperature probes over a time period of 40 ms.

Download Full Size | PDF

6. Model limitations, and future investigations

The model presented here has some limitations in its current form that deserve discussion. The first is that bending losses and bend-induced mode distortion are not accounted for. An effective index gradient may be incorporated in the usual fashion to account for bending [16]. This breaks the azimuthal symmetry of the step-index profile leading to additional terms in the harmonic azimuthal expansion of the index. The optical boundary condition use here is lossless and perfectly reflecting. It is possible that transparent boundary conditions [17] could be adapted to the paraxial scalar beam propagation presented here as was done for a previous similar method [13]. Furthermore, material loss is not accounted for. In the current implementation, the doped region of the core must be assumed to be uniform and circular. While this is true of some step index large mode area fibers, micro-structured fibers often have some additional structure that arises from the stack and draw fiber fabrication method. As was already mentioned above, the presented model does not hold the launched pump perfectly constant for counter-pumped configurations which are practically a necessity for minimizing the effective amplifier length to maximize non-linear thresholds. An efficient method for overcoming this limitation in the transient regime is not immediately apparent. Finally, in its current form, this model requires high performance computing resources to use. It is not suitable for desktop computers.

The model presented here can enable numerous future investigations a few of which are briefly discussed below. First, it is important to verify that in the absence of noise, quantum or otherwise, the amplifier can reach stable thermo-optic equilibrium. While a previous investigation predicted dynamic instability even in the absence of such noise [7], this could have been due to limits in the achievable accuracy of the solutions of the derived equations. This instability was also predicted to occur at power levels much higher than those at which has been observed suggesting that noise plays an important role in the origin of instability.

This model could also be used to analyze the effects of different types of noise on the onset of instability. This could include frequency offsets between modes, fluctuations of pump and seed powers, and other time varying launch conditions. Build up and decay times of instability could also be studied. The potential for increasing the instability threshold through advanced fiber designs could also be studied. For example the model can treat large pitch photonic crystal fibers that rely on de-localization of higher order modes for fundamental mode discrimination.

7. Conclusion

A model of transient modal instability in fiber amplifiers has been presented. This model relies on a time-dependent 3-dimensional treatment of the effects of quantum defect heating on the waveguiding properties of amplifier fibers. This model has confirmed some predictions first made using coupled-mode theory. These include the existence of stable and unstable regions along the length of the fiber, their evolution over time, the increase in the severity of the instability toward the output end of the fiber, and rapid optical intensity variations along the fiber. It has also exhibited several additional aspects of TMI including coupling to un-seeded modes and lower onset threshold powers compared to a prior model that are more in line with experimental observations.

Efforts have been made to reduce the time required to perform the calculations including the realization of a new hybrid finite element, harmonic beam propagation method, a hybrid finite element, harmonic, finite difference thermal solver, and parallel implementation on modern high performance computing architectures using the message passing interface. This method should prove a useful tool in studying and eventually overcoming the technological challenges presented by TMI in high average power fiber amplifiers.

Acknowledgments

This work was made possible by a grant of computing time at the U S. Army Engineer Research and Development Center Department of Defense Supercomputing Resource Center from the U. S. Department of Defense High Performance Computing Modernization Program. This work was also supported by the US Department of Defense High Energy Laser Joint Technology Office under the Multidisciplinary Research Initiative “Fiber Laser Light Engines.” The author would like to thank R. Andrew Motes for helpful discussions. The views expressed in this article are those of the author and do not reflect the official policy or position of the US government or the Department of Defense. Distribution A: Approved for public release, distribution unlimited.

References and links

1. T. Eidam, C. Wirth, C. Jauregui, F. Stutzki, F. Jansen, H. J. Otto, O. Schmidt, T. Schreiber, J. Limpert, and A. Tünnermann, “Experimental observations of the threshold-like onset of mode instabilities in high power fiber amplifiers,” Opt. Express 19(14), 13218–13224 (2011), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-14-13218. [CrossRef] [PubMed]

2. C. Jauregui, T. Eidam, J. Limpert, and A. Tünnermann, “The impact of modal interference on the beam quality of high-power fiber amplifiers,” Opt. Express 19(4), 3258–3271 (2011), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-4-3258. [CrossRef] [PubMed]

3. H. J. Otto, F. Stutzki, F. Jansen, T. Eidam, C. Jauregui, J. Limpert, and A. Tünnermann, “Temporal dynamics of mode instabilities in high-power fiber lasers and amplifiers,” Opt. Express 20(14), 15710–15722 (2012), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-14-15710. [CrossRef] [PubMed]

4. M. Karow, H. Tünnermann, J. Neumann, D. Kracht, and P. Weßels, “Beam quality degradation of a single-frequency Yb-doped photonic crystal fiber amplifier with low mode instability threshold power,” Opt. Lett. 37(20), 4242–4244 (2012), http://www.opticsinfobase.org/ol/abstract.cfm?URI=ol-37-20-4242. [CrossRef] [PubMed]

5. N. Haarlammert, O. de Vries, A. Liem, A. Kliner, T. Peschel, T. Schreiber, R. Eberhardt, and A. Tünnermann, “Build up and decay of mode instability in a high power fiber amplifier,” Opt. Express 20(12), 13274–13283 (2012), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-12-13274. [CrossRef] [PubMed]

6. A. V. Smith and J. J. Smith, “Mode instability in high power fiber amplifiers,” Opt. Express 19(11), 10180–10192 (2011), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-11-10180. [CrossRef] [PubMed]

7. B. Ward, C. Robin, and I. Dajani, “Origin of thermal modal instabilities in large mode area fiber amplifiers,” Opt. Express 20(10), 11407–11422 (2012), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-10-11407. [CrossRef] [PubMed]

8. C. Jauregui, T. Eidam, H. J. Otto, F. Stutzki, F. Jansen, J. Limpert, and A. Tünnermann, “Physical origin of mode instabilities in high-power fiber laser systems,” Opt. Express 20(12), 12912–12925 (2012), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-12-12912. [CrossRef] [PubMed]

9. K. R. Hansen, T. T. Alkeskjold, J. Broeng, and J. Lægsgaard, “Theoretical analysis of mode instability in high-power fiber amplifiers,” Opt. Express 21(2), 1944–1971 (2013), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-21-2-1944. [CrossRef] [PubMed]

10. L. Dong, “Stimulated thermal Rayleigh scattering in optical fibers,” Opt. Express 21(3), 2642–2656 (2013), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-21-3-2642. [CrossRef] [PubMed]

11. A. V. Smith and J. J. Smith, “Influence of pump and seed modulation on the mode instability thresholds of fiber amplifiers,” Opt. Express 20(22), 24545–24558 (2012), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-20-22-24545. [CrossRef] [PubMed]

12. A. V. Smith and J. J. Smith, “Steady-periodic method for modeling mode instability in fiber amplifiers,” Opt. Express 21(3), 2606–2623 (2013), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-21-3-2606. [CrossRef] [PubMed]

13. S. A. Shakir, R. A. Motes, and R. W. Berdine, “Efficient scalar beam propagation method,” IEEE J. Quantum Electron. 47(4), 486–491 (2011), http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=5730173. [CrossRef]

14. P. D. Dragic and B. G. Ward, “Accurate modeling of the intrinsic Brillouin linewidth via finite-element analysis,” IEEE Photon. Technol. Lett. 22(22), 1698–1700 (2010), http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=5590281. [CrossRef]

15. S. Balay, J. Brown, K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang, PETSc Users Manual, (ANL-95/11 - Revision 3.3, Argonne National Laboratory, 2012), http://www.mcs.anl.gov/petsc/.

16. R. T. Schermer and J. H. Cole, “Improved bend loss formula verified for optical fiber by simulation and experiment,” IEEE J. Quantum Electron. 43(10), 899–909 (2007), http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4300920&isnumber=4294077. [CrossRef]

17. G. R. Hadley, “Transparent boundary condition for the beam propagation method,” IEEE J. Quantum Electron. 28, 0.363–370 (1992), http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=119536&isnumber=3419. [CrossRef]

Parameter	Value
core diameter	74 µm
pump cladding diameter	170 µm
outer cladding diameter	400 µm
core numerical aperture	0.03
fundamental mode field area	2750 µm²
beat length LP₀₁-LP₁₁	22.3 mm
fiber length	1.63 m
yb⁺³ doping concentration	3.5 × 10²⁵ m⁻³
signal wavelength	1.064 µm
pump wavelength	0.977 µm
signal power LP₀₁	9.5 W
signal power LP₁₁	0.5 W
pump power	286,321,340 W
output power	241,270,287 W
signal emission cross-section	3.58 × 10⁻²⁵ m²
signal absorption cross-section	6.00 × 10⁻²⁷ m²
pump emission cross-section	1.87 × 10⁻²⁴ m²
pump absorption cross-section	1.53 × 10⁻²⁴ m²
upper state lifetime	850 µs
thermal conductivity	1.38 W/m-K
heat capacity	703 J/kg-K
mass density	2200 kg/m³
thermo-optic coefficient	1.2 × 10⁻⁵ K⁻¹
heat sink temperature	300 K
z grid spacing	4.46 × 10⁻⁴ m
time step	4.0 × 10⁻⁵ s
total simulation time	100 ms
thermal radial points	101
optical radial points	56
doped core radial points	30
optical azimuthal order	4
thermal azimuthal order	3

Modeling of transient modal instability in fiber amplifiers

Abstract

1. Introduction

2. Hybrid finite element-harmonic beam propagation method

3. Gain and heat deposition

4. Computational implementation

5. Example results

6. Model limitations, and future investigations

7. Conclusion

Acknowledgments

References and links

Supplementary Material (2)

Cited By

Figures (6)

Tables (1)

Equations (28)

Optics Express