Transverse Mode Instability (TMI) that results from dynamic nonlinear thermo-optical scattering is the primary limitation to power scaling in high-power fiber lasers and amplifiers. It has been proposed that TMI can be suppressed by exciting multiple modes in a highly multimode fiber. We derive a semi-analytic frequency-domain theory of the threshold for the onset of TMI in narrowband fiber amplifiers under arbitrary multimode input excitation for general fiber geometries. Our detailed model includes the effect of gain saturation, pump depletion, and mode-dependent gain. We show that TMI results from the exponential growth of noise in all the modes at downshifted frequencies due to the thermo-optical coupling. The noise growth rate in each mode is given by the sum of signal powers in various modes weighted by pairwise thermo-optical coupling coefficients. We calculate thermo-optical coupling coefficients for all $\u223c104$ pairs of modes in a standard circular multimode fiber and show that modes with large transverse spatial frequency mismatch are weakly coupled, resulting in a banded coupling matrix. This short-range behavior is due to the diffusive nature of the heat propagation, which mediates the coupling and leads to a lower noise growth rate upon multimode excitation compared to a single mode, resulting in significant TMI suppression. We find that the TMI threshold scales linearly with the number of modes that are excited asymptotically, leading to roughly an order of magnitude increase in the TMI threshold in an 82-mode fiber amplifier.

## I. INTRODUCTION

Fiber lasers based on multi-stage fiber amplifiers (FAs) provide an efficient and compact platform to generate ultra-high laser power.^{1–6} This has enabled potential applications in a wide range of technologies, such as laser welding,^{7} gravitational wave detection,^{8} and directed energy.^{9} To fully realize this potential, further power scaling is needed in the FAs,^{10} which has been limited primarily by nonlinear effects, such as transverse mode instability (TMI)^{11–14} and stimulated Brillouin scattering (SBS).^{15–17} TMI is a dynamic transfer of power between the transverse modes of the fiber caused by thermo-optical scattering and/or inversion fluctuations.^{18–22} During the optical amplification process, heat is generated due to the quantum defect, in an amount proportional to the local optical intensity, creating local temperature fluctuations. The local temperature variations then create refractive index variations, causing significant scattering between the modes. Consequently, when the FA is operated above a certain output power, defined as the TMI threshold, a significant degradation of beam quality occurs,^{23–33} rendering the output unsuitable for many applications. As a result, suppressing TMI or equivalently raising the TMI threshold has been one of the most important and technologically relevant scientific goals in the high-power laser community.^{4,5,11}

Significant efforts have been undertaken to mitigate TMI, such as dynamic seed modulation,^{34} synthesizing materials with low thermo-optic coefficient,^{35,36} modulating the pump beam,^{37,38} increasing the optical loss of higher-order modes (HOMs),^{39–44} and utilizing multicore fibers.^{45–47} Although most of these efforts have had some success, they suffer from one or more drawbacks, such as fluctuating output power in the case of seed or pump modulation, difficulty with mass-manufacturing custom fibers, and increased guidance of HOMs due to fiber heating, rendering the HOM suppression unfeasible. As such, efficient TMI suppression remains a highly active area of research.

A common feature in all of the previous approaches is to excite the fundamental mode (FM) of the fiber as much as possible, hence reducing the excitation and amplification of HOMs.^{34–46,48} This is motivated by the widespread impression that the speckled internal field generated by multimode excitation will necessarily create a poor beam quality at the output.^{49,50} However, recent progress in wavefront shaping, enabled by spatial light modulators (SLMs), has demonstrated that multimode excitation with a narrowband seed laser does not inherently lower the beam quality as long as the light remains spatially coherent.^{51–53} In fact, by using an SLM to wavefront-shape the input light to a fiber, it is possible to obtain a diffraction-limited spot after coherent multimode propagation in both passive^{54} and active fibers,^{51} which can be easily collimated using a lens. This enables a fundamentally novel approach to control nonlinear effects in fibers by utilizing selective multimode excitation.^{55–58} The viability of such an approach has recently been demonstrated for other nonlinear effects in fibers, such as SBS^{59–61} and stimulated Raman scattering (SRS).^{62}

The current authors have recently showed, using numerical simulations along with some initial theoretical results, that sending power in multiple fiber modes robustly raises the TMI threshold.^{63} To investigate this approach thoroughly and in detail, we have developed a theory of TMI for arbitrary multimode input excitations and general fiber geometries, building on significant prior theoretical efforts to model TMI, which have been successful in capturing much of the key physics.^{14,18,19,22,26,31,64–71} It has been shown that TMI can be modeled as the exponential growth of noise in the HOMs due to the thermo-optical scattering of the signal. However, as noted, almost all of the previous efforts have focused on the case when only the FM of the fiber is excited with the presence of noise in a few HOMs. None of these works attempted to model highly multimode excitations, and their formalisms were not suited for such explorations. As mentioned, some key theoretical results for TMI threshold upon multimode excitation were presented in Ref. 63 by the authors of this paper, including results demonstrating linear scaling of the TMI threshold with the number of excited modes. This previous work focused on time-domain numerical simulations, and the detailed theoretical formalism for TMI threshold along with the analysis of pump depletion and gain saturation for fibers with arbitrary geometries and input excitations was deferred to the current work.

In this paper, we derive a general formalism that allows efficient calculation of the TMI threshold for arbitrary multimode input excitations and fiber geometries in narrow-linewidth fiber amplifiers. We derive coupled amplitude equations starting from coupled optical-propagation and heat-diffusion equations and solve them to obtain analytic formulas for the TMI gain and thermo-optical coupling coefficient. A key result of our theory is that the thermo-optical coupling is strong only between the optical modes that have similar transverse spatial frequencies. This is a result of the diffusive nature of the heat propagation, which underlies the thermo-optical coupling. This leads to a banded thermo-optical coupling matrix, due to which the effective TMI gain is significantly lowered upon multimode excitation, resulting in robust TMI suppression. We show that the TMI threshold increases linearly with the number of excited modes, leading to more than an order of magnitude higher TMI threshold in standard multimode fibers. We also show that threshold enhancement can be further increased by tailoring the fiber geometry.

We begin by deriving a semi-analytic solution for the optical and heat equations, coupled via the nonlinear-polarization and quantum-defect-heating terms. We expand the optical fields and the temperature fluctuations in terms of eigenmodes of the optical and heat equations, respectively. From the driven heat equation, the coefficient of each temperature eigenmode is found in terms of optical amplitudes, leading to coupled amplitude equations for the optical modes, resembling those for a four-wave-mixing process.^{72} Next, we solve these coupled amplitude equations for an arbitrary input signal at frequency, *ω*_{0}, along with the presence of small amount of noise at Stokes shifted frequencies, *ω*_{0} − Ω, in each mode. We show that the noise power in various modes grows exponentially and can lead to dynamic spatial fluctuations in the output beam profile after sufficiently high noise growth. The growth rate of the noise power in a given mode is equal to the sum of the signal power in each of the other modes, weighted by pairwise thermo-optical coupling coefficients. For any pair of optical modes, the thermo-optical coupling coefficient is proportional to the overlap integrals of the two optical mode profiles with various eigenmodes of the heat equation. The coupling between the optical modes with a large separation in transverse spatial frequency is mediated by temperature modes with high transverse spatial frequency; these thermal modes necessarily have a very short-lifetime, leading to a quite weak coupling between such mode pairs. To study this explicitly, we calculate the thermo-optical coupling coefficients for all $\u223c104$ mode pairs in a commercially available highly multimode (MM) circular step-index fiber (supporting 82 modes per polarization) and find that the matrix of peak values of pairwise thermo-optical coupling coefficients is a banded, effectively sparse matrix. This sparse/banded nature of the thermo-optical coupling is not affected by the presence of gain saturation or pump depletion. As a result, the noise growth rate goes down linearly with the number of modes excited, if the power is equally (or nearly equally) distributed in all the excited modes. Therefore, the TMI threshold is predicted to increase linearly with the number of equally excited modes, which is verified by explicit calculation of the TMI threshold.

