## Abstract

We present detailed numerical simulations of modal instabilities in high-power Yb-doped fiber amplifiers using a time-dependent temperature solver coupled to the optical fields and population inversion equations. The temperature is computed by solving the heat equation in polar coordinates using a 2D second-order alternating direction implicit method. We show that the higher-order modal content rises dramatically in the vicinity of the threshold and we recover the three power-dependent regions that are characteristic of the transfer of energy. We also investigate the dependence of the threshold on the seed power and the modal content ratio of the seed. The latter has a minimal effect on the threshold while it is shown that for the fiber configuration investigated, the modal instability threshold scales linearly over a wide range with the seed power. In addition, two different gain-tailored core designs are investigated and are shown to have higher thresholds than that of a uniformly doped core. Finally, we show that this full time-dependent model which does not assume a frequency offset between the modes a priori, predicts a reduced threshold when the seed is modulated at the KHz level. This is in agreement with the steady-periodic approach to this phenomenon.

©2013 Optical Society of America

## 1. Introduction

High-power narrow-linewidth laser sources with good beam quality are of considerable interest due to their utility in several applications including coherent beam combination [1] and nonlinear frequency conversion [2]. Diode-pumped rare-earth doped fiber amplifiers are highly desirable as such laser sources as they are blessed with small footprints, superior thermo-optical properties, and high electrical to optical and optical to optical efficiencies. However, fiber amplifiers are limited by the onset of optical nonlinearities including stimulated Brillouin scattering (SBS), stimulated Raman scattering (SRS) and self-phase modulation (SPM). It is well-known that the thresholds of these nonlinearities are generally proportional to the effective area of the fiber core. Consequently, a considerable amount of effort has been devoted over the past decade to increasing the fiber core size while maintaining single-mode operation. Unfortunately, several research groups have reported recently that large mode area (LMA) high-power fiber amplifiers are susceptible to the onset of modal instabilities [3,4]. These instabilities are characterized by a sudden and dramatic increase in the energy within higher-order modes (especially the LP_{11}) as well as rapid temporal oscillations.

There have been several theoretical and experimental investigations of the physics of modal instabilities [3–9]. It is now generally agreed upon that the onset of these instabilities is a result of interplay between thermal effects and modal interference. This interplay was first suggested in the literature by Jauregui et al. [5]. Subsequent work by Smith and Smith [6] emphasized the need to include the temporal dynamics to account for the energy transfer between the fundamental and higher-order modes. In order to achieve this, they introduced a frequency offset between the two modes and a steady-periodic treatment of the heat equation. It was not clear in that work how to accurately determine the frequency offset other than by trial and error. Furthermore, by allowing only one frequency offset, contributions from the other frequency components were neglected. Semi-analytical approaches by Hansen et al. [7,8] and Dong [9] that built on the frequency offset proposition provided valuable insights, but due to the difficulty of finding exact solutions did not resolve the spatial and temporal evolution of laser gain. Alternatively, the work of Ward et al. [10] showed that a full time-dependent treatment of the heat equation can capture much of the physics of the modal instabilities without assuming a frequency offset a priori. Furthermore, this work was able to capture the three regimes of temporal dynamics described in the experimental work of Otto et al. [11].

In this paper, we present extensive theoretical and computational analysis of the temporal dynamics of modal instabilities in fiber amplifiers. Our analytical treatment of laser gain, amplitude evolution equations, and the heat equation provides formal proof of the need for a dynamic approach to the modal instabilities. In this treatment, the thermal load is dictated by population inversion rather than assuming it simply follows the seed intensity profile. The time-dependent heat equation subject to a Dirichlet boundary condition is solved in polar coordinates using the alternating direction implicit (ADI) method. The implementation of this method in polar coordinates is typically avoided due to the occurrence of an artificial singularity at the center of the coordinate system and due to the added complexity of applying the appropriate boundary conditions. However, since the polar coordinate system is more suitable to describe the fiber geometry, we have proceeded with implementation of this technique. The propagation of the electromagnetic waves is described using a coupled-mode theory approach. The resultant system of coupled nonlinear equations is solved using a fourth-order Runge-Kutta solver. Low resolution of the spatial and temporal numerical grid can potentially lead to inaccurate solutions including an artificially induced transfer of energy. Here, every numerical evaluation is verified and quantum efficiency as related to energy conservation has been monitored to track any possible numerical errors for an accurate computation of the field propagation. Numerical simulations exploring dependence on pump power, seed power, and seed modal ratio are presented. In addition, two different gain-tailored fiber designs are simulated and compared. It is shown numerically that both gain-tailoring and increased seed power can mitigate the modal instabilities. Finally, the effect of modulating the seed source is investigated using the time-dependent approach. The results are in general agreement with the predictions of the steady-periodic approach and indicate that the two approaches contain similar underlying physics.

