Theory of Transverse Mode Instability in Fiber Amplifiers with Multimode Excitations

Transverse Mode Instability (TMI) which 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 under arbitrary multimode input excitation for general fiber geometries. We show that TMI results from 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 $\sim$$10^4$ 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 single mode, resulting in significant TMI suppression. We find that the TMI threshold increases linearly with the number of modes that are excited, leading to more than an order of magnitude increase in the TMI threshold in a 82-mode fiber amplifier. Using our theory, we also calculate TMI threshold in fibers with non-circular geometries upon multimode excitation and show the linear scaling of TMI threshold to be a universal property of different fibers.


I. INTRODUCTION
Fiber lasers based on multi-stage fiber amplifiers (FA) 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 realize fully 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][12][13][14] and stimulated Brillouin scattering (SBS) [15][16][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][19][20][21][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][24][25][26][27][28][29][30][31][32][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][40][41][42][43][44] , and utilizing multicore fibers [45][46][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][35][36][37][38][39][40][41][42][43][44][45][46]48 . Thi 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 (SLM), 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][52][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][56][57][58] .The viability of such an approach has recently been demonstrated for other nonlinear effects in fibers such as SBS [59][60][61] and stimulated Raman scattering (SRS) 62 .
The current authors 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][65][66][67][68][69][70][71] .It has been shown that TMI can be modelled as 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 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 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.
Here we develop 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 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 ∼ 10 4 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 the final section of the 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.

A. Coupled Amplitude Equations
TMI is a result of dynamic transfer of power between fiber modes due to the thermo-optical scattering 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 un- dergo optical amplification and generate the output.But 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 in- terfere to create a spatio-temporally varying optical intensity pattern, which due to the quantum-defect heating that accompanies stimulated emission, results in dynamic temperature variations 24 .These spatio-temporal thermal fluctuations result in refractive index variations, which can cause significant transfer of power from signal at ω 0 in one mode to the noise at Stokes-shifted frequencies ω 0 − Ω in other modes.As a re- sult, 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 The signal in each mode m at ω 0 transfers power to the stokes shifted frequency in every other mode n, which causes the noise to grow exponentially with a growth rate depending on the signal power and pairwise thermo-optical coupling coefficients.When the noise at the output becomes a significant fraction of the signal, the output beam becomes unstable with dynamic fluctuations on a millisecond timescale.The output power at the onset of significant fluctuations (typically set at 1%) is defined as TMI threshold.
of 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 : where E is the total electric field , c is the speed of light and n 0 is the linear refractive index of the fiber.The temperaturedependent polarisation is given by: Here, η is the thermo-optic coefficient of the fiber material and ∆T represent the local temperature fluctuation due to quantum-defect heating.∆T satisfies the heat equation with a source term proportional to the local intensity 18 : The heat source Q is proportional to the local optical intensity I, the linear gain coefficient g, and the quantum defect q D = ( λ s λ p − 1), 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 Eq. 1 and Eq. 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 as a sum of the product of transverse mode profile, a rapidly varying longitudinal phase term, and slowly varying amplitude (SVA) for each fiber mode: (5) Here, A m (Ω, z) is the SVA in mode m at a point z along the fiber axis and for a Stokes frequency shift, Ω from the central frequency ω 0 .The Ω = 0 amplitude corresponds to the signal, and Ω = 0 corresponds to the noise.ψ m 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 : The solutions to Eq. 6 correspond to the guided fiber modes, which are confined to the fiber core by total internal reflection.These modes can be calculated using a numerical solver such as COMSOL 75 for an arbitrary fiber cross-section geometry.As an example, in Fig. 2a we have shown electric field profiles for the FMs and HOMs for a circular step-index fiber.
The width of the core in each case in 40 µm and the cladding width is 200 µm.The FM have no node in the electric field profile and follow the symmetries of the cross-section geometry.Each HOM can be labelled as linearly polarized (LP) mode and characterized by two indices (u, v) in the standard notation 73 .The mode LPuv 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 2a shows the profiles of LP01 (FM) and LP44 (HOM) modes.In appendix A, we study fibers with different geometries (D-shaped and Square cross-sections) and the associated optical modes are shown in Fig. A1a.
Normalizing each mode profile ψ 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 |A m (z, Ω)| 2 .The heat source term Q in Eq. 3 can be obtained in terms of optical modal amplitudes using Eq. 4 and Eq. 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 i j = β j − β i for each mode pair {i, j}.Thus, we expand ∆T into a series of temperature profiles corresponding to each heat source term: Here T i j ( r, Ω) is the profile of the temperature grating with longitudinal wavevector q i j and temporal frequency Ω.We further decompose each temperature profile T i j in terms of the eigenmodes of the transverse Laplacian ∇ 2 T : Here, a k i j is the SVA for temperature eigenmode 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 leads directly to a nonlinear thermo-optic coupling which is short-range in the momentum difference between the optical modes.T k denotes the transverse profile of temperature eigenmode k which satisfies the following eigenvalue equation: ∇ 2 T is a self-adjoint operator, thus its eigenmodes { T k } form a complete and orthogonal basis, with −α 2 k as the eigenvalue for eigenmode k.We impose Dirichlet boundary conditions at the outer surface of the fiber, corresponding to the standard assumption 19 that 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. 2b 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 are similar to the optical modes since both are eigenfunctions of same spatial operator ∇ 2 T .An important distinction is that optical modes are localised 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 follow the symmetries of the cladding geometry.The higher order eigenmodes are characterised by u azimuthal and v radial nodes across the entire fiber cross-section (i.e., core plus cladding), with a radial dependence given by (u + 1) th Bessel function of the first kind and angular dependence given by sine or cosine with u zeros.
To complete the solution for ∆T , we need to determine each coefficient a k i j .This is done by simplifying the LHS of Eq. 3 by substituting ∆T from Eq. 7 and Eq. 8 and using the eigenmode equation Eq. 9.The RHS of Eq. 3 can be obtained in terms of optical amplitudes by using Eq. 4 and Eq. 5. Finally, by exploiting orthogonality, we isolate the desired coefficient a k i j by multiplying both sides of simplified Eq. 3 with T k * and integrating it across the fiber cross-section leading to: The amplitude of temperature eigenmode k for the heat source term {i, j} is proportional to the overlap of the dot product of ψ i and ψ j with the temperature mode profile T k 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 Γ 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 Γ i j k ≈ Γ k = α 2 k D, 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.
Next, we use the solution for ∆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 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 simplified Eq. 1 with ψ * m and integrating it over the fiber cross-section: The 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 convolu- tions of three modal amplitudes A n , A * i , 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 the Green's function of the heat equation with the relevant optical mode profiles: The sum over k represents the sum over all the temperature eigenmodes.Thus, G mni j (Ω) denotes the optical modal overlap with the 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 mni j (Ω) 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 mni j vanishes and thus there is no transfer of power between various modes of the fiber at ω 0 .A nonzero Stokes fre- quency 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 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 and ∆T .The fiber properties are taken into account through χ 0 and G mni j via optical mode profiles {ψ m } and temperature mode profiles { T k }.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 1-D 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 the next section.