The generality of our approach allows us to do similar calculations for different fiber cross sections (e.g., square and D-shaped); we find that the linear increase of the threshold is generic and only weakly sensitive to fiber and cladding geometry. Our theory is presented first with a simpler model, which neglects the effects of gain saturation and pump depletion, but highlights the universality of the physics due to the mismatch of spatial and temporal scales (thermal vs optical), which mediates TMI suppression via multimode excitation. In Sec. III of this paper, we present an improved model, which does include these important effects (at a tractable computational cost), and confirms that the asymptotic linear scaling of the TMI threshold is preserved with a significant but modest increase in the calculated TMI threshold. Taken in full, our results quantitatively confirm that highly multimode excitations efficiently suppress TMI, opening up a new platform for robust power scaling in high-power fiber laser amplifiers.

## II. THEORY AND RESULTS

### A. Coupled amplitude equations

^{11,18,19}A schematic of the TMI process is shown in Fig. 1. A signal wave at frequency

*ω*

_{0}is sent through an active fiber to undergo optical amplification and generate the output. However, due to imperfect coupling or other experimental artifacts, there usually exists a small amount of noise at shifted frequencies,

*ω*

_{0}− Ω, in various modes. The optical signal and noise can interfere to create a spatiotemporally varying optical intensity pattern, which, due to the quantum-defect heating that accompanies stimulated emission, results in dynamic temperature variations.

^{24}These spatiotemporal thermal fluctuations result in refractive index variations, which can cause significant transfer of power from the signal at

*ω*

_{0}in one mode to the noise at Stokes-shifted frequencies

*ω*

_{0}− Ω in other modes. As a result, at high-enough output power, the noise in various modes can become a significant fraction of the output power, leading to fluctuations in the output beam profile. This output power level is defined as the TMI threshold

^{18,19}and marks the onset of an unstable regime for the output. To derive a formalism to calculate the TMI threshold, we solve the optical wave equation with a temperature-dependent nonlinear polarization as a source term,

^{19}

^{,}

*c*is the speed of light, and

*n*

_{0}is the linear refractive index of the fiber. The temperature-dependent polarization is given by

*η*is the thermo-optic coefficient of the fiber material and Δ

*T*represent the local temperature fluctuation due to the quantum-defect heating. Δ

*T*satisfies the heat equation with a source term proportional to the local intensity,

^{18}

^{,}

*Q*is proportional to the local optical intensity

*I*, the amplifier gain coefficient

*g*, and the quantum defect $qD=(\lambda s\lambda p\u22121)$, which depends on the difference between the signal wavelength

*λ*

_{s}and pump wavelength

*λ*

_{p}. The thermal conductivity, the specific heat, and the density of the material are given by

*κ*,

*C*, and

*ρ*, respectively. The ratio

*κ*/

*ρC*is equal to the diffusion constant

*D*. Since the thermo-optical polarization depends on the temperature, which in turn depends on the optical intensity, the optical and heat equations are coupled and nonlinear due to the source terms on the right hand side (RHS) of Eqs. (1) and (3). For a translationally invariant system, such as an optical fiber, both optical and heat equations can be solved formally by expanding the fields in terms of the eigenmodes of the linear operator corresponding to each equation. As usual, we decompose the electric field $E\u20d7$ as a sum of the product of a transverse mode profile, a rapidly varying longitudinal phase term, and a slowly varying amplitude (SVA) for each fiber mode,

*A*

_{m}(Ω,

*z*) is the SVA in mode

*m*at a point

*z*along the fiber axis for a Stokes frequency shift Ω from the central frequency

*ω*

_{0}. The Ω = 0 amplitude corresponds to the signal, and Ω ≠ 0 corresponds to the noise. $\psi \u20d7m$ and

*β*

_{m}denote the transverse mode profile and propagation constant for mode

*m*and can be obtained by solving the fiber modal equation,

^{73,74}

^{75}for an arbitrary fiber cross-sectional geometry. As an example, in Fig. 2(a), we have shown the electric field profiles for the FM and a HOM for a circular step-index fiber. The width of the core is 40

*μ*m, and the cladding width is 200

*μ*m. The FM has no node in the electric field profile and follows the symmetries of the cross-sectional geometry. Each HOM can be labeled as a linearly polarized (LP) mode characterized by two indices (

*u*,

*v*) in the standard notation.

^{73}The mode LP

*uv*has

*u*azimuthal nodes and (

*v*− 1) radial nodes, with a radial profile given by the (

*u*+ 1)th Bessel function of the first kind, with (

*v*− 1) zeros within the fiber core and angular profile given by either a sine or cosine function with

*u*zeros. Figure 2(a) shows the profiles of LP01 (FM) and LP44 (HOM) modes. In the Appendix, we study fibers with different geometries (D-shaped and square cross sections), and the associated optical modes are shown in Fig. 7(a) .

*ψ*

_{m}to have unit power, the total power in mode

*m*at point

*z*along the fiber axis and at Stokes frequency Ω is given by $|Am(z,\Omega )|2$. The heat source term

*Q*in Eq. (3) can be obtained in terms of optical modal amplitudes using Eqs. (4) and (5). The optical intensity

*I*contains terms corresponding to the interference of various fiber modes, and thus,

*Q*(∝

*I*) is given by a sum of dynamic heat sources oscillating at Stokes frequency Ω and propagation constant

*q*

_{ij}=

*β*

_{j}−

*β*

_{i}for each mode pair {

*i*,

*j*}. Therefore, we expand Δ

*T*into a series of temperature profiles corresponding to each heat source term,

*q*

_{ij}and temporal frequency Ω. We further decompose each temperature profile

*T*

_{ij}in terms of the eigenmodes of the transverse Laplacian $\u2207T2$,

*k*in temperature grating {

*i*,

*j*}. Temperature eigenmodes of the fiber oscillate in space and relax exponentially in time due to diffusive heat equilibration. Their decay rate is proportional to the thermal diffusion constant and the wavevector of the fluctuations squared. Hence, long wavelength fluctuations dominate the thermal response; as we will see, this directly leads to a nonlinear thermo-optic coupling, which is short-range in the momentum

*difference*between the optical modes. $T\u0303k$ denotes the transverse profile of temperature eigenmode

*k*, which satisfies the following eigenvalue equation:

*k*. We impose Dirichlet boundary conditions at the outer surface of the fiber, corresponding to the standard assumption

^{19}that the outer surface is held at a constant temperature, and all the heat immediately dissipates into a heat bath. The temperature eigenmodes can be solved either analytically or numerically depending on the fiber geometry. We have shown in Fig. 2(b) the fundamental and higher-order temperature eigenmodes for a fiber with circular cladding, which are calculated using the Coefficient-Form-PDE module in COMSOL.

^{75}The spatial structure of the temperature eigenmodes is similar to the optical modes since both are eigenfunctions of the same spatial operator $\u2207T2$. An important distinction is that optical modes are localized in the fiber core, but the temperature eigenmodes spread out over the entire fiber cross section due to the differing boundary conditions. The fundamental mode has no nodes and follows the symmetries of the cladding geometry. The higher-order eigenmodes are characterized by

*u*azimuthal and

*v*radial nodes across the entire fiber cross section (i.e., core plus cladding), with a radial dependence given by the (

*u*+ 1)th Bessel function of the first kind and angular dependence given by sine or cosine with

*u*zeros.

*T*, we need to determine each coefficient $aijk$. This is done by simplifying the left hand side (LHS) of Eq. (3) by substituting Δ

*T*from Eqs. (7) and (8) and using the eigenmode equation [Eq. (9)]. The RHS of Eq. (3) can be obtained in terms of optical amplitudes by using Eqs. (4) and (5). Finally, by exploiting orthogonality, we isolate the desired coefficient $aijk$ by multiplying both sides of simplified Eq. (3) with $T\u0303k*$ and integrating it across the fiber cross section, leading to

*k*for the heat source term {

*i*,

*j*} is proportional to the overlap of the dot product of $\psi \u20d7i$ and $\psi \u20d7j$ with the temperature mode profile $T\u0303k$ integrated over the fiber cross section, denoted by the angular brackets ⟨·⟩. It is also proportional to the convolution of optical amplitudes

*A*

_{i}and

*A*

_{j}. We introduce an inverse mode lifetime for mode

*k*given by $\Gamma kij=((\beta i\u2212\beta j)2+\alpha k2)D$ with units of [

*s*

^{−1}]. The first term in the inverse mode lifetime corresponds to the longitudinal heat diffusion, and the second term corresponds to the transverse heat diffusion. In typical fibers, the transverse heat diffusion dominates, so the longitudinal heat diffusion can be ignored.

