Accurate and Efficient Modeling of the Transverse Mode Instability in High Energy Laser Amplifiers

We study the transverse mode instability (TMI) in the limit where a single higher-order mode (HOM) is present. We demonstrate that when the beat length between the fundamental mode and the HOM is small compared to the length scales on which the pump amplitude and the optical mode amplitudes vary, TMI is a three-wave mixing process in which the two optical modes beat with the phase-matched component of the index of refraction that is induced by the thermal grating. This limit is the usual limit in applications, and in this limit TMI is identified as a stimulated thermal Rayleigh scattering (STRS) process. We demonstrate that a phase-matched model that is based on the three-wave mixing equations can have a large computational advantage over current coupled mode methods that must use longitudinal step sizes that are small compared to the beat length.


Introduction
Ytterbium-doped fiber amplifiers that produce kilowatt output powers have been developed in the past decade [1][2][3][4][5][6][7]. However, the thermal or transverse mode instability (TMI) has become a major barrier to achieving even higher output powers [8][9][10]. Despite almost a decade of work since the original observation of TMI in fiber amplifiers by Jauregui et al. [11] and Eidam et al. [12], it still remains only partially understood, and computationally-efficient methods that are sufficiently accurate for amplifier design have been lacking. It was recognized early by Smith and Smith [13] that the instability could be explained by a thermal grating that is induced by the beating of the fundamental mode of the optical fiber with a higher-order mode at a slightly lower frequency and the quantum defect heating that ensues.
Subsequent work by Jauregui et al. [14], Dong [15], and Smith and Smith [16] identified the instability as a stimulated thermal Rayleigh scattering (STRS) process. In particular, Dong [15] developed a three-wave mixing model that is analogous to models of the Brillouin instability due to stimulated Brillouin scattering (SBS). This identification has remained somewhat controversial, although Kong et al. [17] directly observed the STRS process in a fiber amplifier. Ward et al. [18] and Naderi et al. [19] developed a model of TMI based on a coupled mode method that makes no reference to three-wave interactions. The complexity of TMI has obscured its identification as an STRS-driven, three-wave process. The conditions that are required to treat TMI as a three-wave instability have not been elucidated.
In this work, we demonstrate that the key requirement is that the beat length B = 2 /Δ must be small compared to any other longitudinal scale lengths, where Δ is the difference between the wavenumber of the fundamental mode and any higher-order modes (HOMs) at the same frequency. This condition usually applies in practice. In this limit, we derive the three-wave equations that govern TMI. Terms that are not phase-matched are neglected. This approach is similar to Dong's approach [15], but is more general.
We further demonstrate that this approach has a large computational advantage. Approaches that use the coupled mode method, like the approach of Naderi et al. [19], must use longitudinal steps in their computations that are small compared to the beat length. As a result, Naderi et al. [19] limited their study to a high-gain amplifier with a length of 1.6 m, which is substantially shorter than typical ytterbium-doped fiber amplifiers. Approaches that use the finite-element method to calculate the optical mode profile on each longitudinal step like that of Ward [20] must take steps that are small compared to the optical wavelength and typically require large computational resources. Our approach is only limited by the longitudinal scale lengths over which the amplitudes of the optical modes, the thermal mode, and the pump mode change. These lengths are typically far larger than the beat length. For the examples that we consider in this paper, which correspond to realistic Yb 3+ -doped fiber amplifiers, the computational speedup is more than a factor of 100.

Theoretical Model
We may write the index of refraction as (r ⊥ , , ) = 0 (r ⊥ ) + Δ (r ⊥ , , ), where 0 (r ⊥ ) is the unperturbed index of refraction, and we set Δ 0 , which is always the case (Δ / 0 ∼ 10 −5 ). We will use the slowly varying envelope approximation, which is an excellent approximation in this system due to the large discrepancy between the wavelength and the next-smallest scale, which is the intermodal beat length B . We will also assume that time derivatives of the index of refraction can be ignored when calculating the inter-modal coupling. That is again an excellent approximation, given the large disparity between the light period and the time scale on which either the gain changes (microseconds) or the temperature changes (milliseconds). We assume as well that the only coupling is between modes that are propagating in the forward direction. We start with the expression from coupled mode theory for two coupled modes [21,22], where E 0 and E 1 are the normalized transverse mode profiles for the fundamental mode and the HOM, while 0 and 1 are the corresponding amplitudes. We have set ≡ 0 ≈ 1 since Δ = 0 − 1 0,1 . Equation (1) is valid when only a single HOM is present or when the amplitudes of the other HOMs are small compared to 0 and 1 .
To solve Eq. (1) at any point in time , we must find Δ (r ⊥ , , ). Due to the factors exp(± Δ ) that appear in Eq. (1), it is necessary to determine Δ (r ⊥ , , ) with a computational resolution Δ that is small compared to the beat length B = 2 /Δ even though the field amplitudes vary slowly compared to this length. Since we must determine the transverse dependence of Δ as well, this constraint is computationally demanding. We can bypass this difficulty by replacing the field Δ (r ⊥ , , ) with three fields 0 (r ⊥ , , ), + (r ⊥ , , ), and − (r ⊥ , , ), which are defined so that and all three fields vary slowly compared to the beat length B . When we substitute Eq. (2) into Eq. (1) and keep only the phase-matched terms, we obtain The terms that are not phase-matched are rapidly oscillating and do not contribute significantly to the integrals. All terms in Eq. (3) vary slowly in . It then becomes possible to integrate Eq. (3) with no loss of accuracy while using a resolution in that is far larger than is necessary with Eq. (1).
The form of Eq. (3) makes clear that in the limit where Eq. (3) is valid, TMI is effectively a three-wave process in which two optical fields combine with a density field. The condition for Eq. (3) to be valid is that all terms must vary slowly compared to the beat length. This condition is almost always met in practice. Since the density field is thermally driven, this process is a stimulated Rayleigh scattering process. We discuss this identification in more detail in the Appendix.
The procedure that we use to obtain Δ 0 , Δ + , and Δ − in terms of the field and as a function of time is analogous to the procedure that is used by Naderi et al. [19] that is illustrated in Fig. 1.
We begin by writing the signal intensity (r ⊥ , , ) as all vary slowly compared to . The behavior of a Yb-doped fiber amplifier is accurately described as a two-level system at realistic power levels [2]. Given (r ⊥ , , ) and the pump intensity ( , ), we next compute the upper state density 2 (r ⊥ , , ) of the Yb ions, which is given by 2 = where and are the signal and pump angular frequencies, respectively, is the spontaneous decay time of the upper level, 0 is total density of Yb ions, and (a) , (e) , (a) , and (e) are the absorption and emission cross-sections at the pump and signal frequencies. Because appears in both the numerator and denominator of Eq. (6), the density 2 will have higher harmonics that are proportional to exp(± Δ ) with > 1. We will truncate this expression since the harmonics with > 1 are not phase-matched. Explicitly, we find that Eq. (6) has the form where = Δ + and = ∠( We note that < and < is a consequence of + < 0 , which in turn follows from the Cauchy-Schwartz inequality. We now write where 20 = with = (1 − 2 / 2 ) 1/2 [23]. Equation (10) is an exact truncation, not an expansion. We will show in Sec. 3 that the contribution of higher harmonics with > 1 are negligible. The next stages in the procedure are more straightforward. TMI is generated by the heat deposition due to the quantum defect between the pump and the signal, which in turn leads to a time-delayed temperature response that changes the index of refraction. The temperature response depends linearly on the heat deposition, which in turn depends linearly on the upper state density. From the expression for the heat deposition , we find

20
, Similarly, from the expression for the temperature evolution, where is the density, is the heat capacity, and is the heat diffusivity, we find  , ). Since the temperature tends to a constant room at large radius, the appropriate boundary conditions for both and 0 at large radius are 0 = room , and the appropriate boundary condition for + is + = 0. This integration is where the basic time step occurs, as we show schematically in Fig. 1, and it is this step that is computationally time-consuming.
We can now find the change in the index of refraction. There are two contributions to the index of refraction that we must take into account. The first contribution is from the temperature change, for which Δ = ( / ) ( − room ) and Δ 0 = ( / ) ( 0 − room ), Δ + = ( / ) + in the phase-matched model. The second contribution is from the gain, from which we find It follows that where Although Δ is purely imaginary, the components Δ + and Δ − are not. Naderi et al. [19] have pointed out that the corresponding phase shift does not contribute to the instability, but plays an important role in the energy balance. Finally, we obtain Δ 0 = Δ 0 + Δ 0 , Δ + = Δ + + Δ + , and To complete the model equations, we must obtain the pump intensity ( , ). We use the expression where the overbar indicates that the population densities are averaged over the cross-section of the pump. The sign depends on whether the pump is forward-or backward-propagating. Eq. (18) becomes in the phase-matched model.

Verification, Accuracy, and Timing of the Phase-Matched Model
In this section, we first verify the phase-matched model, [Eqs. (3), (8), (12), (14), (17), (19)] by comparing its predictions to those of the full model, [Eqs. (1), (6), (11), (13), (16), (18)]. We will show that agreement is excellent for a realistic amplifier system similar to the system that Naderi et al. [19] considered, but using a fiber length of 10 m, which is a typical experimental length [7,24]. We then consider in more detail the error as a function of the step size Δ and show that the phase-matched model has a significant computational advantage.

Verification
We show the basic set of parameters that we are considering in Table 1. These parameters are similar to those used in [19], but with a more realistic amplifier length of 10 m. We use the alternating-direct-implicit (ADI) method to integrate the temperature equations, and we used the Runge-Kutta method to carry out the -integration. In all the simulations reported here, we used a 140 × 140 m 2 square grid for the transverse integration with a 2 × 2 m 2 grid spacing. We verified that using a larger grid size of 160 × 160 m 2 with a smaller grid spacing of 1 × 1 m 2 makes a negligible difference in Figs. 2 and 3. We chose the z-step sufficiently small so that the relative error is below 1%. In all the simulations reported here, we used a noise seed ratio at the entry to the optical fiber of 10 −4 between the higher-order mode and the fundamental mode. We verified that using a noise seed ratio of 10 −3 or 10 −5 does not change the agreement between the full and phase-matched models that we present in Figs. 2 and 3. In Fig. 2, we show a comparison of the ratio ( ) = HOM ( )/ total ( ) of the power in the higher-order mode HOM ( ) to the total power total ( ) at the end of the fiber as a function of time. In the case that we show here, the input pump power pump equals 250 W. With this pump power, we find that ( ) rises to a maximum of 17%, shown as a dot in Fig. 2, before returning to a value that is close to the initial higher-order mode seeding. We observe excellent agreement between the full model and the phase-matched model. As the pump power increases, we continue to observe excellent agreement between the two models, although when the pump power is large enough where both models predict chaotic oscillations, the agreement is qualitative rather than quantitative. This behavior is expected since small changes in the seeding ratio also produce large changes in this limit due to the butterfly effect [25].
In Fig. 3, we show max[ ( )] vs. pump and a dotted line that corresponds to a ratio of 1%. We observe excellent agreement between the phase-matched model and the full model. In this work, we define the threshold power as the lowest pump power at which max[ ( )] > 0.01, i.e., the ratio of HOM ( )/ total ( ) exceeds 1% at any time. Beyond this threshold, the beam quality rapidly degrades [11,12,26,27]. This definition of the threshold is consistent with studies of amplifier limits due to SBS. In the case considered here, the threshold power equals 207 W.
In Fig. 4, we show the temperature as a function of longitudinal position and the absolute value of its spatial Fourier transform at a distance of 10 m from the amplifier center and a time = 0.5 ms. We set pump = 450 W. In Fig. 4(a), where we show the temperature as a function of position, we see that the agreement between the full model and the phase-matched model appears excellent. In the inset, where we show the spatial oscillations, the two models appear indistinguishable. However, the subtle difference is visible in Fig. 4   |FT( )| vs. , where is the spatial Fourier transform variable. Agreement for the central harmonic and the harmonics at = ±Δ is excellent. The inset shows excellent agreement for the harmonic at = Δ . However, the harmonics at = ± Δ for ≥ 2 are not present in the phase-matched model. excellent for the central harmonic, as well as the two surrounding harmonics which are located at = ±Δ = ±528 m −1 . However, the phase-matched model has no contribution from the harmonics at ± Δ , where ≥ 2. It is precisely these higher harmonics that we are neglecting.

Accuracy and Timing
In the phase-matched model, the number of dependent variables is almost twice as large as in the full model. In particular, it is necessary to solve the temperature equation, Eq. (14), for both 0 and + instead of just solving the temperature equation, Eq. (13), for the single temperature . As a result, we have found that the computational load on each -step increases by roughly a factor of two. However, it is possible to take significantly larger steps, leading to a large computational advantage. In Fig. 5, we show max[ ( )] for both the full model and phase-matched model as a function of B /Δ when pump = 250 W, so that the pump power is slightly above threshold. As expected, the full model requires an B /Δ > 60 to converge, while the phase-matched model appears to have converged in this case when B /Δ 2. In order to quantify the convergence, we define the relative error, , as the difference between our computation at a given Δ and a four-point Richardson extrapolation [28]. For the full model, we used B /Δ = 80, 40, 20, and 10 for the extrapolation. For the phase-matched model, we used B /Δ = 40, 20, 10, and 5.
In Fig. 6(a), we show the relative error as a function of B /Δ . We see that achieving a relative error of 1% with the full model requires B /Δ 90, while the same relative error can be obtained with the phase-matched model when B /Δ 2. Figure 6(a) also illustrates that the full model is second-order accurate in Δ , so that the relative error decreases proportional to (Δ ) −2 . This result is consistent with the result of Naderi et al. [19], but may be surprising since our integration in is done using the Runge-Kutta method. This result indicates that the global error is dominated by the accumulated error in computing the index of refraction. The variation of the relative error in the phase-matched model is more complex since it depends on the rate at which all the dependent variables change as a function of . A complete error analysis is beyond the scope of this paper, but Fig. 6(a) indicates that it decreases rapidly as Δ increases until it has become quite small.
In Fig. 6(b), we show the runtime of both the full model and the phase-matched model using 16 cores in a shared memory system that consists of dual Intel E5-2695 V4 processors. We observe that the runtime for the phase-matched model is approximately twice as the runtime for the full model, which is consistent with the greater computational load per step in the phase-matched model. When we compare the runtime corresponding to a relative error of 1%, shown with dots, we found a runtime of 163 hours for the full model and 1.17 hours for the phase-matched model, indicating that the phase-matched model runs a factor of 139 faster in this case. While we obtain a speedup of 139 in our study, this number will vary with different simulation parameters. Nevertheless, the advantage of using phase-matched model is clear.

Conclusions and Discussion
In this work, we derived the three-wave mixing equations that govern TMI in the limit where a single higher-order mode is present, and the longitudinal rate of change of all quantities is slow compared to the beat length. This limit normally applies in practice. In this limit, TMI can be identified as an STRS process, and we reviewed the theoretical justification for this identification in the appendix, where we also discussed similarities and differences between STRS and SBS in the Appendix. There, we verified the accuracy of the phase-matched model in a Yb 3+ -doped fiber amplifier with a relatively simple step index profile. The amplifier that we considered is like that of Naderi et al. [19], but has a more realistic 10-m length. We demonstrated that this model reproduces the nonlinear saturation of the higher-order mode and the instability threshold that are predicted by the full model. We demonstrated a computational speedup that is more than a factor of 100.
We derived the three-wave mixing equations in the case that a single higher-order mode is present, but we expect this result to be more broadly applicable when several higher-order modes are present. In the linear limit below threshold, each higher-order mode will only interact with the fundamental mode. In that case, the three-wave mixing equations can be extended by adding a new set of equations for the index of refraction, Yb 3+ population density, heat, temperature, and optical mode amplitude for each of the higher-order modes. The computational complexity scales proportional to , where is the number of modes. More generally, we anticipate that the three-wave mixing equations can be extended to include a coupling between all the modes as long as none of the beat lengths between any of mode pairs becomes large enough to be comparable to the scale length on which any of the amplitudes change. However, the computational complexity grows proportional to 2 , and higher-order nonlinear interactions with a slowly varying amplitude could invalidate this approach.