B. Phase-Matched Noise Growth
In this section, we derive an approximate solution to the coupled amplitude equations to study noise growth due to thermo-optical scattering when signal power is launched into multiple fiber modes.We utilize two key approximations: (1) we retain only the phase-matched terms 18,19,31,71 , since phasemismatched 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 {A s m (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: Similarly, the total amplitude in mode m at any point z can be decomposed into the signal (Ω = 0) and noise (Ω = 0) amplitudes: As the light propagates down the fiber, both the signal and noise grow exponentially due to the linear optical gain g.Below the TMI threshold, signal amplitudes roughly grow independently of the noise and in the case of no gain saturation are given by: Here, φ NL m is the nonlinear phase evolution in mode 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: We have utilized m = j, i = n as a solution to the phasematching 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 i j ≈ 0), thus we do not retain these terms.We have also assumed 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 needs to be modified slightly.For a detailed discussion see Supplementary Information Section I. Next, we substitute the ansatz for A m from Eq. 14 and Eq. 15 in Eq. 16 and simplify the RHS to obtain the equation growth equation for B m (z, Ω): All the terms which are quadratic or higher order in noise amplitudes are ignored.Each noise amplitude grows exponentially, independently of other noise amplitudes, and has contributions in growth rate from both linear amplification and nonlinear power transfer from the signal.P s m represents the signal power in mode m and is given by |A s m | 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. S.2.We can convert noise amplitude growth equations for mode m into growth equations for noise power (given by |B m | 2 ), by multiplying with B * m on both sides and adding the complex conjugate term: Here, P N m denotes the noise power in mode m.The first term on the RHS in Eq.S.2 corresponds to 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 .Eq. 18 can be solved analytically, resulting in exponential growth in noise power in each mode m: P N m (0, Ω) represents noise power in mode m at the input end of the fiber.P N m (L, Ω) 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: The z integral can be simplified by using the signal growth formula given in Eq. 15, L 0 gP s n (z)dz = P s n (L)− P s n (0) ≈ ∆P Ps n , where ∆P is the total signal power extracted from the amplifier and Ps n is the fraction of signal power in mode n (∑ n Ps n = 1), which is determined by the input excitation.Thus, TMI gain in mode m is given by: This formula suggests a relatively straightforward interpretation.The TMI gain for any mode depends on the total extracted signal power and effective thermo-optical coupling coefficient χm (Ω), which is given by the weighted sum of thermo-optical coupling coefficients with all the other modes, with weights depending on the input excitation.For a twomode fiber with FM-only excitation ( Ps 1 = 1), the above formula reduces to the well known formula for TMI gain 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 the next section.

C. Multimode TMI Threshold
The TMI threshold is defined as the total output power when the noise power becomes a significant fraction (ξ > 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 Eq. 19 and Eq.21: For each optical mode pair only a few temperature eigenmodes have significant overlap (dotted bars), the ones which have spatial frequency similar to the frequency of interference pattern of optical modes within the fiber core.For neighboring optical mode pair (LP01-LP11), significant overlap occur with lower order temperature eigenmodes (with high mode lifetime) and for mode pair with relatively larger separation (LP01-LP12), significant overlap occurs with relatively higher order temperature eigenmodes (with short mode lifetimes).Note that the thermo-optical coupling coefficients are decreased when the mode lifetimes 4πΓ −1 k are smaller.Consequently, optical mode pairs with larger separation have weaker overall coupling as indicated when the thermo-optical overlaps are weighted by lifetime (solid bars), leading to the lower peak coupling and higher peak frequency as shown in (a).
P N (L) = P N (0) e gL e (P th −P s (0)) χ = ξ P s (0)e gL , which can be rearranged as: where, P th is the TMI threshold, L is the length of the fiber, P s (0) and P N (0) are the input signal and noise powers, respectively.We have introduced an overall thermo-optical coupling coefficient χ which is equal to the maximum value of the effective thermo-optical coupling coefficients across all the modes and Stokes frequencies: A key insight from the multimode TMI threshold formula in Eq. 23 is that TMI threshold is roughly inversely proportional to the overall thermo-optical coupling coefficient χ, which depends on both the fiber properties through χ mn and the input power distribution, Ps n .It also depends weakly (logarithmically) on the input noise power and the noise fraction ξ at which the threshold is set.Also, any increase in 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 Supplementary Information Sec.II.It can be easily verified that when only the fundamental mode is excited ( P1 = 1, Pn =1 = 0), the overall thermo-optical coupling coefficient, χ simplifies to χ 21 , reducing our multimode TMI threshold formula to a previously derived formula in studies which 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 the next section.

D. Thermo-optical coupling
The thermo-optical coupling coefficient between any modes m and n is proportional to the imaginary part of a particular component of the 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 ψ m and ψ n with the temperature mode profile T k .Each con- tribution 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πΓ −1 k .To illustrate the key properties of the thermo-optical coupling coefficient we calculate χ mn (Ω) for all the mode pairs of a cir- cular step-index fiber with core diameter of 40 µm, cladding diameter of 200 µm and 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 xdirection.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. 2a).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.
In Fig. S.3a we have shown the thermo-optical coupling coefficient between the FM (LP01) and three HOMs in increasing mode order (LP11, LP02, 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 increasing transverse spatial frequency mismatch between the modes (green and orange curves).The peak frequency and linewidth are higher for 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. S.3b.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πΓ −1 k (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. S.3b).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. S.3b).The thermal mode lifetime is shorter for temperature eigenmodes with higher transverse spatial frequency (as shown by red curve in Fig. S.3b), so the contribution of higher-order temperature eigenmodes is damped.Such intrinsic dampening of high spatial frequency contributions to the temperature are 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, peak thermo-optical coupling coefficient decreases with 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 our use of thermal eigenmodes (each with a definite spatial frequency) for expressing the temperature fluctuations.This is a generic result for all optical mode-pairs.To explicitly show this, we calculate the peak value of χ mn (Ω) for all ∼ 10 4 mode pairs in circular fiber .The resulting coupling matrix is shown in Fig. 4a 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 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 Each element m, n is a positive entry given by the peak value of thermo-optical coupling coefficient χ mn (ω).Self-coupling (m = n) is not considered for the reasons discussed in the phase-matching subsection (Sec.IIB).Only a few entries close to the diagonal have significant value, giving a banded matrix.(b) TMI threshold scaling in circular fiber with the number of equally excited modes M. Threshold increase factor is defined by taking the ratio of TMI threshold with the threshold for FM-only excitation.As the number of excited modes are increased, the TMI threshold increases linearly and can reach upto 13 times higher when all 82 modes are excited.In comparison, the TMI threshold increase upon single HOM excitation (red curve) is significantly smaller.
the beam profile fluctuates dynamically rendering it useless for many applications 11 .According to Eq. 23 and Eq. 24, the TMI threshold is inversely proportional to the overall thermooptical coupling coefficient χ 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 ( Ps 1 = 1) and noise is present in HOMs, and (2) Equal mode excitation : the input power is divided equally in M modes of the fiber ( Ps n = 1/M).We use Eq.24 to calculate χ 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 signal power in the FM (m = 1), giving χ = χ 21 .Since χ 21 is the highest coupling out of all mode pairs (Fig. 4a), 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., χ = ∑ n =2 χ 2n /M.Due to the banded nature of χ mn only a few elements in the sum are significant, leading to χ ≈ sχ 21 /M.Here, s is used to denote the average number of significant elements in the any row of the thermo-optical coupling matrix, which does not scale with M and is roughly equal to 6 for 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.Fig. 4b shows the values of TIF as M is varied.As predicted, the TMI threshold increases linearly with 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 is equally efficient for considering fibers with non-circular cross-sections, such as the D-fiber with more complex spatial structures of its modes, due to the underlying chaotic ray dynamics, as discussed and depicted in Fig. 3.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 A, where we find linear scaling of the TMI threshold with the number of modes excited for the D-fiber.The slope of the scaling line does depend weakly 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 compo-sition, geometry (see above) and the precise distribution of excited modes.The level of enhancement will depend mainly 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 1-D cross-section waveguides 63 .Additionally, several recent experimental studies on the effect of fiber bending on TMI threshold in fewmode fibers provide evidence for our predictions [80][81][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 bendloss 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 the previous section, 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 semianalytical 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 the previous section (II E) 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 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 the previous section, 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
We consider a co-pumped fiber amplifier with ytterbium (Yb) doped gain medium 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, A s m , can be written as 13 : where, g mn is the gain in mode m due to amplitude in mode n and is given by: Here, the integral is over the doped area of the fiber crosssection.The new intensity-dependent factor in the denominator modifies the linear gain by taking into account the reduction in the inversion (gain saturation) due to removal of the pump energy by the signal amplification 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 Supplementary Information Section III).Notice that in addition to self-gain terms (m = n), we also must 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 specklelike 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.
As the pump beam excites the gain medium it loses energy, leading to a decreasing pump power along the fiber axis, referred to as pump depletion 13 .We take this into account quantitatively by introducing a standard equation for the evolution of pump power, which is given by 13 : The pump power decays exponentially with a longitudinally varying loss coefficient g p which depends on the local inversion and is given by 13 : Here, the integral is over the doped area of the fiber crosssection.σ e p and σ a p are the emission and absorption cross-sections of the pump respectively, N Yb is the density of doped ytterbium atoms and A cl is the area of pump core which is typically 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 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 σ e s and σ a s 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 S.1) coupled with the evolution of the pump power (Eqs.28 and 29) and upper level population (Eq.S.2) 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 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 value of g p , g mn and n u are updated directly from Eqs. S.1, 29 and S.2.To validate our numerical model, we first simulated the 20/400 fiber amplifier studied in Smith et al. 64 where single mode excitation was considered and obtained excellent agreement with their results.We also simulated the two-mode excitation case discussed in Li et al. 82 and found good agreement.
Next we study the highly multimode step-index fiber amplifier discussed in section II D. All the parameters for Yb doped gain medium are taken from Table I in Smith et al. 64 .The radius of the Yb doped area was considered to be equal to the fiber core radius.The unsaturated gain coefficient g 0 ≈ 4m −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 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 converted to the signal power.Fig. 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 ( λ p λ s = 0.9457).Notice that the signal power in different modes grow 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 the next subsection, 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
The gain saturation impacts the local inversion and therefore also changes the amount of heat generated due to the quantum defect.The gain-saturated heat source is given by 66,71 : where, 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,t).At any point, the signal intensity can be written as a sum of a static contribution I 0 ( r) and a dynamic contribution Ĩ( r,t) i.e.I( r,t) = I 0 ( r)+ Ĩ( r,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 a static and a dynamic contribution, Q( r,t) = Q 0 ( r) + Q( r,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 by using the Taylor expansion in terms of Ĩ/(I 0 + I sat ) and keeping the leading order term, giving the following expression 66 (for a detailed derivation, see Supplementary Information Section III): The dynamic heat source is proportional to the time varying intensity Ĩ 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 Ĩ. 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 et al. 66 and later by several other authors 69,71,86 .
Comparing the form of the relevant heat load in Eqs. 4 and S.9, we can directly derive a modified formula for the thermooptical coupling χ s mn , taking into account the gain saturation: The key modification in the formula for thermo-optical coupling upon gain saturation is the presence of a saturation term in the overlap integral associated with the heat equation.Here, we have written the transverse integrals explicitly instead of using .as we did in Eq. 25 to highlight the space-dependent nature of the saturation term due to the spatial-hole burning.In addition, since the local intensity 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 χ s mn inside the z integral: Here, P s n is the signal power in mode n calculated with gain saturation in the first step discussed in the previous subsection.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 χ s mn .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: Here the denominator is for normalization, chosen such that it is equal to L 0 dzg s P s n (z), ensuring that in the absence of gain saturation the 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, ) 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 the previous section everywhere along the fiber when all 82modes were equally excited and calculated the z-integrated coupling matrix in the presence of gain saturation.The results are shown in Fig. 6a.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.6a and S.3a 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 the next subsection.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
As just noted, gain saturation affects the local heat load and reduces the resultant thermo-optical coupling, leading to an increased value of TMI threshold, a fact that has been demonstrated both theoretically and experimentally for single mode excitations 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 χ 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 χ, which depends strongly on the input excitation via signal power distribution in various modes { Ps n }.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 Smith et al. 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 W 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 are increased, the TMI threshold also increases linearly especially for large number of modes and reaches close to 5700 W when all 82 modes are equally excited (Fig. 6b).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 W to 800 W for various number 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 which 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 low number of modes causing deviations from strict linear scaling at the lower end in Fig. 6b.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 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 spatially coherent output that can be focused to a diffraction-limited spot or collimated to a Gaussian beam.By contrast, our approach allows the generation of 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 mis-takenly assumed that the presence of multiple spatial modes necessarily leads to a poor output beam quality and a high value of M 2 . 49,50However, this belief is not correct, as previous work on the manipulation of coherent optical fields has 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-meter-long silica fiber with an NA of 0.1, the spectral correlation width is approximately 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 depends mostly on the pump power.Since focusing to diffraction-limited spot in the near field necessarily leads to 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 was recently demonstrated experimentally for the case of Stimulated Brillouin Scattering (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 thermooptical 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 5 mode excitations in 1-D 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][81][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 which 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-section geometry.In Appendix A 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][60][61] , stimulated Raman scattering (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 multimode TMI threshold formula and justify the use of only the Stokes mode with maximum gain at peak frequency.In section III, we provide a derivation of the saturated gain coefficient along with the formulae 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.
qualitatively as shown in the main text.A reason to study these particular shapes is because square and D-shaped crosssections support modes with statistically different profiles.While a square shape results in an integrable transverse "cavity", supporting modes with regular spatial structure 101 (as does the circular fiber); in contrast, the D-shaped crosssection 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 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 centre.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 slightly smaller core area compared to 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 thermooptical 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. 8a.As expected, the optical modes for the square fiber are highly structured whereas optical modes for 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 same as the respective core shape.Example of a temperature mode for each fiber is shown in Fig. 7b.The temperature eigenmodes have similar spatial properties to optical modes except 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 Fig. 8a and Fig. 8b.As expected both fibers produce banded coupling matrices with a small number of significant elements.In each case χ 21 is the largest element, suggesting FM-only excitation have 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 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. 5c, where we we have also plotted results for circular fiber as 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.shown in Refs.? and ?that such terms along with other nonphase matched terms can actually lead to a behavior similar to modulation instability, even at zero Stokes frequency shift, at really high powers.But the threshold for such a process is typically much higher than that due to the phase-matched dynamic transfer of power leading to TMI in fibers with diameters smaller than 80 µm ? .

C. Impact of degeneracies
While retaining the phase-matched terms in Eq. [16] in the main text, we also made an assumption that none of the fiber modes are exactly degenerate.This is justified because most high-power fiber amplifiers utilize step-index fibers and have imperfections and spatially inhomogeneous temperature distribution which typically break the exact circular geometry of the refractive index distribution in the fiber, leading to a lifting of mode degeneracies.However in the case when exact or nearly exact degeneracies are present, our formalism needs to take these into account, which can be achieved rather straightforwardly.We can replace Eq. [17] in the main text with a slightly modified equation for the Stokes growth in each mode: 2) Here we have introduced the index σ for labelling modes within a degenerate subspace, in addition to the mode index n that identifies various non-degenerate mode groups.For instance in perfect circular step-index fibers, the LPmn modes (for m = 0) are two-fold degenerate, due to two solutions in the azimuthal direction (sine and cosine) for the same propagation constant, hence we have σ = 1, 2. The key consequence of modal degeneracy is that the equations for the growth of the Stokes power in the modes within a degenerate sub-group become coupled.The stokes growth rate for each sub-group is given by the eigenvalue of the growth matrix whose elements are given by the weighted sum of the thermooptical coupling (given by G σ σ σ ′ σ ′ mnnm ) with weights equal to the signal amplitudes in various modes.Note that the equations for the Stokes growth in various non-degenerate mode groups are still uncoupled.This leads to an overall growth matrix for the Stokes amplitudes which has a block diagonal form with relatively small block sizes and can be easily diagonalized.In the case of sine/cosine degeneracy in LP modes, the size of the blocks is 2 × 2. For the graded-index fibers, the role of degeneracies can be more important since the size of degenerate sub-groups increases with order of the modes.
An interesting consequence of the coupling of the Stokes growth equations for degenerate modes is that the Stokes growth rate typically becomes larger, as it is mainly determined by the maximum eigenvalue, which is increased by the presence of cross-couplings, while the smallest eigen-value is suppressed.For 2-fold LP mode degeneracy, the TMI gain would be increased by a factor of 2. However, since such an effect would impact the TMI gain equally for both fundamental-mode excitation, where an LP01 mode couples with two degenerate LP11 modes, and for multimode excitations; hence it does not change the enhancement factor of TMI threshold upon multimode excitation.Note that the coupling to exactly degenerate HOM modes has never been considered in previous theories for TMI which assume single mode excitation, and some of these theories have shown good agreement to experimental results.A possible explanation for this is the likely breaking of exact degeneracies in realistic fiber amplifiers, as mentioned earlier

S.2. MULTIMODE TMI THRESHOLD
In section II.C of the main text we have defined the multimode TMI threshold as the total output signal power at which the Stokes power reaches a significant fraction (> 1%) of the signal power.The noise power at the input is assumed to be seeded either by relative intensity noise due to signal linewidth or other experimental imperfections.In absence of these effects, the noise power has a lower limit given by the quantum noise and is approximately 10 −16 W per stokes mode in a KHz frequency window.We approximate the total noise power at the output by the power in the mode with maximum Stokes gain at the peak Stokes gain frequency.Although the noise power is present in multiple modes and frequencies, the approximation of considering only the mode and frequency with peak gain is still justified due to the exponential nature of the Stokes growth.Here we confirm the validity of this approximation both mathematically and with a concrete example.At the TMI threshold, the total noise power at the output is given by: dΩ Pm (0, Ω)e gL e χm (Ω)(P th −P(0)) . (S.1) Here χm (Ω) is the effective thermo-optical coupling for mode m, defined in Eq. [21] in the main text.We assume a uniform density of input noise power ( Pm (0, Ω) = PN (0)) across a relevant frequency window in all fiber modes.Since the integrand in Eq.S.1 is an exponential function we utilize saddle-point method to approximate the integral, leading to: 0)e gL 2π χ′′ m (Ω p )(P th − P(0)) e χm (Ω p )(P th −P(0)) .
(S.2) The total noise power depends exponentially on the total extracted signal power (P th − P(0)) and the effective thermooptical coupling coefficient at the peak gain frequency Ω p .Here PN (0) is the input noise power density with the units of [W/Hz].The square root term in the prefactor depends on curvature (given by the second-order derivative) of the TMI gain curve and has units of [Hz].Further simplifications of the prefactor can be achieved by approximating χm (Ω) Scaling of TMI threshold with the number of equally excited modes in a multimode fiber The approximate model (solid black curve) which only considers the Stokes power in the mode with maximum growth rate at peak gain frequency produces a very similar scaling as the full model (dashed red curve) which considers the total Stokes power.This is a result of the exponential nature of the Stokes power growth starting from a very small seed power.
with a parabolic approximation around ω p giving χ′′ m (Ω p ) = χm (Ω p )/(∆Ω) 2 , where ∆Ω is the relevant frequency bandwidth for the thermo-optical coupling.This leads to the more familiar form for the Stokes power used in the main text: P N (L) ≈ ∑ m P N (0)γe gL e χm (Ω p )(P th −P(0)) . (S.3) Here, P N (0) = PN (0)∆Ω is the total input noise power per mode in the relevant frequency window, and we have defined a unitless quantity γ = 2π/ χm (Ω p )(P th − P(0)), which is roughly fixed at the TMI threshold with a value on the order of unity.For our fiber, γ ≈ 0.5.In any case, the prefactor in the above equation plays little role in determining the TMI threshold, which primarily depends on the argument of the exponential.Further, the sum over the Stokes modes can be approximated by the largest contribution or a few largest contributions in the case that multiple Stokes modes have the same gain, leading directly to Eq. [22] in the main text.This approximation is equivaluent to replacing the LogSumExp function (also called RealSoftMax) by the max function.It is known that LogSumExp is tightly bounded above and below by the max function up to small logarithmic corrections, which rigorously justifies our approximation ? .
To validate the approximation of the total Stokes power by the Stokes power in the mode with maximum Stokes gain at the peak Stokes gain frequency, we computed the scaling of TMI threshold using both the exact expression for total Stokes power Eq.S.1 and the approximate expression used in the main text Eq. [22].The results are shown in Fig. S.1.It can be seen that the scaling of the TMI threshold computed by the two expressions is nearly identical.Approximating the Stokes power with the power in the mode with maximum Stokes gain at peak frequency is quite useful, as it allows a straightforward formula for TMI threshold in terms of a single axiallyintegrated thermo-optical coefficient, which linearly depends on the signal power in various modes.