^{18,19}In that case, we can define $\Gamma kij\u2248\Gamma k=\alpha k2D$, which we use in the rest of this paper. However, in fibers where longitudinal diffusion is significant, both contributions need to be included in the inverse mode lifetime. This completes the formal solution for Δ

*T*in terms of the optical mode profiles and amplitudes.

*T*to evaluate the source terms in Eq. (1), which can be used to derive coupled self-consistent equations for optical amplitudes. We use the ansatz for $E\u20d7$ from Eq. (5) to simplify the left hand side (LHS) of Eq. (1). Finally, we use the orthogonality of the optical modes to isolate the equation of each optical amplitude {

*A*

_{m}} by multiplying both sides of the simplified equation [Eq. (1)] with $\psi \u20d7m*$ and integrating it over the fiber cross section,

*z*derivative of the modal amplitude

*A*

_{m}has contributions from both the linear amplification term with growth rate

*g*/2 due to stimulated emission and a nonlinear four-wave-mixing term

^{72}due to the thermo-optical interaction, which is proportional to a susceptibility constant

*χ*

_{0}and the sum of convolutions of three modal amplitudes

*A*

_{n}, $Ai*$, and

*A*

_{j}. The exponential term tracks the phase mismatch for each contribution, and the strength of each term is proportional to the overlap of Green’s function of the heat equation with the relevant optical mode profiles,

*k*represents the sum over all the temperature eigenmodes. Thus,

*G*

_{mnij}(Ω) denotes the optical modal overlap with Green’s function of the heat equation written in the spectral representation.

^{76}The susceptibility constant

*χ*

_{0}involves a combination of various material and optical constants involved and has dimensions of [

*W*

^{−1}].

*G*

_{mnij}(Ω) has both real and imaginary parts, in general. The real part is responsible for the nonlinear phase evolution of

*A*

_{m}, while the imaginary part is responsible for dynamic transfer of power between various modes. Clearly, for Ω = 0, the imaginary part of

*G*

_{mnij}vanishes, and thus, there is no transfer of power between various modes of the fiber at

*ω*

_{0}. A non-zero Stokes frequency shift is needed for TMI, which is in agreement with the well-known result that a moving intensity grating (Ω ≠ 0) is required for TMI.

^{11,18,24,42}Typically, there is amplitude noise in the fiber, either due to experimental imperfections or due to spontaneous emission, having multiple frequency components. The noise frequency at which the imaginary part of Green’s function peaks undergoes the highest amplification and is typically a few kHz, leading to dynamic fluctuations in the beam profile on the order of milliseconds.

^{11,18,23}

For any given set of input amplitudes in each mode, {*A*_{m}(*z* = 0, Ω)}, Eq. (11) can be solved to obtain modal amplitudes {*A*_{m}(*z*, Ω)} everywhere, fully determining both $E\u20d7$ and Δ*T*. The fiber properties are taken into account through *χ*_{0} and *G*_{mnij} via optical mode profiles {*ψ*_{m}} and temperature mode profiles ${T\u0303k}$. Equation (11) is a set of highly coupled nonlinear differential equations and, in general, cannot be solved analytically. Numerical solutions are possible by using finite difference methods^{77} to discretize the derivative operator and evaluating the convolutions by either built-in or custom operations. It is expected that numerically solving Eq. (11) would be more computationally efficient than solving the original coupled optical and heat equations. This is because Eq. (11) is a set of 1D equations, where transverse degrees of freedom are accounted through the fiber modes, which exploit the longitudinal translation invariance and only need to be solved once for a given fiber. However, despite this simplification, when a large number of modes are present, the computational complexity can quickly become very high. Additionally, for studying TMI suppression using multimode excitation, the input modal content needs to be parametrically tuned, requiring a fast and efficient solution. Fortunately, an approximate solution to the coupled amplitude equations can be obtained that captures the essential features of the onset of TMI and can be used to calculate the TMI threshold, which is derived in Sec. II B.

### B. Phase-matched noise growth

^{18,19,31,71}since phase-mismatched terms become insignificant over long enough length scales, and (2) we assume that the noise power is significantly lower than the signal power and the change in the signal due to the noise growth is ignored.

^{18,19,69}This is typically valid below the TMI threshold, which marks the onset of significant beam fluctuations, and can be used to calculate the threshold. We assume that at the fiber input, the signal (or, seed) is injected at frequency

*ω*

_{0}in various modes with complex amplitudes ${Ams(z=0)}$. In addition, there is noise present at Stokes shifted frequencies

*ω*

_{0}− Ω in all the fiber modes, denoted by complex amplitudes {

*B*

_{m}(

*z*= 0, Ω)}. The total input amplitude in mode

*m*can be written as

*m*at any point

*z*can be decomposed into the signal (Ω = 0) and noise (Ω ≠ 0) amplitudes,

*g*. Below the TMI threshold, signal amplitudes roughly grow independently of the noise and in the case of no gain saturation are given by

*m*. More importantly, noise also grows due to the transfer of power from the signal due to the thermo-optical scattering. The noise growth can be obtained by solving Eq. (11), which is first simplified by retaining only the phase-matched terms, given by the condition

*β*

_{m}−

*β*

_{n}+

*β*

_{i}−

*β*

_{j}= 0, leading to

*m*=

*j*,

*i*=

*n*as a solution to the phase-matching condition, reducing the triple sum in the second term of Eq. (11) to a single sum in Eq. (16). Note that

*m*=

*n*,

*i*=

*j*is also a valid solution to the phase-matching condition, but such terms correspond to only a cross-phase modulation, but no transfer of power as these terms result in only a uniform temperature shift (

*q*

_{ij}≈ 0); thus, we do not retain these terms. We have also assumed the absence of any exact degeneracies. This assumption can break down especially in graded index fibers where significant degeneracies are present in which case the above equations need to be modified slightly. For a detailed discussion, see the supplementary material, Sec. I. Next, we substitute the ansatz for

*A*

_{m}from Eqs. (14) and (15) in Eq. (16) and simplify the RHS to obtain the growth equation for

*B*

_{m}(

*z*, Ω),

*m*and is given by $|Ams|2$, and

*G*

_{mnnm}represents the contribution to noise growth in mode

*m*by signal power in mode

*n*. Note that here we have considered a signal at a single frequency without any linewidth broadening since our goal is to study TMI in narrowband fiber amplifiers with a focus on coherent excitations. However, in some applications of fiber amplifiers, signal linewidth broadening is utilized. In that case, our theory can be straightforwardly generalized by directly solving Eq. (16) instead of simplifying it to Eq. (17). We can convert noise amplitude growth equations for mode

*m*into growth equations for noise power (given by |

*B*

_{m}|

^{2}) by multiplying with $Bm*$ on both sides and adding the complex conjugate term,

*m*. The first term on the RHS in Eq. (17) corresponds to the linear optical gain, and each term in the second term corresponds to a transfer of power from a particular signal mode

*n*to noise in mode

*m*due to the thermo-optical scattering. We have ignored the self-coupling term

^{18}(corresponding to

*G*

_{mmmm}) since the grating formed by interference of mode

*m*with itself at a Stokes shifted frequency has a grating period [2

*π*/(

*q*

_{mm}≈ 0)] much larger than the length of the fiber. For each mode pair (

*m*,

*n*), we have defined a thermo-optical coupling coefficient

*χ*

_{mn}(Ω), which is equal to −2 times the imaginary part of

*G*

_{mnnm}multiplied with

*χ*

_{0}. Equation (18) can be solved analytically, resulting in exponential growth in noise power in each mode

*m*,

*m*at the input end of the fiber. $PmN(L,\Omega )$ is the noise power at the output upon exponential growth both due to the linear gain

*g*and the TMI gain, which depends on signal power distribution,

*z*integral can be simplified by using the signal growth formula given in Eq. (15), $\u222b0LgPns(z)dz=Pns(L)\u2212Pns(0)\u2248\Delta PP\u0303ns$, where Δ

*P*is the total signal power extracted from the amplifier and $P\u0303ns$ is the fraction of the signal power in mode

*n*$(\u2211nP\u0303ns=1)$, which is determined by the input excitation. Thus, TMI gain in mode

*m*is given by