The paper is organized as follows: section 2 describes the governing equations. In section 3, we describe the analytical and numerical approaches to this phenomenon. The results of our simulations are presented in section 4 and a summary is provided in section 5.

## 2. Governing equations

Since we are examining the spatial and temporal evolution of the modes, we begin by considering the real part of the index of refraction to be a function of space and time, i.e. $n(\overrightarrow{r},t)$. The deviation of $n(\overrightarrow{r},t)$from the background index of refraction,${n}_{0}$, is attributed to thermal effects. The contribution of laser action to$n(\overrightarrow{r},t)$, which is dictated by the Kramer-Krönig relations, is negligible as argued in [12]. The pertinent wave equation for an isotropic medium derived from Maxwell’s equations is then given by:

In the analysis that follows, we consider a weakly guiding, polarization-maintaining (PM) fiber with the light launched into the fiber polarized along one of its axis; chosen here to be the x-direction. Efficient coupling to the other polarization state through modal interference is not possible and thus the other polarization state is not considered here. The longitudinal component of the electric field,${E}_{z}$, is small and is neglected in the analysis.

Since the laser gain and the variation in the index of refraction are small perturbations, we can construct the solution to the system as a superposition of the modes in the “ideal” fiber:

The variation in refractive index $\delta n=n-{n}_{0}$ is directly dependent on temperature difference, $\Delta T$, by $\delta n=\Delta T(dn/dT)$where $dn/dT=1.29\times {10}^{-5}$/°C is the thermo-optic coefficient for silica. Typically, for high-power amplifiers,$\Delta T$within the core is on the order of 100 °C. Therefore, ${(\delta n)}^{2}~O({10}^{-6})$while${n}_{0}~O(1)$, and the following approximation is valid: ${n}^{2}={({n}_{0}{}^{2}+\delta n)}^{2}={n}_{0}{}^{2}+2{n}_{0}\delta n+{(\delta n)}^{2}~{n}_{0}{}^{2}+2{n}_{0}\delta n$. The polarization in the x-direction due to laser gain is given by${P}_{LG,x}\approx 2i{\epsilon}_{0}{n}_{I}{n}_{0}{E}_{x}$, where ${n}_{I}$ is the imaginary index of refraction and ${n}_{I}\ll {n}_{0}$. We then proceed by relating ${n}_{I}$ to the imaginary part of the wave number,${k}_{I}$, which corresponds to the spatial growth (or decay) of the field due to the population inversion of the two-level system:

We proceed by using Eq. (3) in the wave equation for ${E}_{x}$ and applying the differential operators. The set of ${\phi}_{k}$ is comprised of the eigenfunctions of the “ideal” fiber and thus:

The differential operators on the left hand side (LHS) of Eq. (2) will also yield the following terms that need further examination to determine if they are retained in the first-order system of equations: $-(2{A}_{k}{\phi}_{k}/{n}_{0}){\partial}_{x}{}^{2}\delta n$, $-(2{A}_{k}/{n}_{0})({\partial}_{x}{\phi}_{k})({\partial}_{x}\delta n)$, $(4i{n}_{0}\omega /{c}^{2}){A}_{k}{\phi}_{k}{\partial}_{t}\delta n$, $-(4{n}_{0}/{c}^{2}){\phi}_{k}({\partial}_{t}{A}_{k})({\partial}_{t}\delta n)$, and $-(4{n}_{0}/{c}^{2}){\phi}_{k}{A}_{k}{\partial}_{t}{}^{2}\delta n$. The last three terms are clearly negligible as ${\partial}_{t}\delta n$is dictated by the thermal response time which is much greater than the optical period (i.e. ${\partial}_{t}\delta n\ll \omega \delta n$). Of the first two terms, one would expect the second term to provide the greater contribution. Even this term is typically negligible as $({\partial}_{x}{\phi}_{k})({\partial}_{x}\delta n)\ll {\beta}_{k}{}^{2}{\phi}_{k}\delta n$. This can be argued by noting that the length scale of the transverse profile is typically larger than the optical wavelength (i.e. ${\partial}_{x}{\phi}_{k}<{\beta}_{k}{\phi}_{k}$) and that ${\partial}_{x}\delta n$ which is dictated by the temperature profile in the transverse direction is even smaller as compared to ${\beta}_{k}\delta n$. For some of the simulations presented below, we have compared $({\partial}_{x}{\phi}_{k})({\partial}_{x}\delta n)/({\beta}_{k}{}^{2}{\phi}_{k}\delta n)$and determined it to be <$O({10}^{-4})$.