S.3. GAIN SATURATION A. Gain Saturation parameters
In Eq. [27] in the main text, the following formula is used for g mn , the gain in mode m due to complex amplitude of field in mode n in presence of gain saturation: Here we derive this formula along with the formulas for unsaturated gain coefficient g 0 and signal saturation intensity I sat from the steady-state solution to the rate equations for Yb doped gain medium.In the steady state, the fraction of the population in the upper lasing level n u is given by 13,87 : Here, ω s is the signal frequency, ω p is the pump frequency, σ e s and σ a s are the emission and absorption cross-sections of the signal, σ e p and σ a p are the emission and absorption cross sections of the pump, P p is the local pump power, and I 0 is the local signal intensity.This formula can be rewritten as follows: where P sat p is the pump saturation power and I 0 sat is signal saturation intensity at zero pump power and these are defined as follows: .
(S.4) Since σ a s ≪ σ e s , the second term in the numerator in Eq.S.3 can be ignored, especially when significant pump power is present in the fiber.Next, we divide both the numerator and denominator in Eq.S.3 by 1 + P p P sat p to obtain a familiar form for the upper level population: where n u is the maximum upper level population (occurring when the signal intensity is zero) and I sat is the pumpdependent signal saturation intensity and these are given by: Note that the peak fractional upper level population, which defines the unsaturated gain coefficient g 0 depends on both the pump power, and the properties of the gain medium.At low pump powers, n 0 u is directly proportional to the pump power but for P p ≫ P sat p , n 0 u is roughly constant.The formula for signal saturation intensity also shows a dependence on the pump power along with various properties of the gain medium.At small pump power (P p ≪ P sat p ), I sat is roughly constant and is equal to I 0 sat .However in high-power fiber amplifiers typically the starting pump power is significantly higher than P sat p , resulting in I sat varying linearly with the pump power.Hence pump depletion and gain saturation are interconnected in such amplifiers.The upper level population can be used to calculate signal gain as follows: (S.7)Here N Yb is the density of doped Ytterbium ions.The term involving the lower level population n l is ignored since σ a s ≪ σ e s .Substituting the expression for n u from Eq. S.5 into Eq.S.7, we obtain the formula for g mn given in Eq.S.1, and the unsaturated gain coefficient g 0 is equal to: In Eq. [32] in the main text, we provide the following formula for the saturated dynamic heat load: Q( r,t) ≈ g 0 q D Ĩ( r,t) (1 + I 0 ( r)/I sat ) 2 . (S.9) The quadratic power of the saturation term in the denominator was first pointed out by Hansen et al. 66 and is responsible for the increased TMI threshold due to the gain saturation.Here we provide a detailed derivation for this formula.The total gain-saturated heat load is given by 66 : Q( r,t) = g s q D I( r,t) = g 0 q D I( r,t) 1 + I( r,t)/I sat , (S.10) where, g s and g 0 are the saturated and unsaturated gain coefficients and q D is the quantum defect.At any point, the signal intensity can be written as a sum of a static contribution I 0 ( r) and a dynamic contribution Ĩ( r,t), i.e., I( r,t) = I 0 ( r) + Ĩ( r,t).Substituting this in Eq.S.10 and rearranging the denominator we obtain: Q( r,t) = g 0 q D (I 0 + Ĩ) (1 + I 0 /I sat )(1 + Ĩ/(I 0 + I sat )) = g 0 q D (1 + I 0 /I sat ) (I 0 + Ĩ)(1 + Ĩ/(I 0 + I sat )) −1 .(S.11) Here we have suppressed the arguments of all the functions for compactness.Below the TMI threshold the dynamic signal intensity Ĩ is much smaller than I 0 and I sat , therefore we can Taylor-expand the last term in the numerator: Q( r,t) = g 0 q D (1 + I 0 /I sat ) (I 0 + Ĩ)(1 − Ĩ/(I 0 + I sat ) + O( Ĩ2 )).
(S.12) Multiplying the terms in the numerator and ignoring the terms which are O( Ĩ2 ) and higher, we get: Q( r,t) ≈ g 0 q D (1 + I 0 /I sat ) (I 0 + Ĩ(1 − I 0 /(I 0 + I sat ))) = g 0 q D I 0 (1 + I 0 /I sat ) + g 0 q D Ĩ (1 + I 0 /I sat ) 2 (S.13) The first term on the right hand side of the above equation describes the static heat load, and the last term gives the dynamic heat load, which has an exact form used in the main text and shown in Eq.S.9.S.2.Saturated thermo-optical coupling χ′ between the fundamental mode and first higher order mode (HOM) computed using our TMI model (solid curves) and compared to detailed simulations in Ref. 65 (dotted curves).Fibers with two different cladding radii (with different gain saturation amount) are considered.Our results are in close agreement with those in Ref. 65 .

S.4. COMPARISON WITH PREVIOUS STUDIES
In this section, we perform feasible comparisons/validations to experiment and previous theories as described below.
Using our model, we calculated the TMI threshold for the co-pumped 50/400 fiber discussed in Ref. 86 operating at 1064 nm signal wavelength and obtained a value of 680 W which closely matches with the experimentally obtained value of 658 W shown in Table II in Ref. 86 .For the noise seed power, we utilize the value of relative intensity noise (RIN) provided in Table I (R N = 10 −10 Hz −1 ) in Ref. 86 .The axially-integrated TMI gain to reach the TMI threshold is roughly 13.8.Since there are few existing experimental or detailed numerical studies of TMI for highly multimode excitations prior to our work we have validated our model against existing studies on TMI for single mode and two-mode excitations.Two examples of detailed validation of our model are presented below: (2) Single mode excitation: Ref. 65 .
One of the earliest studies on the impact of gain saturation on the TMI threshold was published by Smith et al. 65 Using detailed numerical simulations, they studied the thermooptical coupling χ′ between the fundamental mode and the first higher order mode in a fiber with core diameter of 50 µm.It was shown that in the presence of gain saturation the effective thermo-optical coupling coefficient is lowered due to the reduction of dynamic heat load.The results from Fig. 5 in Ref. 65 are reproduced here as dotted curves in Fig. S2.We use our model to calculate the effective thermo-optical coupling for a co-pumped fiber with same parameters as used in Ref. 65 .Our results (solid curves in Fig. S2) are in excellent agreement with the detailed simulation results in Ref. 65 .When a fiber with larger cladding (pump core) diameter is used (Fig. S2b) the amount of gain saturation is higher as the pump intensity in the gain core is lowered.This leads to an increased reduction in the thermo-optical coupling and a higher TMI threshold.We also calculate the TMI threshold for the both the fibers (d clad = 100µm and d clad = 400µm) using our model and the values (488 W and 1101 W) match closely with the values (495 W and 1100 W) in Ref. 65 .For noise seeding, we consider quantum noise (P N = 10 −16 W ) to match the parameters given in Table I in Ref. 65 .The axially-integrated TMI gain to reach the TMI threshold is roughly 34.5.
(3) Two mode excitation: Ref. 82 A recent study which utilized detailed simulations to calculate the TMI threshold upon two-mode excitation was presented in Li et al. in Ref. 82 .The gain fiber was operated at a signal wavelength of 1080 nm and a pump wavelength of 1018 nm.By varying the higher order mode (HOM) loss coefficient α they were able to achieve a two-mode excitation with variable HOM content.Fig. S3a shows the HOM content as α is increased calculated using our model for the co-pumped case and that presented in Fig. S3b in Ref. 82 .As α is increased HOM content decreases and the values calculated with our model match reasonably with Ref. 82 .For each value of α the TMI threshold is obtained with our model, and it agrees with the prior result.For noise seeding, we consider quantum noise (P N ≈ 4 × 10 −16 W ) to match the parameters given in Ref. 82 .The axially-integrated TMI gain to reach the TMI threshold is roughly 34.5.Figure S3b shows that for smaller α (higher HOM content) the TMI threshold is higher as the effective number of modes are higher.This is explicitly shown in Fig. S3c where we plot the TMI threshold as a function of HOM content which shows linear scaling with the HOM content and our results match closely with those in Ref. 82 .Note that the linear increase is expected to work only until the HOM content is less than 0.5 in the two-mode case after which the effective number of modes would start to decrease.

FIG. 1 .
FIG.1.Schematic of TMI in a multimode fiber amplifier.At the fiber input, a signal wave is launched in multiple modes at the same frequency, ω 0 , with a small amount of noise at the shifted Stokes frequency, ω 0 − Ω.As the signal propagates in the fiber, it undergoes amplification due to the stimulated emission, generating heat in the process and causing temperature fluctuations (shown in red in the inset).These temperature variations result in refractive index variations due to the thermo-optical effect, which cause scattering of power between various modes.(Note:longitudinal and transverse dimensions in lower panel are not to scale).(b) The signal in each mode m at ω 0 transfers power to the stokes shifted frequency in every other mode n, which causes the noise to grow exponentially with a growth rate depending on the signal power and pairwise thermo-optical coupling coefficients.When the noise at the output becomes a significant fraction of the signal, the output beam becomes unstable with dynamic fluctuations on a millisecond timescale.The output power at the onset of significant fluctuations (typically set at 1%) is defined as TMI threshold.
FIG. 2. (a)Amplitude profiles of guided optical fundamental mode (FM) and a higher order mode (HOM) for a circular step-index fiber.The optical modes are guided in the core of the fiber.The FM (LP01) has no nodes whereas each HOM denoted as LPuv is characterized by u azimuthal nodes and v−1 radial nodes.Here LP44 is shown.(b) Amplitude profiles for thermal modes, which are the spatial eigenmodes of the heat equation with constant temperature at the outer cladding.Each thermal mode fills the entire fiber cross section.Similar to optical modes, FM has no nodes whereas each HOM denoted is characterized by an integer number of radial and azimuthal nodes.Cladding and core sizes are not to scale.
FIG. 3. (a)Thermo-optical coupling coefficient between the FM (LP01) and three HOMs in increasing mode order (LP11, LP02, LP12) for circular fiber .Relevant mode profiles are shown in the insets.The thermo-optical coupling coefficient vanishes for zero Stokes shift (Ω = 0) and peaks at frequencies on the order of few KHz.LP01 couples most strongly with LP11 (blue curve) and the peak coupling decreases with increasing transverse spatial frequency mismatch between the modes (green and orange curves).The peak frequency and linewidth also increases with mode order of HOMs.(b) thermo-optical mode overlap for two different optical mode pairs: LP01-LP11 (blue) and LP01-LP12 (orange) with various temperature eigenmodes.The temperature modes are ordered by increasing eigenvalues and decreasing mode lifetimes, 2πΓ −1 k (shown by red curve).For each optical mode pair only a few temperature eigenmodes have significant overlap (dotted bars), the ones which have spatial frequency similar to the frequency of interference pattern of optical modes within the fiber core.For neighboring optical mode pair (LP01-LP11), significant overlap occur with lower order temperature eigenmodes (with high mode lifetime) and for mode pair with relatively larger separation (LP01-LP12), significant overlap occurs with relatively higher order temperature eigenmodes (with short mode lifetimes).Note that the thermo-optical coupling coefficients are decreased when the mode lifetimes 4πΓ −1 k are smaller.Consequently, optical mode pairs with larger separation have weaker overall coupling as indicated when the thermo-optical overlaps are weighted by lifetime (solid bars), leading to the lower peak coupling and higher peak frequency as shown in (a).

FIG
FIG. 4. (a)Thermo-optical coupling matrix for circular fiber (Circular cross-section).Each element m, n is a positive entry given by the peak value of thermo-optical coupling coefficient χ mn (ω).Self-coupling (m = n) is not considered for the reasons discussed in the phase-matching subsection (Sec.IIB).Only a few entries close to the diagonal have significant value, giving a banded matrix.(b) TMI threshold scaling in circular fiber with the number of equally excited modes M. Threshold increase factor is defined by taking the ratio of TMI threshold with the threshold for FM-only excitation.As the number of excited modes are increased, the TMI threshold increases linearly and can reach upto 13 times higher when all 82 modes are excited.In comparison, the TMI threshold increase upon single HOM excitation (red curve) is significantly smaller.

FIG. 5 .
FIG.5.Pump power (solid red line, scale on the right y-axis) gets depleted as it is absorbed along the fiber creating inversion.This inversion leads to the growth of total signal power (shown by dotted red line, scale on the right y-axis) which grows exponentially at first and then linearly due to the gain saturation and eventually flattens out as most of the pump is depleted.All the other curves with various colors show the signal power in individual modes (scale on left yaxis) displaying the mode dependent gain due to spatial hole burning and differential overlap with the gain medium.
FIG.6.(a) z-integrated thermo-optical coupling matrix for fiber with circular cross-section including the effects of gain saturation and pump depletion.Self-coupling (m = n) is not considered for the reasons discussed in the phase-matching subsection (Sec.IIB).Similar to the unsaturated gain case, only a few entries close to the diagonal have significant value, giving a banded matrix.Saturation of the heat load leads to a reduction in overall thermo-optical coupling as displayed by the color-scale.(b) TMI threshold scaling in circular step-index fiber with the number of equally excited modes M with (blue circles) and without (black crosses) gain saturation.In both the cases the TMI threshold increases linearly with the number of modes.The dotted lines show best fit lines with slope = 54 W/mode without gain saturation and slope = 62 W/mode with gain saturation.When gain saturation is included the TMI threshold and the slope of linear scaling is higher.The difference in threshold in these two cases does not appear very dramatic because the y-axis spans a large range (0-6 kW) of values due to a large impact of multimode excitation.

FIG. 7 .
FIG. 7. (a) Amplitude profiles of guided optical modes for D-shaped and square cross-section step-index fibers.The optical modes are guided in the core of the fiber.(b) Amplitude profiles for thermal modes for fibers with D-shaped and square shaped cladding.The thermal modes are the spatial eigenmodes of the heat equation with constant temperature at the cladding boundary.Each thermal mode fills the entire fiber cross section.Both optical and thermal modes in the square fiber are structured and have a particular number of nodes along x and y axes.The modes in the D-shaped fiber have a random structure as D-shaped cavities are wave-chaotic.

FIG. 8 .
FIG.8.Thermo-optical coupling matrix for (a) square fiber (b) D-shaped fiber.Self-coupling is not considered (m = n).Similar to circular fiber, the thermo-optical coupling matrices for both square fiber and C are also banded.(b) As a result TMI threshold in square fiber (red) and D-shaped fiber (blue) also increases linearly with the number of equally excited modes M. The slope of the scaling is highest for D-shaped fiber and lowest for square fiber.The scaling for circular fiber is reproduced as the dotted curve for reference.
FIG.S.1.Scaling of TMI threshold with the number of equally excited modes in a multimode fiber The approximate model (solid black curve) which only considers the Stokes power in the mode with maximum growth rate at peak gain frequency produces a very similar scaling as the full model (dashed red curve) which considers the total Stokes power.This is a result of the exponential nature of the Stokes power growth starting from a very small seed power.
FIG. S.2.Saturated thermo-optical coupling χ′ between the fundamental mode and first higher order mode (HOM) computed using our TMI model (solid curves) and compared to detailed simulations in Ref.65 (dotted curves).Fibers with two different cladding radii (with different gain saturation amount) are considered.Our results are in close agreement with those in Ref.65 .
FIG. S.3.Comparison between results of our TMI model and Ref. 82 by Li et al. for TMI threshold upon variable two-mode excitation.(a) As the HOM loss α is increased, the relative HOM content decreases.(b) TMI threshold increases with decreasing α as the excitation becomes more multimode (c) TMI threshold scales roughly linearly with the HOM content.Overall the results from Li et al. are in close agreement with our TMI model.

TABLE I
. Detailed parameters for circular fiber.
Yb σ e s n 0 u ≈ N Yb σ e B. Dynamic heat load saturation