^{18,19}where it is given by the product of extracted signal power and the thermo-optical coupling between the two modes

*χ*

_{21}. More generally, Eq. (21) can be used to derive the formula for the TMI threshold under arbitrary multimode excitations, as shown in Sec. II C.

### C. Multimode TMI threshold

*ξ*> 1%) of the signal power, which leads to the onset of a fluctuating beam profile.

^{18,19}For a given amount of input noise power, this occurs when the highest TMI gain across all frequencies and modes becomes large enough. Therefore, a quantitative condition for the threshold can be derived from Eqs. (19) and (21),

*P*

_{th}is the TMI threshold,

*L*is the length of the fiber, and

*P*

^{s}(0) and

*P*

^{N}(0) are the input signal and noise powers, respectively. We have introduced an overall thermo-optical coupling coefficient $\chi \u0304$, which is equal to the maximum value of the effective thermo-optical coupling coefficients across all the modes and Stokes frequencies,

*χ*

_{mn}and the input power distribution, $Pns\u0303$. It also weakly depends (logarithmically) on the input noise power and the noise fraction

*ξ*at which the threshold is set. In addition, any increase in the input signal power (or, seed power) leads to a corresponding increase in the TMI threshold. Typically, the input power is significantly smaller than the TMI threshold; thus, the overall relative change due to a seed power increase is small.

^{70,71}Note that in Eq. (22), we have approximated the total Stokes power by the Stokes power in the mode with maximum Stokes gain at Stokes frequency with peak gain. This is justified by the exponential nature of the Stokes growth. For a detailed derivation, see the supplementary material, Sec. II. It can be easily verified that when only the fundamental mode is excited $(P\u03031=1,P\u0303n\u22601=0)$, the overall thermo-optical coupling coefficient, $\chi \u0304$ simplifies to

*χ*

_{21}, reducing our multimode TMI threshold formula to a previously derived formula in studies that only consider single mode excitation.

^{19}More generally, multimode excitation provides a parameter space to control the overall thermo-optical coupling coefficient and the TMI threshold, even for a fixed fiber. The TMI threshold is highest for the signal power distribution, which minimizes the maximum effective thermo-optical coupling coefficient. The amount of tunability therefore depends on the properties of the pairwise thermo-optical coupling coefficients, which we discuss in detail in Sec. II D.

### D. Thermo-optical coupling

*m*and

*n*is proportional to the imaginary part of a particular component of Green’s function,

*G*

_{mnnm}, and is equal to

*χ*

_{mn}has contributions from each temperature eigenmode

*k*. The strength of each contribution is proportional to the integrated overlap of the dot product of the optical mode profiles $\psi \u20d7m$ and $\psi \u20d7n$ with the temperature mode profile $T\u0303k$. Each contribution has a characteristic frequency curve, which peaks at a frequency given by inverse mode lifetime Γ

_{k}, and the peak value is proportional to the thermal mode lifetime $4\pi \Gamma k\u22121$. To illustrate the key properties of the thermo-optical coupling coefficient, we calculate

*χ*

_{mn}(Ω) for all the mode pairs of a circular step-index fiber with a core diameter of 40

*μ*m, a cladding diameter of 200

*μ*m, and a numerical aperture NA = 0.15, supporting 82 modes per polarization. Detailed fiber or optical parameters used for calculation are provided in Table I. We consider optical modes linearly polarized (LP) along the

*x*-direction. Each LP mode is designated by two indices [

*u*,

*v*], which correspond to the number of azimuthal and radial nodes, respectively

^{73,78}[see Fig. 2(a)]. The modes are arranged in the order of decreasing longitudinal propagation constants. For instance, first five modes (excluding rotations) are LP01 (FM), LP11, LP21, LP02, and LP12. In our notation, therefore,

*χ*

_{12}≡

*χ*

_{LP01,LP11}. Note that LP modes are only approximate modes of the fiber, which become inaccurate for fibers with long lengths, in which case the exact vector fiber modes (HE and EH modes) need to be considered.

^{73}For the purpose of our discussion, we assume that LP modes are accurate enough for the length of the fiber considered.

Parameter . | Value . |
---|---|

Core shape | Circular |

Core diameter (μm) | 40 |

Cladding diameter (μm) | 200 |

Core refractive index | 1.458 |

Cladding refractive index | 1.45 |

Signal wavelength, λ_{s} (nm) | 1032 |

Pump wavelength, λ_{p} (nm) | 976 |

Number of optical modes | 82 |

Thermo-optic coefficient, η (K^{−1}) | 3.5 × 10^{−5} |

Thermal conductivity, κ (W m^{−1} K^{−1}) | 1.38 |

Diffusion constant D (m^{2} s^{−1}) | 8.46 × 10^{−7} |

Parameter . | Value . |
---|---|

Core shape | Circular |

Core diameter (μm) | 40 |

Cladding diameter (μm) | 200 |

Core refractive index | 1.458 |

Cladding refractive index | 1.45 |

Signal wavelength, λ_{s} (nm) | 1032 |

Pump wavelength, λ_{p} (nm) | 976 |

Number of optical modes | 82 |

Thermo-optic coefficient, η (K^{−1}) | 3.5 × 10^{−5} |

Thermal conductivity, κ (W m^{−1} K^{−1}) | 1.38 |

Diffusion constant D (m^{2} s^{−1}) | 8.46 × 10^{−7} |

In Fig. 3(a), we have shown the thermo-optical coupling coefficient between the FM (LP01) and three HOMs in increasing mode order (LP11, LP02, and LP12) for this fiber. Relevant mode profiles are shown in the inset. The thermo-optical coupling coefficient in all three cases is zero at Ω = 0 and attains its peak value at frequencies on the order of few kHz, consistent with the millisecond timescale of TMI.^{11,18,19} The FM couples most strongly with the lowest HOM (blue curve), and the peak coupling decreases with the increasing transverse spatial frequency mismatch between the modes (green and orange curves). The peak frequency and linewidth are higher for the larger spatial frequency mismatch between the modes. These results can be understood by investigating the specific temperature eigenmodes, facilitating the coupling for different mode pairs, as shown in Fig. 3(b). The contribution of each temperature eigenmode to the coupling between a pair of optical modes is proportional to (1) the thermo-optical modal overlap integral and (2) the corresponding thermal mode lifetime $4\pi \Gamma k\u22121$ [Eq. (25)]. Overall, only a few temperature eigenmodes, which match the transverse intensity profile of the interference between the optical modes, have a significant overlap integral. For LP01−LP11 coupling, the optical interference pattern has relatively lower spatial frequencies and thus has significant overlap with relatively lower-order temperature eigenmodes [shown by blue dotted bars in Fig. 3(b)]. On the other hand, the LP01–LP12 mode pair has a significant overlap with relatively higher-order temperature eigenmodes [shown by orange bars in Fig. 3(b)]. The thermal mode lifetime is shorter for temperature eigenmodes with higher transverse spatial frequency [as shown by the red curve in Fig. 3(b)], so the contribution of higher-order temperature eigenmodes is damped. Such intrinsic dampening of high spatial frequency contributions to the temperature is characteristic of the heat propagation being a diffusive process.^{79} Consequently, the lifetime weighted thermo-optical mode overlap for neighboring optical modes (blue solid bars) is much higher than for modes with larger separation of transverse spatial frequency (orange solid bars). As a result, the peak thermo-optical coupling coefficient decreases with the increasing spatial frequency mismatch between the optical modes. The quantitative relation between the spatial frequency mismatch and the lowered thermo-optical coupling coefficient is made possible by the use of thermal eigenmodes (each with a definite spatial frequency) for expressing the temperature fluctuations.

Weak thermo-optical coupling between modes with large spatial frquency mismatch is a generic result for all optical mode-pairs. To explicitly show this, we calculate the peak value of *χ*_{mn}(Ω) for all $\u223c104$ mode pairs in circular fiber. The resulting coupling matrix is shown in Fig. 4(a) as a false color image. We have omitted the self-coupling (*m* = *n*) since it is not responsible for intermodal power transfer between the modes.^{18} Clearly, every mode couples strongly with only a few modes, resulting in a highly banded/sparse coupling matrix. The strong coupling occurs when the spatial frequencies of two modes are closely matched, i.e., similar number of radial and azimuthal nodes (Δ*u*, Δ*v* ≤ 1). As we will see below, this banded nature of the coupling matrix is what allows for a significant suppression of TMI upon multimode excitation.