On the right hand side (RHS) of Eq. (2), only the second derivative in time term generating${\omega}^{2}$is retained. Terms containing ${\partial}_{t}({\sigma}_{s}^{(e)}{N}_{2}-{\sigma}_{s}^{(a)}{N}_{1})$ or ${\partial}_{t}^{2}({\sigma}_{s}^{(e)}{N}_{2}-{\sigma}_{s}^{(a)}{N}_{1})$are neglected as the laser gain response time is much greater than the optical period. The following equation that describes the evolution of the laser field is then obtained:

Equation (7) represents a system of coupled partial differential equations (PDEs) whose size is dictated by the number of modes. The standard approach to solving this system is to use the method of characteristics. In this case, the characteristics are defined by$dt/dz={n}_{0}/c$; leading to the characteristic lines given by $t-{t}_{0}=(z{n}_{0}/c)$ where ${t}_{0}$ initializes each set of characteristics. The right hand side (RHS) of Eq. (7) depends on time through$\delta n$,${N}_{2}$and _{${N}_{1}$}. The time scales for these quantities are much longer than the traverse time of the signal, $L{n}_{0}/c$ (where $L$is the length of the fiber). To be certain, for high-power amplifier experiments utilizing LMA fiber, the thermal response which dictates the variation in $\delta n$ in time is on the order of milliseconds and the population inversion occurs on a time scale of 10 *µ*sec. In contrast, the traverse time is on the order of 10 nsec. Therefore, we can estimate $t~{t}_{0}$ which is equivalent to freezing time as each set of light rays (characteristics) traverses the fiber. Therefore the equation describing the spatial evolution of ${A}_{j}$ at a time ${t}_{0}$ is given by:

The variation of $\delta n$ in time is captured by solving the time-dependent heat equation which is solved in conjunction with the coupled amplitude equations. Neglecting longitudinal heat flow, the heat equation takes the following form in polar coordinates:

Finally, the evolution of the pump intensity is given by:

## 3. Analytical and numerical approach

We simulate a co-pumped 1.6 m LMA Yb-doped fiber amplifier with 50 *μ*m core diameter, 250 *μ*m pump cladding diameter, and 500 *μ*m outer cladding diameter. The Yb doping concentration in the core is taken to be of 6 × 10^{25} /m^{3}. It is assumed that the fiber is “cold” at the beginning of the simulation and therefore evaluation starts immediately after the introduction of pump light into the fiber. Various seed powers ratios, and pump powers are examined and their effects on modal instabilities are reported in section 4. In this section we discuss the numerical and analytical details which are essential in understanding this phenomena as well as conducting accurate evaluations.

#### 3.1 The need for time-dependent treatment

One of the requirements for efficient transfer of energy is that the phase-matching condition is met. This is realized through the oscillations created in the temperature (and consequently the index of refraction) due to modal interference [5,6]. It was pointed out based on general arguments that, in addition, a phase lag between the signal intensity oscillations and the index of refraction oscillations is required [6]. We will take a different approach by showing mathematically using the population inversion and amplitude propagation equations that a steady state approach to the heat equation will largely lead to a phase shift of the amplitudes of the fields and no transfer of energy. Therefore, in order to simulate the modal instabilities and observe the energy transfer, the dynamical behavior of this phenomenon needs to be investigated.

In Eq. (11) for the generated heat, the term containing ${N}_{0}$is a DC term that does not provide for the required phase matching for energy transfer among the modes. The generated heat is also a function of the upper state population inversion,${N}_{2}$, which results in the dependence of the temperature and index of refraction variation on${N}_{2}$. This term will provide for oscillations with the proper periodicity. We restrict our analysis to the fundamental LP_{01} mode and the LP_{11} mode. In this case, ${I}_{s}$is given by:

_{01}mode, ${A}_{1}$and ${\phi}_{1}$are the amplitude and modal profile of the LP

_{11}mode, and $\Delta {\beta}_{01}={\beta}_{0}-{\beta}_{1}$. Substituting Eq. (15) into Eq. (12), one obtains the following expression for${N}_{2}$:

For brevity, we define$\xi ={A}_{0}{A}_{1}^{*}{\phi}_{0}{\phi}_{1}^{*}{e}^{i\Delta {\beta}_{01}z}+{A}_{0}^{*}{A}_{1}{\phi}_{0}^{*}{\phi}_{1}{e}^{-i\Delta {\beta}_{01}z}$. It can readily be shown that $\xi \le {\left|{A}_{0}\right|}^{2}{\left|{\phi}_{0}\right|}^{2}+{\left|{A}_{1}\right|}^{2}{\left|{\phi}_{1}\right|}^{2}$ leading to${\eta}_{2}\xi /{\alpha}_{2}<1$. By applying the binomial expansion, the following expression is obtained:

Equation (18) indicates that the most dominant term for the oscillations in ${N}_{2}$along $z$ is $\left(-\frac{{\alpha}_{1}{\eta}_{2}}{{\alpha}_{2}{}^{2}}+\frac{{\eta}_{1}}{{\alpha}_{2}}\right)\xi $ which oscillates with the period of $2\pi /\Delta \beta $. The term containing ${\xi}^{2}$does not oscillate at the required period with the next term providing this period being${\xi}^{3}$.

Note that $Q$ oscillates as ${N}_{2}$ does and is given by:

The term containing $\xi $ in the expression above accounts for the oscillation in the temperature. Consequently, the period of the oscillatory part of the average temperature is$2\pi /\Delta \beta $. The rest of the terms account for quantum defect heating in the absence of modal interference. It should be pointed out, however, that the temperature at the center of the core oscillates at ½ the period, i.e. $\pi /\Delta \beta $. This can be deduced by noting that $Q(r,\theta ,z)\approx Q(r,\theta +\pi ,z+\frac{\pi}{\Delta \beta})$and that from symmetry, the temperature at the center should be the same regardless if the higher value of $Q$at a given point in $z$ lies in the upper half ($0<\theta <\pi $) or the lower half ($0<\theta <\pi $) of the fiber.

We now proceed by examining the coupled system as described by Eq. (9). We examine, only the phase-matched terms as they are the only terms of relevance. Also, SPM due to the coupling coefficients ${\kappa}_{\delta n,0\to 0}$and ${\kappa}_{\delta n,1\to 1}$does not contribute directly to a transfer of energy and is neglected in this theoretical analysis. The function ${\phi}_{0}$is chosen to be real and based on symmetry arguments there is no loss of generality in also choosing in our analysis${\phi}_{1}$to be real. Since $\delta n$is a real quantity, ${\kappa}_{\delta n,1\to 0}$and ${\kappa}_{\delta n,0\to 1}$ are real and equal to each other. In general, one can show that the following expression holds:

Note that contribution of the term containing ${A}_{0}^{*}{A}_{1}{\phi}_{0}{\phi}_{1}{e}^{-i\Delta {\beta}_{01}z}$in the coupling coefficient as described by Eq. (20) is not included in the above system of equations as it does not provide phase matching. Also, note that ${g}_{l,0}$and ${g}_{l,1}$represent the laser gain for each mode in the absence of modal interference.

The terms containing $\varsigma $in Eqs. (21) are much smaller than the other terms. Therefore, the solution to system of Eqs. (21) can be approximated as:

Thus, the intensity of each mode is only affected by laser gain and there is no transfer of energy between them due to$\delta n$; i.e. only cross-phase modulation (XPM) happens. Equations (24) can be used in Eqs. (23) to obtain more explicit expressions for the amplitudes:

Consequently, a time-independent treatment of this phenomenon can only lead to a phase modulation of the amplitudes. In order for a significant energy transfer to occur, a phase lag between temperature and intensity is required which forms in time.

To illustrate how this phase lag develops, Fig. 1(a) displays the difference between core temperature and room temperature,$\Delta T$, along the fiber after 40 ms, 60 ms, and 80 ms for the fiber whose parameters are given above and at a pump power of 800 W. The numerical techniques used for the simulations are described below in subsection 3.2. The temperature inside the core is plotted at a fixed point away from the center and the close up is displayed in Fig. 1(b) to indicate the oscillations possess a period of approximately$2\pi /\Delta \beta $. The plot indicates that a phase lag occurs among the thermal gratings sampled at various times. In the inset to Fig. 1(b), we show that the oscillations at the center possess a period of $\pi /\Delta \beta $and in agreement with the conclusions given in [13].

#### 3.2 Numerical evaluation of the optical fields and temperature

We consider only the LP_{01} and LP_{11} in our numerical simulations. In order to solve the system of nonlinear ordinary differential Eqs. (9), we have applied the fourth-order Runge-Kutta method to compute the evolution of the optical fields along the fiber. By using a higher-order method, the solution converges with a lesser number of computational points. Note that the longitudinal length scale is dependent on the beat length; that is enough resolution is required to capture the beating between the two modes. For these computations, we used 110 computational points per beat length. The overlap integrals along the transverse direction required to solve the coupled system were evaluated to second-order using the trapezoidal rule. We have found that tracking the quantum efficiency of our system offered an effective tool to validate the numerical resolution of the grid in the longitudinal direction. Low resolution in the longitudinal direction led to a noticeable departure from the expected quantum efficiency of the system and consequently a lack of energy conservation. Unlike the analytical approach presented above, all terms in the nonlinear system of equations were retained: i.e. SPM, XPM, energy transfer as captured by${\kappa}_{\delta n,1\to 0}$, and direct coupling of the modes due to laser gain (i.e. term containing $\varsigma $). We noted that retention of the latter term is important in retaining high fidelity for conservation of energy.

The ADI method in polar coordinate was used to solve the heat equation. Due to the cylindrical shape of the fiber, it is more suitable to use polar coordinates to impose the correct boundary conditions which will result in smaller numerical error. The ADI method is typically implemented in rectangular domains to study the modal instabilities using a semi-periodic treatment of the heat equation and rectangular boundary conditions. The reason that this method is avoided in polar coordinates is the added complexity due to the occurrence of an artificial singularity at the center of the coordinate system. Nevertheless we have successfully implemented this method in polar coordinates.

In some studies of the modal instabilities, the boundary condition was artificially placed well inside the fiber. However, we believe the spatial mesh in the transverse direction must encompass a significant portion of the fiber for an accurate computation of the threshold of the modal instabilities. To illustrate this, we examine the temperature profile in a fiber of core radius of ${r}_{core}$ and outer radius of${r}_{out}$and a constant heat deposition${Q}_{0}$ inside the core. It is straightforward to determine the solution for a Dirichlet boundary condition to be:

We have chosen the spatial mesh in $r$ and $\theta $ so it minimized the computational cost while keeping the error in determining the temperature under 0.5%. The adaptive stepping along the radius was chosen by selecting a fine mesh inside the core and a larger spatial spacing for the pump cladding and outer cladding. For our simulations, 45 computational points in the radial direction and 90 in the azimuthal direction were used which summed up to 60 million computational points at each time step.

In order to avoid the singularity at the center of coordinate system, a grid was chosen so that the grid points are half-integrated in the radial direction and integrated in the azimuthal direction as given below:

**.**Using the periodicity of$\theta $, we also used $T(r,{\theta}_{0})=T(r,{\theta}_{m})$ and $T(r,{\theta}_{1})=T(r,{\theta}_{m+1})$. The implementation of half-integrating the mesh in the radial direction resulted in cancelling the terms with ${r}_{1}$ which provided us with all the boundary conditions required for solving the system. Note that the matrices formed here are nonsingular and the inverses can be pre-computed.

Figure 2 displays the convergence of this method for a test problem with a known solution. This example uses$Q=(r\mathrm{cos}\theta +r\mathrm{sin}\theta -2{t}^{2})\mathrm{exp}\left(rt(\mathrm{sin}\theta +\mathrm{cos}\theta )\right)$. Here $n=8,16,32$, $m=2n$, and ${k}_{t}=80n$ where ${k}_{t}$ is the time step. Table 1 and Fig. 2 confirm that the method converges with second-order accuracy as expected.

## 4. Numerical results

Our numerical and analytical work has confirmed the thermal origin of modal instabilities as well as offering robust explanation on dynamical dependence of temperature effects. We believe the longitudinal dependence of heat is negligible due to the comparison of longitudinal length scale versus transverse length scale. The variation of temperature in the transverse direction is much larger than its variation in a beat length. Phase shift occurs in time without considering longitudinal heat. We have not considered any frequency offset between the two modes. Instead the time-dependent evolution has been studied whereby the Doppler shift is automatically considered. Furthermore, the polar coordinate has been used in our model where correct boundary condition has been implemented. Our numerical results are presented below. Various seed powers and ratios are examined and their effects on modal instabilities are reported. In addition, two different gain-tailored designs for modal instabilities suppression are simulated and compared. Furthermore, we investigate the effect of signal modulation on the threshold.

#### 4.1 Power-dependent behavior

An important characteristic of the modal instabilities is power dependence and as reported by Otto et al. [11] consists of three regimes: stable, transition, and chaotic. Figure 3(a) displays the total intensity in the two modes at the output end of the fiber for various times for a seed power of 30 W of which 5% is in the LP_{11} mode. Also shown, is the power in the two modes at the output end of the fiber as a function of time (up to 80 ms). The amplifier is co-pumped with a pump power of 720 W (below threshold). Temporal variation of intensity at the output end as well as temporal variation in modal powers along the propagation direction,$z$, is shown in Media 1. Changes occur inside the fiber below the power threshold but they are not strong enough to impact the fundamental mode after a certain time. The energy transfers from the fundamental mode to the higher-order mode back and forth, but after the temperature approaches a quasi-equilibrium phase, energy remains largely in the fundamental mode. A similar behavior is shown in Ward et al. [10]. Figure 3(b) and Media 2 show the corresponding plots for a pump power of 820 W (just past threshold). In this case, the modes keep transferring energy and no steady state is reached. Increasing power further leads to an increase in the frequency of “deep” oscillations between the two modes and this instability enters a “chaotic” regime. In this regime, the power in each mode, on average, is ~50%. Figure 3(c) and Media 3 display the corresponding plots for a pump power of 1040 W. It should be pointed out that, in all three cases, there are transient effects occurring in the first 10-20 ms related to the abrupt turn on of the pump power.

We have computed the ratio of the average power over a “long time” range in the LP_{11} mode to the total signal power as a function of pump power. The “long time” range in our simulations was 60-80 ms. The higher-order mode content vs. pump power is shown in Fig. 4. The modal instabilities threshold represents the pump power where a sudden increase in the average power of the higher-order mode is observed. The region below the threshold is referred to in the figure as region I. Past the threshold, there exists a pump power range (region II) whereby quasi-regular temporal oscillations in the modal content occur and the average power in the higher-order mode can exceed 50% of the total signal power. In region III (chaotic regime), rapid temporal oscillations occur with almost the entire power in one mode transferring to the other mode and back and forth. These numerical simulations are in general agreement with the experimental results reported by Otto et al. [11] (refer to Fig. 4 of reference).

#### 4.2 Seed-power dependence

The work of Wirth et al. reports a 50% increase in modal instabilities threshold using a seed power of 200 W instead of 27 W [14]. To investigate this observation numerically, we varied the seed power from 25 to 175 W in increments of 25 W. The power in the LP_{11} was set at 5% of the total seed power. We estimated the threshold by examining the modal content as a function of pump power and determining the point at which the higher-order mode power rises rapidly. Figure 5 displays the threshold as a function of seed power. We have found a linear relation over the range of seed powers investigated between the modal instabilities threshold and the seed power. The increase in threshold for a seed power changing from 25 W to 175 W is ~20%. The experimental work in [14] reports a more dramatic effect. However, the fiber used in the experiments is different than the fiber in our simulations. Also, as reported in that reference, the degree of polarization in their experiments varies as the output power changes. The effect of polarization can be one of the reasons that they observe an improvement of ~50%.

Furthermore, it has been proposed that the seed power responsible for the modal instability can be as small as quantum noise [6]. We studied the change in threshold by using different seed powers ratios in each mode in order to investigate how sensitive the threshold is to the amount of power injected into the higher-order mode. Note that here the seed power was fixed to 30 W and only the ratio of power in each mode was varied. The power injected in the higher-order mode was varied between 10^{-7%} of the seed power and all the way up to 10%. Threshold for all cases was determined to lie in a narrow range ~750-770 W. Therefore, the threshold is insensitive to the seeding ratio; i.e. injecting more power into the higher-order mode has little effect. Figure 6 shows the modal content at the output of the fiber vs. time for a pump power of 760 W, and a seed power of 30 W with 10^{-7%} injected in the LP_{11} mode. At this pump power, we observe the beginning of an appreciable amount of energy being transferred to the higher-order mode. By pumping at slightly higher power, the amplifier will then firmly enter region II as displayed in Fig. 4.

#### 4.3 Gain-tailored design

The work of Robin et al. [15] reports on the observation of modal instabilities in a segmented acoustically tailored (SAT) photonic crystal fiber (PCF) amplifier. The purpose of the acoustic tailoring was to suppress SBS by creating two acoustic regions in the core: one at the center of the core and another in the six hexagons surrounding it [16]. While this design was effective in suppressing SBS, the threshold for the modal instabilities was encountered at ~530 W. A second generation SAT fiber whereby the Yb-ion concentration was tailored to have preferential overlap with the fundamental mode and reduced pump absorption was designed and fabricated [15]. The gain tailoring was achieved by ensuring that three non-adjacent hexagons lying in the outer ring were Yb-free as shown in Fig. 7. Since the fabrication of PCFs is performed using the draw and stack technique, this gain-tailored design is well-suited for PCFs.

The results from the gain-tailored SAT fiber were promising and modal instabilities were mitigated appreciably. The modal instabilities threshold is increased to >990 W . The numerical model presented here is used to investigate the effect of gain-tailored designs. Two different designs have been investigated numerically and been compared to the uniformly Yb-doped core fiber. For both cases, the core and cladding dimensions are the same as the fiber described above. The first design, design A, is a similar design to the one shown in Fig. 7. The design on the numerical grid of the core is shown in Fig. 8(a). White area displays the Yb-free parts of the core. The ytterbium is removed from $\theta =\pi /3$to$\theta =2\pi /3$, $\theta =\pi $to$\theta =4\pi /3$, and $\theta =5\pi /3$to $\theta =2\pi $for $R>{R}_{1}$ where ${R}_{1}$is the radius of the inner core segment. Therefore, each doped and un-doped region in the outer ring of the core has the same area. The lengths of the gain tailored fibers were adjusted to provide for the same level of pump absorption.

Numerically, ~20% increase in threshold is observed for design A. Therefore, the modal instabilities were mitigated using this design but this suppression is less than the increase in the threshold observed for the gain-tailored SAT fiber used in the experiments described in [15]. We note that the simulated fiber and the fiber used in the experiments have different dimensions. In addition, the Yb doping concentration is much lower at 3.5 × 10^{25}/m^{3} for the latter fiber than the 6 × 10^{25}/m^{3} concentration used throughout our simulations. Furthermore, in the experiments, the threshold of modal instabilities is a difficult parameter to characterize due to large dependency of this threshold to system parameters which can lead to some discrepancy between experimental and theoretical results.

For design B, Yb is not present in a ring on the outer edge of the core. This design is suitable for fabrication through modified chemical vapor deposition (MCVD) process. A total of 70% of core area is doped in this design. This design is shown in Fig. 8(b) and provides a promising result of ~30% suppression of the modal instabilities threshold. Figure 9 displays the modal content of the higher-order mode averaged between 60 to 80 ms as the pump power increases for the three cases: uniformly doped core, design A, and design B. Overall, gain tailoring is an effective technique to suppress the instabilities and further investigation of various designs to optimize the suppression is warranted.

#### 4.4 Seed modulation

We have studied the effect of seed modulation on the modal instabilities threshold and compared our time-dependent results to the work of Hansen et al. [8]. One major difference is that we do not assume in our computations a frequency offset as is assumed by the quasi-periodic models [6–9]. We assume a sinusoidal modulation such that the signal power,${P}_{s}$, is given by:

where $a$ is the modulation depth, $\Omega $ is the modulation angular frequency, and ${P}_{s}{}^{0}$is the average power over the modulation frequency. Note in comparing to [8], the modulation depth here is ½ that used in the reference. In this set of simulations, we have taken the same uniformly doped fiber parameters used in the results presented above. At a modulation frequency of 1 KHz, three modulation depths were considered: $a=1/3000$,$a=1/30000$, and $a=1/300000$. The corresponding thresholds were computed to be 350W, 400W, and 450, respectively. The plot of higher-order modal content vs. pump power for the three cases is shown in Fig. 10. These results are in general agreement with results presented in [8] and reveal that a full time-dependent approach, without the assumption of a frequency offset a priori, recovers one major conclusion of the quasi-periodic approach.We also examined for the case $a=1/3000$, the modal instabilities threshold as a function of modulation frequency. Four cases were numerically investigated: $\Omega =$0.5 kHz, 1 kHz, 1.5 kHz, and 2 kHz. We found that the highest gain for the instability for the frequencies investigated occurred at $\Omega =$1 kHz, while the lowest occurred at$\Omega =$2 kHz. This would indicate that gain peak lies in the range between 0.5 kHz and 1.5 kHz. Figure 11 displays the higher-order mode content vs. pump power for the four modulation frequencies investigated.

## 5. Summary

Using an analytical approach, we showed that a steady state solution to the heat equation could not lead to any appreciable energy transfer and a time-dependent study is required. We presented detailed numerical simulations of modal instabilities in Yb-doped fiber amplifiers using a time-dependent temperature solver in polar coordinates coupled into the optical fields and population inversion equations. Numerical convergence of the solver was presented and the importance of correct implementation of boundary conditions and computational window size was discussed. The accuracy of the spatial grid for the integration of the optical fields along the direction of propagation was confirmed trough numerical convergence and by tracking the quantum efficiency of the system as related to energy conservation. The model captures the three power-dependent regions that are characteristic of the dynamics of modal instabilities. Using this model, we found a linear relationship between threshold and seed power over a wide range. Also, we showed that the threshold is insensitive to the input modal content. Mitigation of the modal instabilities through gain tailoring was investigated by simulating two different fiber designs including a design that has been implemented experimentally in a PCF. Furthermore, we have shown that this fully time-dependent approach agrees with the quasi-periodic approach regarding the effect of signal modulation on the threshold. Finally, although this time domain approach is computationally intensive, it captures: 1) the whole range of offset frequencies that can contribute to the transfer of energy, 2) laser gain both spatially and temporally. The former has yet to be included in the quasi-periodic fully numerical approach (this approach assumes one offset frequency at a time), while the latter is not included in the semi-analytical approach that integrates over the frequency range.

## Acknowledgments

The authors would like to thank the High Energy Laser Joint Technology Office (HEL-JTO) for partial funding of this work. This work was supported in part by a grant of computer time from the DoD High Performance Computing Modernization Program. Dr. Naderi performed this research while serving as a postdoctoral fellow under the National Research Council Research Associateship program.

## References and links

**1. **C. X. Yu, S. J. Augst, S. M. Redmond, K. C. Goldizen, D. V. Murphy, A. Sanchez, and T. Y. Fan, “Coherent combining of a 4 kW, eight-element fiber amplifier array,” Opt. Lett. **36**(14), 2686–2688 (2011). [CrossRef] [PubMed]

**2. **F. J. Kontur, I. Dajani, Y. Lu, and R. J. Knize, “Frequency-doubling of a CW fiber laser using PPKTP, PPMgSLT, and PPMgLN,” Opt. Express **15**(20), 12882–12889 (2007). [CrossRef] [PubMed]

**3. **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). [CrossRef] [PubMed]

**4. **F. Stutzki, H.-J. Otto, F. Jansen, C. Gaida, C. Jauregui, J. Limpert, and A. Tünnermann, “High-speed modal decomposition of mode instabilities in high-power fiber lasers,” Opt. Lett. **36**(23), 4572–4574 (2011). [CrossRef] [PubMed]

**5. **C. Jauregui, T. Eidam, H.-J. Otto, F. Stutzki, F. Jansen, J. Limpert, and A. Tünnermann, “Temperature-induced index gratings and their impact on mode instabilities in high-power fiber laser systems,” Opt. Express **20**(1), 440–451 (2012). [CrossRef] [PubMed]

**6. **A. V. Smith and J. J. Smith, “Mode instability in high power fiber amplifiers,” Opt. Express **19**(11), 10180–10192 (2011). [CrossRef] [PubMed]

**7. **K. R. Hansen, T. T. Alkeskjold, J. Broeng, and J. Lægsgaard, “Thermally induced mode coupling in rare-earth doped fiber amplifiers,” Opt. Lett. **37**(12), 2382–2384 (2012). [CrossRef] [PubMed]

**8. **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). [CrossRef] [PubMed]

**9. **L. Dong, “Stimulated thermal Rayleigh scattering in optical fibers,” Opt. Express **21**(3), 2642–2656 (2013). [CrossRef] [PubMed]

**10. **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). [CrossRef] [PubMed]

**11. **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). [CrossRef] [PubMed]

**12. **K. R. Hansen, T. T. Alkeskjold, J. Broeng, and J. Lægsgaard, “Thermo-optical effects in high-power ytterbium-doped fiber amplifiers,” Opt. Express **19**(24), 23965–23980 (2011). [CrossRef] [PubMed]

**13. **C. Jauregui, T. Eidam, H.-J. Otto, F. Stutzki, F. Jansen, J. Limpert, and A. Tünnermann, “Temperature-induced index gratings and their impact on mode instabilities in high-power fiber laser systems,” Opt. Express **20**(1), 440–451 (2012). [CrossRef] [PubMed]

**14. **C. Wirth, T. Schreiber, M. Rekas, I. Tsybin, T. Peschel, R. Eberhardt, and A. Tünnermann, “High-power linear-polarized narrow linewidth photonic crystal fiber amplifier,” Proc. SPIE **7580**, 75801H, 75801H-6 (2010). [CrossRef]

**15. **C. Robin, I. Dajani, C. Zeringue, B. Ward, and A. Lanari, “Gain-tailored SBS suppressing photonic crystal fibers for high power applications,” Proc. SPIE **8237**, 82371D, 82371D-10 (2012). [CrossRef]

**16. **C. Robin and I. Dajani, “Acoustically segmented photonic crystal fiber for single-frequency high-power laser applications,” Opt. Lett. **36**(14), 2641–2643 (2011). [CrossRef] [PubMed]