## Abstract

A model of transient modal instability in fiber amplifiers is presented. This model combines an optical beam propagation method that incorporates laser gain through local solution of the rate equations and refractive index perturbations caused by the thermo-optic effect with a time-dependent thermal solver with a quantum defect heating source term. This model predicts modal instability a fiber amplifier operating at 241, 270, and 287 Watts of output power characterized by power coupling to un-seeded modes, the presence of stable and unstable regions within the fiber, and rapid intensity variations along the fiber. The instability becomes more severe as the power is increased.

©2013 Optical Society of America

## 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

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)

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

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)

## 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^{5} for each of the $\left(2Q+1\right)$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.

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}_{\text{opt}}\sim 10-20$are appropriate while for step index fibers, ${Q}_{\text{opt}}$can be chosen to include the highest azimuthal order propagating modes supported by the fiber. For typical large mode area fibers ${Q}_{\text{opt}}\sim 2-4$suffices. Likewise, the averaging effect resulting from heat diffusion means that ${Q}_{\text{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}_{\text{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.

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_{01} 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].

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^{2} and at the output end of the heated fiber, the area decreased to 2100 µm^{2} 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.

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_{02} 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.

## 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]