### E. Threshold scaling

Recall that the TMI threshold is defined as the output signal power at which the noise power in any mode becomes a significant fraction $(>1%)$ of the output signal power,^{19} and the beam profile fluctuates dynamically rendering it useless for many applications.^{11} According to Eqs. (23) and (24), the TMI threshold is inversely proportional to the overall thermo-optical coupling coefficient $\chi \u0304$, which is a sum of *χ*_{mn} weighted by the fraction of signal power in each mode, which depends on the input excitation. Most previous TMI suppression efforts consider exciting only the fundamental mode and avoid sending power in HOMs, which our results indicate is not the optimal approach. To show this explicitly, we consider two special cases of input excitation—(1) FM-only excitation: all of the signal power is present in the FM $(P\u03031s=1)$ and noise is present in HOMs, and (2) equal mode excitation: the input power is divided equally in *M* modes of the fiber $(P\u0303ns=1/M)$. We use Eq. (24) to calculate $\chi \u0304$ in the two cases and use it to compare the TMI threshold. In the first case, the noise in the first HOM (*m* = 2) has the highest growth rate due to the signal power in the FM (*m* = 1), giving $\chi \u0304=\chi 21$. Since *χ*_{21} is the highest coupling out of all mode pairs [Fig. 4(a)], FM-only excitation actually leads to the lowest TMI threshold. In the second case, all the modes contribute to the overall coupling to a given mode (say, *m* = 2) but with weights reduced by a factor of *M*, i.e., $\chi \u0304=\u2211n\u22602\chi 2n/M$. Due to the banded nature of *χ*_{mn}, only a few elements in the sum are significant, leading to $\chi \u0304\u2248s\chi 21/M$. Here, *s* is used to denote the average number of significant elements in any row of the thermo-optical coupling matrix, which does not scale with *M* and is roughly equal to 6 for our circular fiber. Thus, for a highly multimode excitation (*M* ≫ *s*), the TMI threshold can be significantly higher than the FM-only excitation. In fact, our reasoning predicts that the TMI threshold will increase linearly with the number of equally excited modes, *M*. To verify this explicitly, we calculate the TMI threshold as *M* is increased. We define the threshold increase factor (TIF) as the ratio of TMI threshold for *M*-mode excitation with that for FM-only (*M* ≈ 1) excitation. Figure 4(b) shows the values of TIF as *M* is varied. As predicted, the TMI threshold increases linearly with the number of excited modes with a slope roughly equal to 1/*s* (≈0.16). When all 82 modes in the circular fiber are excited equally, we find a remarkable 13× enhancement of the TMI threshold over FM-only excitation. This confirms that highly multimode excitation in a multimode fiber amplifier can be a great approach for achieving robust TMI suppression.

Our model and formalism are equally efficient for considering fibers with non-circular cross sections, such as the D-shaped fiber with more complex spatial structures of its modes, due to the underlying chaotic ray dynamics, as discussed and depicted in Fig. 7. The physical argument for the suppression of TMI via highly multimode excitation does not rely on any special features of the fiber transverse modal geometry and should be valid for general step-index fiber geometries. Confirmation of this expectation is shown in the Appendix, where we find linear scaling of the TMI threshold with the number of modes excited for the D-shaped fiber. The slope of the scaling line does depend on both the geometry of the core and the cladding, but rather weakly, as discussed in the Appendix.

It should be noted that although we have only considered equal-mode excitation to achieve a higher TMI threshold, it is by no means a strict condition. Our theory predicts that multimode excitation will generically lead to a higher threshold than the FM-only excitation due to the sparse nature of the thermo-optical coupling matrix. Therefore, it is expected to be a universal result, independent of the details of fiber composition, geometry (see above) and the precise distribution of excited modes. The level of enhancement will mainly depend on the effective number of excited modes. While most previous studies of TMI consider FM-only excitation, an increase in TMI threshold due to multimode excitation has been observed in some cases. A recent study by the current authors demonstrated an increase in TMI threshold upon multimode excitation using explicit time-domain numerical simulations of coupled optical and heat equations in 1D cross-sectional waveguides.^{63} Additionally, several recent experimental studies on the effect of fiber bending on TMI threshold in few-mode fibers provide evidence for our predictions.^{80–82} They have observed that upon launching light in a fiber that supported multiple modes, a higher TMI threshold is obtained for a larger bend diameter (i.e., a loosely coiled fiber). The bending loss of HOMs decreases when the bend diameter is larger, and thus, the effective number of excited modes is higher, which leads to a higher TMI threshold in accordance with our predictions. It should be noted that these experimental findings are opposite to the predictions of previous theories, which only consider FM excitation and predict that decreased HOM loss upon increasing bend diameter should lead to lowering of the TMI threshold.^{44,83} The authors of these studies^{80,81} pointed out this anomaly and tried explaining these results by arguing that lower mode mixing leads to a higher threshold. In our opinion, this discrepancy is a result of ignoring the signal power in HOMs, which is fully taken into account in our formalism and can straightforwardly account for their results.

## III. NONLINEAR MODEL: GAIN SATURATION AND PUMP DEPLETION

In Sec. II, we approximated the signal growth by considering a simplified treatment of optical gain due to stimulated emission by assuming a constant gain coefficient *g*_{0}. This model was treated first to highlight the physics of multimode excitation and its impact on the thermo-optical coupling and TMI threshold. However, such a model of fiber gain neglects important effects, such as gain saturation,^{65,84} mode-dependent gain/loss,^{85} and pump depletion,^{71} which are present in any real high power fiber amplifier. In this section, we utilize a more realistic model of gain in multimode fiber amplifiers, including all the above-mentioned effects, and again are able to derive a more complicated, but still semi-analytical formula for thermo-optical coupling, which gives the TMI threshold for arbitrary multimode excitations. This model should be used in order to make quantitative predictions to be compared with experiments.

We will show that the banded nature of the thermo-optical coupling matrix and the linear scaling of TMI threshold upon multimode excitation found for the simpler model in Sec. II are also found in the presence of gain saturation, mode-dependent gain, and pump depletion. This was expected, as the scaling of TMI threshold upon multimode excitation is a result of the diffusive nature of heat propagation, leading to weak thermo-optical coupling between modes with a large transverse spatial frequency mismatch. This mismatch of the spatial frequencies of the thermal and optical fluctuations is qualitatively unchanged by gain saturation and the other effects. We do find that including gain saturation impacts the exact value of TMI threshold, increasing it for all input excitations due to reduced dynamic heat load; this result is in agreement with the findings of previous studies of TMI in the case of single-mode excitations, which include these effects.^{21,26,65,66,70,71,82,86}

To set up our new model, we generalize the approach used by several previous studies to obtain TMI threshold upon single mode excitation^{21,26,64,66,70,71,82,86} for the case of multimode excitation. This involves two key steps. First, the signal amplitudes are obtained in each mode and along the entire fiber axis by numerically solving the saturated signal amplification equations coupled with the evolution of pump and the upper-level population in the gain medium. In this step, any effect of Stokes wave (due to noise) and the thermo-optical coupling on the signal are neglected. This is consistent with the undepleted signal approximation utilized in Sec. II B since as argued there, the backaction of the noise on the signal is negligible below the TMI threshold. In the second step, the signal amplitudes and upper-level population are used to calculate relevant dynamic heat load and the resultant thermo-optical coupling is used to obtain the growth of Stokes power in each mode, which determines the TMI threshold.

### A. Signal amplification

^{87}in the fiber core, with a pump beam propagating in the pump core exciting the electrons to the upper level and creating inversion. As a result, the signal beam propagating in the fiber core undergoes amplification due to stimulated emission. The growth of signal amplitudes in each mode, $Ams$, can be written as

^{13}

^{,}

*g*

_{mn}is the gain in mode

*m*due to amplitude in mode

*n*and is given by

^{84,87}

*g*

_{0}is the unsaturated gain coefficient. The amount of gain saturation depends on the ratio of the local signal intensity

*I*

_{0}and a saturation intensity

*I*

_{sat}. Both

*g*

_{0}and

*I*

_{sat}depend on various properties of the gain medium and the local value of pump power

*P*

_{p}, and their formulas are obtained from steady state rate equations

^{13,88}(for more details, see the supplementary material, Sec. III). Notice that in addition to self-gain terms (

*m*=

*n*), we must also consider cross gain terms (

*m*≠

*n*), which are now non-zero due to the spatially varying nature of the gain saturation term (referred to as spatial hole burning) and leads to non-linear mode coupling due to gain saturation.

^{13,84,88}Some previous studies of TMI, which consider irradiance based models,

^{66,86}ignore the spatial hole burning retaining only the self-gain terms. For studying multimode excitations, including the hole burning is essential since the multimode signal intensity will typically have rapid speckle-like spatial variations. Note that we have ignored any variation in the real part of the refractive index due to the gain saturation as this tends to be significantly smaller than the imaginary part for relevant wavelengths.

^{13}We take this into account quantitatively by introducing a standard equation for the evolution of pump power, which is given by

^{13}

^{,}

*g*

_{p}, which depends on the local inversion and is given by

^{13}

^{,}

*N*

_{Yb}is the density of doped ytterbium atoms, and

*A*

_{cl}is the area of pump core, which is typically the same as the core plus the first cladding. Note that we are able to consider a single equation for pump power instead of a different one for each pump mode as in typical experiments the pump is incoherent and fills the fiber cross section uniformly on average.

^{11,69,87}The Yb-doped gain medium we consider is a quasi-three-level medium. The pump photon causes a transition from the lower lasing level to the uppermost level after which the electron population immediately transitions downward to the upper lasing level causing heat generation in the process. This is known as quantum defect heating.

*n*

_{u}and

*n*

_{l}are the fraction of population in the upper and lower lasing levels, respectively, satisfying

*n*

_{l}= 1 −

*n*

_{u}, where

*n*

_{u}is obtained by solving steady-state rate equations, leading to

^{13,26,70,71,87}

Here, *ω*_{s} and *ω*_{p} are the signal and pump frequencies and $\sigma se$ and $\sigma sa$ are the emission and absorption cross sections of the signal, respectively. Each individual term involving scattering cross sections is equal to either the rate of stimulated emission or absorption of the signal or the pump and has units of [s^{−1}]. *τ* is the upper state lifetime. *n*_{u} has a strong spatial dependence due to both the pump power *P*_{p} and the signal intensity *I*_{0}.

A solution to the signal amplification equations [Eqs. (26) and (27)] coupled with the evolution of the pump power [Eqs. (28) and (29)] and upper-level population [Eq. (30)] can be obtained by standard finite difference based numerical methods.^{89,90} We utilize a centered difference approximation for all the *z* derivatives and use it to iteratively update the values of the signal amplitudes and the pump power at the *z* + *δz* position from the values of respective gain coefficients at *z* and values of the variables at *z* − *δz* positions. At each *z*, the values of *g*_{p}, *g*_{mn}, and *n*_{u} are updated directly from Eqs. (27), (29), and (30). To validate our numerical model, we first simulated the 20/400 fiber amplifier studied by Smith and Smith^{64} where single mode excitation was considered and obtained excellent agreement with their results. We also simulated the two-mode excitation case discussed by Li *et al.*^{82} and found good agreement.

Next, we study the highly multimode step-index fiber amplifier discussed in Sec. II D. All the parameters for Yb-doped gain medium are taken from Table I in the work of Smith and Smith.^{64} The radius of the Yb doped area was considered to be equal to the fiber core radius. The unsaturated gain coefficient *g*_{0} ≈ 4 m^{−1}, dopant concentration *N*_{Yb} = 3.25 × 10^{25} m^{−3}, and length of the fiber *L* = 2 m. Initial pump power depends on the TMI threshold for different excitations. For any given multimode input excitation, we launch 10 W of seed power divided appropriately in the signal modes and a variable amount of pump power determined self-consistently such that the total output signal power is equal to the TMI threshold. The length of the fiber is chosen to be such that $>95%$ pump power is absorbed and converted to the signal power. Figure 5 shows the results for the case when all 82 modes are excited equally at the input. 6 kW of pump power is launched in a 2 m long fiber, and it decreases monotonically along the fiber axis as it excites the gain medium creating inversion. The signal power in each mode grows at first exponentially when the gain saturation is low and linearly when gain saturation becomes high and eventually flattening out when nearly all the pump is depleted. The final output power is roughly 5675 W, leading to a conversion efficiency of 0.945, which is close to the quantum efficiency limit $(\lambda p\lambda s=0.9457)$. Notice that the signal power in different modes grows with a different growth rate displaying mode dependent gain. This is a result of both the differential overlap with the fiber core where the gain medium is present and spatial hole burning induced variations in the fiber gain. At this stage, thermo-optical effects have not yet been calculated. In Subsection III B, having obtained the saturated signal power along the entire fiber axis, we show how the gain-saturated thermo-optical coupling can be utilized to calculate the TMI threshold for arbitrary multimode excitations.

### B. Thermo-optical coupling with gain saturation

^{66,71}

*g*

_{s}and

*g*

_{0}are the saturated and unsaturated gain coefficients and

*q*

_{D}is the quantum defect. The denominator comes directly from the saturation term in the gain coefficient and depends on the local signal intensity $I(r\u20d7,t)$. At any point, the signal intensity can be written as a sum of a static contribution $I0(r\u20d7)$ and a dynamic contribution $I\u0303(r\u20d7,t)$, i.e., $I(r\u20d7,t)=I0(r\u20d7)+I\u0303(r\u20d7,t)$. The first term is a result of self-interference of the electric field at the signal frequency

*ω*

_{0}, whereas the second term results from the interference between the signal at

*ω*

_{0}and the noise at Stokes shifted frequencies

*ω*

_{0}− Ω. Consequently, the heat profile has both static and dynamic contributions, $Q(r\u20d7,t)=Q0(r\u20d7)+Q\u0303(r\u20d7,t)$. It is the latter term that results in the dynamic power transfer responsible for TMI. Below the TMI threshold, the Stokes power is much smaller than the signal power, and hence, the dynamic part of the intensity fluctuation is significantly smaller than the static part. Therefore, we can simplify the expression for $Q\u0303$ by using the Taylor expansion in terms of $I\u0303/(I0+Isat)$ and keeping the leading-order term, giving the following expression

^{66}(for a detailed derivation, see the supplementary material, Sec. III):

The dynamic heat source is proportional to the time varying intensity $I\u0303$ and is inversely proportional to the square of the saturation term, which depends on local static intensity *I*_{0}. This quadratic behavior of the saturation term is a result of the extra contribution to the gain saturation from the dynamic part of the intensity, which comes with a negative sign upon the Taylor expansion, and partially cancels the direct term proportional to $I\u0303$.^{66} Such a quadratic saturation of the dynamic heat load explains the well-known increase in the TMI threshold due to the gain saturation^{13,82} and was first pointed out by Hansen and Lægsgaard^{66} and later by several other authors.^{69,71,86}

*I*

_{0}in the denominator depends on

*z*, now, the thermo-optical coupling varies along the fiber axis. This modifies the form of the TMI gain given in Eq. (20) by shifting $\chi mns$ inside the

*z*integral,

*n*calculated with gain saturation in the first step discussed in Subsection III A. Note that in the above formula, the unsaturated gain coefficient

*g*

_{0}is used, as the effect of saturation on the heat load is already considered in the thermo-optical coupling $\chi mns$. For a direct comparison with the form of the TMI gain given in Eq. (21), we define a

*z*-integrated thermo-optical coupling as follows:

*z*-integrated

*χ*becomes equal to the value of

*χ*everywhere. This allows us to write an expression for TMI gain similar to the one in Eq. (21),

The TMI gain in any mode is given by a product of the total extracted power Δ*P* and the weighted sum of *z*-integrated thermo-optical coupling with all the other modes. The weight $P\u0303ns$ is the fraction of signal power extracted by mode *n* $(P\u0303ns=Pns(L)\u2212Pns(0)\Delta P)$ and can be controlled by the input excitation.

Using the formulas in Eqs. (33) and (35), we were able to calculate the thermo-optical coupling for all the mode pairs in the 82-mode circular step-index fiber amplifier discussed in Sec. II everywhere along the fiber when all 82-modes were equally excited and calculated the *z*-integrated coupling matrix in the presence of gain saturation. The results are shown in Fig. 6(a). It can be seen that the thermo-optical coupling matrix remains sparse and approximately banded, even with gain saturation. This is as we expected from the physical arguments above. Although the structure of the matrix is unchanged, the exact values of the couplings are certainly impacted by the gain saturation. It can be seen by comparing the scale in Figs. 3(a) and 6(a) that the *z*-integrated thermo-optical coupling is lowered due to the gain saturation for all the mode pairs by roughly 20%. This leads to a higher value of TMI threshold upon including the gain saturation as we will see in Subsection III C. The amount of coupling reduction depends on the ratio of the size of the pump core and the gain core as that determines the degree of gain saturation. For a larger pump core, the pump intensity in the gain core is lower, decreasing the saturation intensity and increasing the amount of gain saturation, leading to a smaller value of z-integrated thermo-optical coupling. Note that the reduction in overall thermo-optical coupling also depends on the number of modes excited albeit weakly.

### C. TMI threshold scaling

^{21,26,64,66,70,71,82,86}To calculate the TMI threshold for multimode excitations, we utilize the formula for TMI gain in Eq. (36), which produces a similar formula for TMI threshold as given in Eq. (23). The only modification is that the overall thermo-optical coupling coefficient $\chi \u0304$ is now given by a weighted sum of

*z*-integrated thermo-optical coupling between various mode pairs,

The TMI threshold is inversely proportional to the overall thermo-optical coupling coefficient $\chi \u0304$, which strongly depends on the input excitation via signal power distribution in various modes ${P\u0303ns}$. As in the simpler model, due to the sparse nature of the saturated thermo-optical coupling matrix, we expect power division to increase the TMI threshold. To show this quantitatively, we first calculate the TMI threshold when only the fundamental mode (FM) is excited and subsequently when various number of modes are equally excited. In each instance, the amount of pump power is chosen such that the total output signal power is equal to the TMI threshold. When only the FM is excited, we obtain a TMI threshold equal to 775 W. To validate this result, we compare it to the predicted value of the TMI threshold in the work of Smith and Smith^{65} for a fiber with the same ratio of the diameter of the gain core and pump core (50/250) as in our fiber (40/200). A value of 785 W is obtained for the TMI threshold for such a fiber from Table II in Ref. 65, which matches closely with our prediction. Without the gain saturation, the TMI threshold is found to be 375 W, which is also in close agreement with prediction in Ref. 65 (345 W). In the previous studies on narrowband amplifiers involving large mode area (LMA) fibers (both step-index and photonic-crystal fibers), the TMI threshold ranged from 200 to 800 W.^{42,91,92} The TMI threshold upon FM-only excitation in this fiber is therefore comparable to LMA fibers. As the number of modes excited is increased, the TMI threshold also increases linearly especially for a large number of modes and reaches close to 5700 W when all 82 modes are equally excited [Fig. 6(b)]. The origin of this linear scaling is the same as in the simpler model. The slope of the linear increase in TMI threshold is higher with saturation (62 W/mode) than without gain saturation (54 W/mode) since the starting point of the curve (fundamental mode threshold) is higher when gain saturation is included. The value of TMI threshold with gain saturation is higher by an amount ranging from 100 to 800 W for various numbers of excited modes. For a fiber with a larger pump core or a smaller gain core, this difference is expected to be even higher. As mentioned above, the reduction in thermo-optical coupling due to gain saturation that leads to increased threshold weakly depends on the number of modes excited. When gain saturation is taken into account, this effect can compete with linear enhancement in the TMI threshold for a low number of modes causing deviations from strict linear scaling at the lower end in Fig. 6(b). For a large number of excited modes, linear enhancement dominates giving asymptotically linear scaling of TMI threshold with the number of modes excited.

## IV. DISCUSSION AND CONCLUSION

In this work, we have developed a theory of TMI, which can be used for efficient calculation of TMI threshold in narrowband fiber amplifiers for arbitrary multimode excitations and fiber geometries. A key result obtained from this theory is that TMI threshold increases linearly with an effective number of excited modes. The scaling is a result of the diffusive nature of the heat propagation underlying the thermo-optical mode coupling. As such, it is expected to be applicable to a broad range of fibers. This opens up the use of highly multimode fibers as a promising avenue for instability-free power scaling in high-power fiber amplifiers.

The existing high-power multimode fiber amplifiers do not produce a spatially coherent output that can be focused to a diffraction-limited spot or collimated to a Gaussian beam. By contrast, our approach allows for the generation of a spatially coherent beam out of a highly multimode fiber amplifier. This capacity will expand the potential applications of multimode fiber amplifiers, e.g., for a large-scale laser interferometer and coherent beam combining. It is sometimes mistakenly assumed that the presence of multiple spatial modes necessarily leads to a poor output beam quality and a high value of *M*^{2}.^{49,50} However, this belief is not correct, as previous work on the manipulation of coherent optical fields has been shown.^{51,52,93} Indeed, the fields in multimode fibers can maintain a high output beam quality as long as the light remains sufficiently coherent,^{94,95} i.e., the signal linewidth is narrower than the spectral correlation width of the fiber. In a typical scenario involving a 10-m-long silica fiber with an NA of 0.1, the spectral correlation width is ∼2 GHz for a laser wavelength of 1032 nm. In fact, by using an SLM to wavefront-shape the input light to a fiber, it is possible to obtain a diffraction-limited spot after coherent multimode propagation in both passive^{54} and active fibers,^{51} which can be easily collimated using a lens. It is noteworthy that the SLM may introduce a few percent of optical loss due to diffraction by pixel edges into higher orders and absorption by the various layers in the SLM. However, since the SLM modulates only the input seed, its loss is negligible, as long as the amplifier operates in the gain-saturation regime where the output power mostly depends on the pump power. Since focusing to a diffraction-limited spot in the near field necessarily leads to the excitation of multiple fiber modes, our theory predicts that it will lead to a higher TMI threshold than FM-only excitation, with a scaling proportional to the effective number of fiber modes. Hence, the method just discussed can be used to increase the TMI threshold while maintaining good beam quality. A similar result has recently been demonstrated experimentally for the case of SBS in multimode fibers; the threshold for the onset of SBS was demonstrated to increase significantly by focusing the output light of a multimode fiber to a diffraction-limited spot.^{61}

Above, we derived computationally tractable semi-analytic formulas for the TMI threshold in the presence of multimode excitation: First with a model that neglects gain saturation and related effects but illustrates the fundamental mechanism of suppression of TMI through the banded nature of the thermo-optical coupling matrix, leading to a linear threshold increase with the number of excited modes in the signal. Within this model, all calculations are effectively linearized. Second, we improved the model by calculating the saturated signal and pump depletion with a non-linear computational method, which is then used in a generalization of the first model to calculate the thermo-optical coupling matrix and the TMI threshold. The model again finds a sparse coupling matrix and a linear scaling of the threshold with the number of excited modes. The absolute threshold is increased compared to the simpler model, as expected, but the difference is not dramatic. This model incorporates all of the major effects missing in the first model but present in realistic fiber experiments: gain saturation and hole-burning of the signal, mode-dependent gain/loss, and depletion of the pump; hence, it is appropriate for quantitative comparisons with experimental data.

We have still ignored any random linear mode coupling^{85} in the fiber. This assumption can be relaxed relatively straightforwardly; the multimode TMI threshold is typically increased by the presence of mode mixing since it promotes equipartition of signal power in various modes.^{63} Note that our results are valid only under the assumption that the random linear mode coupling does not have a strong temporal variation on the time scales faster than what can be compensated by the use of spatial light modulators. This is typically true for most high-quality multimode fibers.

Experimental validation of the theoretical results provided in this paper would be an important next step. Time-domain numerical simulations of optical and heat equations for up to five mode excitations in 1D cross-section waveguides are provided in Ref. 63 and are in good agreement with our theory. Additionally, several recent experimental studies investigating the effect of fiber bending on the TMI threshold in few-mode fibers provide evidence for our predictions.^{80–82} It has been observed that in a few-mode fiber, when the bend diameter of the fiber is increased (fiber being loosely coiled), a higher TMI threshold is obtained. For a large bend diameter, the bending loss of the HOMs decreases; thus, the effective number of excited modes is higher, which leads to a higher TMI threshold in accordance with our predictions. In fact, the previous theories that only consider FM-only excitation predict that decreasing the HOM loss by increasing bend diameter should lead to lowering of the TMI threshold, in contrast to the experimental findings.^{44,83} This discrepancy is a result of ignoring the signal power in HOMs, which is fully taken into account in our theory. More systematic experimental studies are needed to investigate our predictions quantitatively.

In recent years, there has been a significant progress in fabricating low-loss fibers with non-circular cross sections.^{96,97} These fibers have been proposed as a way to manipulate the strength of nonlinearities.^{98} As noted above, the theory derived here can be used to calculate the TMI threshold for any fiber cross-sectional geometry. In the Appendix, we have utilized it to demonstrate that a D-shaped fiber performs better than a standard circular fiber in raising the TMI threshold via multimode excitation. This interplay between input excitation and fiber geometry can provide a number of avenues for manipulating the strength of nonlinearities in the fiber. We believe that our theory can be utilized as a framework for customizing the strength of the thermo-optical nonlinearity using optimized fiber designs.^{99} Utilizing the spatial degrees of freedom of the fields to control nonlinear effects is becoming a standard tool and has been demonstrated for diverse effects, such as SBS,^{59–61} SRS,^{62} and Kerr nonlinearity.^{56,100} Our work contributes to this exciting area by bringing the thermo-optical nonlinearity, with its different physical origin, into this category and at the same time providing a solution to the practical challenge of TMI suppression.

## SUPPLEMENTARY MATERIAL

A supplementary material is attached for additional information. In the first section, we provide a detailed discussion on the phase matching for nonlinear thermo-optical scattering in multimode fibers. In the second section, we provide more details regarding the derivation of the multimode TMI threshold formula and justify the use of only the Stokes mode with maximum gain at peak frequency. In Sec. III, we provide a derivation of the saturated gain coefficient along with the formulas for various gain saturation parameters. In addition, we derive the formula for saturated dynamic heat load. Finally, a quantitative comparison of our model is provided with previous studies on TMI in single mode and few mode fiber amplifiers.

## ACKNOWLEDGMENTS

We thank Ori Henderson-Sapir, Heike Ebendorff-Heidepriem, and David Ottaway at The University of Adelaide; Stephen Warren-Smith and Linh Viet Nguyen at the University of South Australia; and Peyman Ahmadi at Coherent for stimulating discussions. We acknowledge the computational resources provided by the Yale High Performance Computing Cluster (Yale HPC). This work was supported by the Air Force Office of Scientific Research (AFOSR) under Grant No. FA9550-20-1-0129. We also acknowledge the support of Simons Collaboration on Extreme Wave Phenomena Based on Symmetries.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Kabish Wisal**: Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Chun-Wei Chen**: Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Validation (supporting); Visualization (supporting); Writing – review & editing (equal). **Hui Cao**: Conceptualization (lead); Formal analysis (supporting); Funding acquisition (lead); Investigation (supporting); Project administration (equal); Supervision (supporting); Writing – review & editing (equal). **A. Douglas Stone**: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (lead); Investigation (supporting); Project administration (equal); Supervision (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: IMPACT OF FIBER GEOMETRY

In the main text, we considered a standard step-index multimode fiber with a circular cross section, which shows that TMI threshold increases linearly with the number of excited modes owing to the banded nature of the thermo-optical coupling matrix. This banded nature is a fundamental consequence of the diffusive nature of the heat propagation, which underlies the thermo-optical scattering due to the intrinsic damping of high-spatial-frequency features. As a result, we expect the banded nature of the thermo-optical coupling matrix to be maintained even in fibers with non-standard geometries,^{96,97} while the number of significant elements in the coupling matrix and their relative values can depend on the particular fiber geometry. As such, we expect the linear scaling of the TMI threshold to be maintained, but the slope to differ for fibers with different cross-sectional geometries. Thus, studying different fiber geometries to determine which yields a higher slope of threshold increase with *M* may be useful in designing fibers with enhanced TMI thresholds.

As our formalism does not assume any particular fiber cross-sectional geometry, unlike most previous approaches,^{18,19,69} it can be utilized straightforwardly for calculating the TMI threshold for different fiber geometries under multimode excitation. We consider two additional fiber cross-sectional geometries other than circular: square fiber and D-shaped fiber (referring to a circular shape truncated by removing a section bounded by a chord). For simplicity, we ignore gain saturation and pump depletion in this section as these effects do not change the TMI threshold scaling qualitatively as shown in the main text. A reason to study these particular shapes is because square and D-shaped cross sections support modes with statistically different profiles, while a square shape results in an integrable transverse “cavity,” supporting modes with a regular spatial structure^{101} (as does the circular fiber); in contrast, the D-shaped cross section leads to wave-chaotic behavior in the transverse dimensions, with highly irregular and ergodic modes.^{102} Due to this property, D-shaped micro-cavities have found applications in suppressing instabilities and speckle-free imaging in 2D semiconductor micro-lasers.^{103,104} Similar to the circular fiber, the square fiber is chosen to have a core width of 40 *μ*m and a cladding width of 200 *μ*m. For the square fiber, we consider a slightly smaller NA of 0.135 as it has a slightly larger core area compared to the circular fiber such that it also supports 82 modes per polarization. In the D-shaped fiber, the core shape is formed by slicing a circle of diameter 40 *μ*m with a line at a distance 10 *μ*m (half the circle radius) from the center. Similarly, the cladding is obtained by slicing a circle of diameter 200 *μ*m at the same relative distance. For the D-shaped fiber, we consider a slightly larger NA (0.17) as it has a slightly smaller core area compared to the circular fiber such that it also supports 82 modes per polarization. All other relevant parameters for calculation are kept the same for all three fibers and are given in Table I.

Similar to the circular fiber, we first calculate the thermo-optical coupling coefficients for all mode pairs for both the square fiber and D-shaped fiber by using the formula in Eq. (25). We calculate the optical modes in each case by using the wave-optics module in COMSOL.^{75} Example of an optical mode for each fiber is shown in Fig. 8(a). As expected, the optical modes for the square fiber are highly structured, whereas optical modes for the D-shaped fiber are irregular. We also calculate first 1000 temperature eigenmodes for each fiber by using the Coefficient-Form-PDE module in COMSOL.^{75} The temperature eigenmodes depend on the cladding shape, which in this case is chosen to be the same as the respective core shape. An example of a temperature mode for each fiber is shown in Fig. 7(b). The temperature eigenmodes have similar spatial properties to optical modes, except that they are spread out over both the core and cladding. We numerically evaluate the overlap integrals of the optical and temperature modes to find *χ*_{mn}(Ω). We take the peak values over frequency for each mode pair to obtain effective coupling matrices, which are shown in Figs. 8(a) and 8(b). As expected, both fibers produce banded coupling matrices with a small number of significant elements. In each case, *χ*_{21} is the largest element, suggesting that FM-only excitation has the lowest TMI threshold and multimode excitation will lead to a higher TMI threshold.

To verify this, we calculate the TMI threshold for both the square and D-shaped fibers for FM-only excitation and equal mode excitation with the increasing number of modes. A threshold increase factor (TIF) is defined by taking the ratio between the TMI thresholds for multimode excitation and FM-only excitation. The results are shown in Fig. 5(c), where we have also plotted results for the circular fiber as the dotted curve for comparison. For both the square and D-shaped fibers, the TMI threshold increases linearly with the number of excited modes, similar to the circular fiber. This is in accordance with our reasoning based on the banded nature of the coupling matrix. The slope of the threshold scaling is highest for the D-shaped fiber, leading to more than a 14× higher TMI threshold, when all 82 modes are equally excited. The square shape leads to a slightly lower enhancement compared to the standard circular fiber. The superiority of the D-shaped fiber can be attributed to a lack of regular structure in the higher-order modes. As a result, unlike the more regular shapes, here, no single temperature eigenmode has a particularly strong overlap with a given mode pair; instead, many temperature eigenmodes with different eigenvalues participate in the coupling, leading to a broader TMI gain spectrum with a lower peak value.

## REFERENCES

*Fiber Lasers: Basics, Technology, and Applications*

^{3+}-doped few-mode fiber amplifiers

*Heat Conduction Using Greens Functions*

*Finite Difference Methods for Nonlinear Evolution Equations*

*Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems*