1. Introduction

Due to their large surface area, fiber lasers can reach output powers in the kilowatt range while maintaining a small footprint and high reliability. As fiber technology has matured, fiber lasers have gained an ever-increasing market share as manufacturing tools, and they are commonly used in the growing field of directed energy weapons. After their initial demonstration, the output power of fiber lasers increased exponentially for a few decades, quickly reaching diffraction-limited time-averaged output powers that surpassed traditional lasers that used bulk, solid-state lasing media [1,2]. Although the geometry of fiber lasers makes them uniquely suited for efficient heat dissipation, the drive for ever-increasing output powers eventually gave rise to thermal nonlinear parasitic effects, posing a serious obstacle for further power scaling.

The first observations of transverse mode instability (TMI) were published over a decade ago, and there have been no new single-mode average output power records since [2,3]. TMI is a mode coupling effect that impacts high-power fiber laser systems. It results in a seemingly chaotic oscillation of the output beam mode structure between the fundamental mode (FM) and one or more higher order modes (HOMs), drastically reducing beam quality. It is widely accepted that TMI arises from the oscillatory heat pattern generated by the interference beating between the FM and one or more HOMs [1]. As the temperature profile evolves, this creates a thermo-optic long-period fiber grating that can transfer additional power from the FM to the HOMs. Since the HOM interference produces the temperature "wave" that couples more power into the HOM, it can be classified as a stimulated thermal Rayleigh scattering (STRS) process [4].

TMI is not the only power-limiting nonlinear effect in fiber amplifiers. Many experimental fiber laser systems are power limited by either TMI or stimulated Raman Scattering (SRS), and theoretical studies suggest that the ultimate nonlinear limit for standard master oscillator power amplifier (MOPA) systems is determined by the intersection of TMI and SRS [5,6]. Which nonlinear effect dominates will depend on the design of the laser and its operating conditions. In order to overcome SRS, many fiber designs used a wide core and a short length. Unfortunately, experimental evidence and theoretical arguments based on stability analysis show that the TMI threshold tends to decrease with increasing core diameter [5,7]. The reason for this decrease and the exact relationship between TMI threshold and core diameter has been debated for some time [8]. Recent modelling suggests that gain saturation is a primary factor in determining the TMI threshold, giving the threshold an inverse dependence on core area to first order [6,9].

Within the past decade, TMI mitigation has been an active field of research, and there are already several promising methods to reduce its effects [2,10–12]. However, some basic design considerations for TMI reduction are still poorly understood. Because TMI is an inherently high-power effect that generally will not be observed at fiber output powers below 100W, it is difficult and expensive to rapidly iterate through a wide range of fiber amplifier designs, making it difficult to quantify the dependence of the TMI threshold on fiber design parameters. This difficulty is exemplified by the long-term lack of consensus regarding the TMI threshold’s dependence on core diameter. Progress towards understanding the relationship between TMI and core diameter has already required multiple advancements in theory and modelling [4,6,9]. Similarly, future advancements in TMI mitigation will require simulation methods that are accurate, efficient, and representative of a wide range of fiber designs and operating conditions.

Many simulation models have been proposed and tested [2,9,13–18]. Some models operate in the time domain [2,15–17], others operate in the frequency domain [9,13,14], and a few assume steady-periodic operation [18]. Each of these methods comes with its own set of limitations. Steady-periodic methods assume the mode coupling is periodic, so they are unable to model transient effects or any active modulation that is not at a harmonic of the mode coupling frequency. Frequency-domain models are forced to either ignore or approximate the rare-earth dopant inversion dynamics. Because dopant inversion is a nonlinear function of the signal and pump intensities, there is no straightforward way to calculate exact, frequency-domain expressions for gain and heat generation. To get around this, some frequency-domain models have assumed the gain is constant [4,13], while some newer models use a small signal linearization to approximate the dopant inversion [9,14]. These approximations reduce accuracy and limit the ability of frequency-domain models to account for effects like spatial hole burning. In theory, time-domain models could account for the full physics of the fiber amplifier, and they allow for the modelling of transient effects and high-frequency modulation that would often be impractical or impossible to model with frequency-domain methods. Unfortunately, many time-domain models are too computationally intensive for practical simulation of fiber amplifiers with realistic dimensions [15]. More recently, simplified time-domain models capable of calculating TMI evolution with reasonable runtime and modest computational resources have been developed [2,17], but they still tend to be much more computationally intensive than their frequency-domain counterparts.

This paper presents a time-domain model that uses a minimal set of fiber thermal modes to further improve the computational efficiency of TMI simulation in $\text {Yb}^{3+}$ doped fiber amplifiers, significantly reducing both runtime and memory requirements. This method describes the optical power coupling as a three-wave mixing between two optical modes and an oscillating temperature field, which is represented using the finite set of fiber thermal modes. As such, it provides a clear description of the time-dependent TMI evolution as an STRS process.

2. Temperature eigenfunctions

Because TMI is a thermo-optic nonlinear process, the heart of any TMI simulation is a calculation of temperature evolution in the fiber as the signal amplification process generates heat. For time-domain simulation, this is most often done using a finite-difference time-domain (FDTD) solution to the heat transport equation [15,17]. Depending on implementation details, tracking this temperature evolution will likely be the most computationally intensive portion of the simulation. While FDTD makes sense for a completely general thermal evolution problem, TMI can only be induced by temperature distributions with certain features or symmetries; therefore, any portion of a distribution without these features can be ignored without loss of accuracy. It would be prudent, then, to find a computational method of tracking temperature that is able to ignore unnecessary thermal information and that has simple methods for calculating thermal evolution.

Rather than tracking temperature as a single function that varies in both space and time, it can be decomposed into a set of orthogonal temperature eigenfunctions (TEs). These are the temperature patterns in the cylindrical fiber that maintain their shape as they decay exponentially over time. They can be found by searching for temperature distributions of the form

Assuming a separable form for $\Delta T_{mn}$ in the azimuthal angle, $\phi$, and the radial distance from the fiber center, $r_\perp$, gives the normalized eigenfunction solution shown below

Then $\gamma _{mn}$ is found by applying an appropriate condition at the thermal boundary shown in Fig. 1. A convectively cooled boundary can accurately represent many practical scenarios.

 

Fig. 1. Diagram of thermal and optical boundary conditions used in TMI simulation.

Download Full Size | PDF

Solutions of Eqs. (3) and (5) can be used as an expansion to accurately represent the temperature in the fiber amplifier during operation:

The TEs described above are the standard thermal modes of a convectively cooled uniform cylinder. They have been used to model thermal effects in optical fiber as far back as 1998 [20], and they are the same basis used in the frequency-domain models [4,9,14] and the time-domain model [16]. The TE decomposition is also conceptually similar to the eigenvector method used in [2]. As described in the next section, the contribution of this work is the definition and use of a minimal set of TEs in an efficient, time-domain STRS model.

3. Eigenfunctions in time-domain simulation

For maximum efficiency, the model described below drops any coupling terms that are not phase matched. This simplification reflects identification of TMI as a three-wave mixing STRS process. This type of simplification has been previously used in the frequency domain by Hansen et al. and in the time domain by Menyuk et al. [13,17]. The time-domain model described below generally follows the notation of Menyuk et al. presented in [17] but substitutes TE expansion rather than the usual thermal FDTD calculation. It assumes operation at or below the TMI threshold, meaning that the dominant coupling is between the FM and the first HOM and that coupling to higher-order modes can be ignored. For a polarized signal beam, the field is then a simple expansion in these first two modes:

In Eqs. (8) and (9), $\omega$ is the angular frequency of the signal beam, $\beta = \beta _0 \approx \beta _1$ is the propagation constant of the FM, $\Delta \beta = \beta _0 - \beta _1$ is the difference between the FM and HOM propagation constants, and $n_0$ is the fiber index without perturbation. In a simple TMI model, the index perturbation, $\Delta n$, includes at least the thermo-optic effect and the gain. Following the contribution of Menyuk et al., these equations can be simplified using phase matching [17].

The equations above have dropped all terms that are not phase matched, since, under normal operating conditions, their contribution to mode coupling vanishes. They use a truncated expansion $\Delta n = \Delta n_0 + (1/2) [\Delta n_+ e^{{i} \Delta \beta z} + \Delta n_- e^{-{i} \Delta \beta z}]$. All $\Delta n_j$ are slowly varying fields. The $\Delta n_{\pm }$ describe the oscillating index perturbation that arises from the beating between the FM and the HOM. Terms with higher order z-dependence are dropped since they make no phase-matched contribution to mode coupling. Temperature, heat, and dopant inversion are also written in this slowly varying, truncated form.

Once the propagation coefficients $c_{ij}(z,t)$ are known, the CMT equations are easily solved numerically. Figure 2 is a flowchart for calculating the propagation coefficients as they evolve over time. Since the novel contribution of this work is the use of TEs rather than the standard heat diffusion FDTD, only the thermal evolution and the calculation of the $c_{ij}$ coefficients are described here. For expressions of the dopant inversion ratio ($n_2$), the heat generation, and the gain see the work of Menyuk et al. [17].

 

Fig. 2. Causal flow chart in TMI simulation.

Download Full Size | PDF

The heat generated at any given time, $Q = Q_0 + (1/2) [Q_+ e^{{i} \Delta \beta z} + Q_- e^{-{i} \Delta \beta z}]$, can be calculated as a function of the pump and signal intensities. This heat generation determines the evolution of the temperature, $\Delta T = \Delta T_0 + (1/2) [\Delta T_+ e^{{i} \Delta \beta z} + \Delta T_- e^{-{i} \Delta \beta z}]$. The $Q_{\pm }$ and $\Delta T_{\pm }$ terms represent the oscillatory heat and temperature patterns that arise as a result of intermodal beating. Both heat and temperature are real functions, meaning that $Q_+ = Q_-^*$ and $\Delta T_+ = \Delta T_-^*$. Since the heat diffusion equation is linear, the slowly varying temperature fields can be calculated by splitting the heat diffusion equation into two $\Delta \beta z$ harmonic components and solving for $\Delta T_0$ and $\Delta T_-$ separately [17].

Rather than solving these equations explicitly in the time domain, the temperature evolution can be expressed using the TE decomposition, and it is then solved using simple $e^{-\gamma _{mn} t}$ decay. So, substituting the expansion in Eq. (6) for either slowly varying temperature field in (12):

Equation of the form of (14) and (15) provide a simple expression for the evolution of the TE coefficients, but an expansion in the form of (6) still requires an infinite number of coefficients. This is obviously not computationally feasible, but many of the $\{\Delta T_{ab}^{c,s} \}$ terms make no contribution to the $c_{ij}$ coupling coefficients in the CMT equations. This can be seen by plugging (6) into (10) and (11) and using the thermo-optic index perturbation given by $\eta _T \Delta T$, where $\eta _T$ is the thermo-optic coefficient. Consider the thermo-optic contribution to the $c_{10}$ coefficient.

Recalling orthogonality relations of trigonometric functions, the thermo-optic contribution to the $c_{10}$ coefficient can be simplified by dropping all but the ${a}_{1n}$ terms in Eq. (6).

One of the major computational advantages of the TE method is that the transverse integrals shown above can be calculated outside the simulation loop. Simulation methods that use spatially resolved temperature distributions are required to compute the transverse integrals shown in (10) and (11) across the entire fiber cross section at every grid point along the fiber length and at every time step. In the TE method, the only transverse integrals that must be computed in the simulation loop are described by (15), and $Q_j$ is generally nonzero only in the actively doped fiber core, which is a much smaller area than the entire fiber cross section.

Although there are still infinite sums in (18) and (19), they can be accurately approximated by the first few dozen terms. Because the $\gamma _{mn}$ decay rates are calculated from progressively higher zeros of Eq. (5), the $\gamma _{mn}$ get larger with increasing $n$, corresponding to a smaller thermal decay time constants for higher-order TEs. In standard FDTD heat diffusion calculations, simulations often have a time step proportional to the spatial grid spacing, which must be made fine enough to represent the smallest relevant features of the temperature distribution. A reasonable approach would be to drop any terms in Eqs. (18) and (19) with time constants lower than the simulation time step of the corresponding FDTD simulation. This sum truncation can be quantitatively justified by calculating $k^{ab}_{mn} a_{mn}^j$ for a range of maximum $n$ and truncating when the product is vanishingly small. In practice, it is simplest to truncate to the smallest number of TE terms that gives simulation convergence for a given fiber design.

It is also worth noting that Eq. (4) shows continuous, analytical expressions for the TEs. Since the heat generation and decomposition is handled on a discrete spatial grid, the TE must also be evaluated on the same grid in order to compute the integrals in (15). Although the TEs are orthogonal in their analytical form, they will no longer be exactly orthogonal in their discretized form. To maintain the accuracy of the thermal decomposition, the discretized TEs can be re-orthogonalized using a Gram-Schmidt process [21].

4. Computational performance

While it is difficult to claim that the model described above is computationally superior to every possible model that uses spatially resolved FDTD for the heat diffusion, comparison of convergence and runtime between the TE model and a basic thermal FDTD model should be able to give a fairly accurate assessment of the TE model’s performance. The runtime and convergence tests of this paper used the fiber design and physical parameters shown in Table 1. The entirety of the high index core is actively doped, and $N_{Yb}$ is the active doping density. The $\sigma _{i}^{(j)}$ terms represent the emission and absorption cross sections of the $\text {Yb}^{3+}$ ions at the pump and signal wavelengths.

Table 1. Fiber design and simulation parameters

View Table

For both the thermal FDTD and TE methods, the fiber region was discretized with a 1cm spacing in the z-direction. Given that the intermodal beat length of this fiber is $\sim$2.55cm, a 1cm spacing roughly matches the resolution requirement found by Menyuk et al. [17]. Also, previous simulations using the TE model with a range of longitudinal resolutions showed that the simulated TMI threshold converges to a stable value for grid spacing below $\sim$5cm but falls off rapidly for grid spacing above $\sim$7.5cm, indicating that a 1cm resolution is sufficient for accurate simulation.

Time step resolution and transverse spatial resolution were varied in the simulation performance tests. The transverse spatial resolution was defined using a cartesian grid with the same grid spacing in x and y. For simplicity, these simulations all assume a well-cooled fiber boundary (equivalent to $h_c \rightarrow \infty$), giving a Dirichlet boundary condition with the outer surface of the fiber cladding held at room temperature. The number of eigenfunctions used in the TE method is represented by the minimum relevant thermal decay time constant, $\tau _{min}$. The thermal time constant is $\tau _{ab} = 1/\gamma _{ab}$. For a simulation with a given $\tau _{min}$, the TEs with $\tau _{ab} > \tau _{min}$ are used to model the thermal evolution, while TEs with smaller time constants are ignored. Figure 3 shows the number of eigenfunctions used by the TE model as a function of $\tau _{min}$.

 

Fig. 3. Number of eigenfunctions used by the TE model as a function of $\tau _{min}$ for the fiber described above.

Download Full Size | PDF

Before running any simulations, it is already possible to see one of the advantages of the TE method. Because the TE method ignores irrelevant data, it requires much less memory than the full spatially resolved model. Using the discretized fiber model described above and assuming a $2 \mu \text {m} \times 2 \mu \text {m}$ transverse spatial resolution, the thermal FDTD uses $\sim$2500 transverse grid points ($\sim$10,000 if not utilizing the symmetry of $LP_{01}$ and $LP_{11}$ coupling). This is multiplied by $\sim$1000 grid points in the z-direction and 8 bytes for each double-precision floating point temperature value, giving a memory requirement on the order of 20MB (note that this is doubled for the phase-matched model, which requires two separate temperature fields). For a single time step, this is a fairly modest memory requirement, but, at times, it may be helpful to examine and analyze the full evolution of the thermo-optic system, which would require large amounts of storage. For instance, saving the spatially and temporally resolved temperature distribution for $10\text {ms}$ of simulation time with $1 \mu \text {s}$ time steps would require $\sim$200GB of storage, which is by no means impossible but does make it difficult to store and analyze the thermal evolution. The TE method drastically reduces the storage required to record the full thermo-optic evolution. Simulation of the fiber described above with a $\tau _{min}$ of $1 \mu \text {s}$ uses 118 TEs, meaning that the simulation needs only $\sim$1MB in memory, and storage of the full thermo-optic evolution as described above requires only $\sim$10GB.

Understanding runtime advantage using TEs is a little less straightforward. A thermal FDTD simulation method must directly solve Eq. (12) discretely and then evaluate the integrals in (10) and (11), while a simulation using TEs must solve Eqs. (14) and (15) for each TE, requiring a larger number of numerical integrals but only over the area of the actively doped core. A theoretical comparison would be clunky and would be unlikely to be accurate for the full range of possible implementations of the FDTD heat diffusion. This work chooses to test the TE method by a direct comparison to a fairly standard thermal FDTD model that evaluates the spatially resolved temperature evolution using an alternating direction implicit (ADI) implementation of Eq. (12) [22], which offers high speed and numerical stability. One of the primary purposes of this work is to improve TMI simulation efficiency so that it becomes practical to model TMI using modest computational resources, so these simulations were run on hardware that is roughly equivalent to a consumer grade personal computer. For consistency, each simulation was run using four cores on an AMD 7763 processor with 16GB of allocated RAM.

Since a primary purpose of using time-domain TMI simulation (as opposed to frequency-domain or steady-periodic methods) is modelling of transient behavior and/or high-frequency modulation, the simulation convergence was evaluated using a simple, transient coupling. Each of the simulations was initiated with a previously calculated steady state temperature distribution for a 75W co-propagating pump and a pure FM 15W seed, meaning that there is no mode coupling in the initial state. This distribution was used as the starting point for the thermal FDTD simulations and was decomposed into the appropriate number of TE coefficients for the eigenfunction simulations. At the beginning of the simulation, $1{\% }$ of the seed power was transferred from the FM to the LP$_{11}$ HOM. An example of the resulting evolution of the output HOM power is shown in Fig. 4.

 

Fig. 4. Example of transient mode coupling data; the FDTD simulation uses a time step of $0.1 \mu \text {s}$ and transverse spatial resolution of $1 \mu \text {m}$ and the TE simulation uses time step of $16 \mu \text {s}$, transverse resolution of $9.5 \mu \text {m}$, and $\tau _{min} = 8 \mu \text {s}$.

Download Full Size | PDF

At these power levels, the system is well below the TMI threshold, but even a small transfer of power from the FM to the HOM moves some of the heat load from near the center of the doped core to the outer regions of the doped core, where it overlaps with the LP$_{11}$ mode. This causes an intermodal phase modulation, producing the phase shift between the modal interference pattern and the thermo-optic refractive index grating that is needed for power transfer from the FM to the HOM. This produces a transient power coupling from the FM to the HOM even though the system is below the TMI threshold. The experimental work in [23] showed a similar coupling below the TMI threshold. The two-mode TMI model cannot claim to be accurate if the FM becomes depleted. Since the transient coupling can be quite strong and it increases significantly with increased pump power, 75W of pump power strikes a good balance, showing nontrivial coupling behavior and maintaining FM dominance.

The coupling was computed for 10ms of simulation time with varying transverse spatial resolution, time step size, and $\tau _{min}$ values. This work uses the root mean squared error of the fraction of output signal power in the HOM to quantify convergence.

 

Fig. 5. RMSE convergence data for TE simulation with varying time step and $\tau _{min}$ using a TE simulation with $\tau _{min}$ of 0.1$\mu$s and time step of 0.1$\mu$s as the benchmark.

Download Full Size | PDF

The fraction of output power in the HOM will often be the most important result of a TMI simulation. It directly impacts the $\text {M}^2$ factor of the output beam, and it is often used in definitions of the TMI threshold. Also, for the fiber simulated in this work, the total output power varies by only a fraction of a watt over the 10ms of simulation time, meaning that $\zeta$ captures nearly all relevant information about the fiber output. For reference, the simulations depicted in Fig. 4 have $\zeta$ ranging from about $0.7{\% }$ to $15{\% }$.

Before checking the performance of the TE simulation method against the thermal FDTD method, it is important to determine an appropriate truncation for Eqs. (18) and (19), which in this work is represented by the choice of $\tau _{min}$. This should be chosen to provide a balance of simulation speed and accuracy. Figure 5 shows the accuracy of the TE simulation method for a range of time step and $\tau _{min}$ values. These simulations use a transverse resolution of $1 \mu \text {m}$. Equation (21) was evaluated using a TE simulation with $\tau _{min}$ of $0.1 \mu \text {s}$ and time step of $0.1 \mu \text {s}$ as the benchmark data. For reference, the minimum time constant values for these results correspond to 380 TEs at $\tau _{min} = 100\text {ns}$ but only 8 TEs at $\tau _{min} = 128\mu \text {s}$.

Observing Fig. 5, it is quite easy to see that the convergence error increases with increasing time step. The RMSE is roughly constant for $\tau _{min}$ below $16 \mu \text {s}$. For $\tau _{min}>32 \mu \text {s}$ the error increases rapidly. This justifies a truncation of the sums in Eqs. (18) and (19), and it implies that mode coupling dynamics are dominated by eigenfunctions with decay times greater than $16 \mu \text {s}$. Using the information in Fig. 5, a conservative $\tau _{min}$ of $8 \mu \text {s}$ is chosen for subsequent performance tests, which corresponds to 40 TE terms total and only 13 TEs with azimuthal order $m=1$, which is the only order capable of coupling power between the LP$_{01}$ and LP$_{11}$ modes. Simulations with this choice for $\tau _{min}$ run almost twice as fast as simulations using $\tau _{min}=100 \text {ns}$. Any errors introduced by truncating at $\tau _{min}=8 \mu \text {s}$ are considered negligible. Differences between the calculated RMSEs of simulations using $\tau _{min}=100 \text {ns}$ and simulations using $\tau _{min}=8 \mu \text {s}$ are on the order of $10^{-8}$, seven orders of magnitude lower than the range of $\zeta$ shown in Fig. 4.

 

Fig. 6. RMSE error for both FDTD and TE simulations using an FDTD simulation with 0.5$\mu$s time step and 0.5$\mu$m transverse spatial resolution as the benchmark.

Download Full Size | PDF

Using the appropriate choice of $\tau _{min}=8 \mu \text {s}$, the TE model can be compared directly to the FDTD thermal diffusion model. Runtime and convergence for both simulation methods were compared on a range of spatial and temporal resolutions, and the results are shown in Figs. 6(a-b). For both plots, Eq. (21) is calculated using reference data from an FDTD simulation with 0.5$\mu$s time step and 0.5$\mu$m transverse spatial resolution. The RMSE of the benchmark simulation is zero by definition, so it is not included as a data point in Fig. 6(a). It should also be noted that Fig. 6(a) is lacking convergence data for the highest resolution FDTD simulations because their runtime is impractically long for most applications. The "missing" data correspond to high-resolution simulations with transverse grid spacing below 1$\mu$m and temporal resolution below 2$\mu$s. Most prior work in simulating fiber with comparable design parameters uses transverse grid spacing of at least 2$\mu$m [17,24].

As expected, the convergence data clearly show that the accuracy of both models decreases with decreasing resolution. Visually comparing 6(a) and 6(b) it is easy to see that the convergence behavior of the two simulation methods is similar throughout most of the parameter sweep, with the TE method having notably superior convergence for the lowest resolution simulations. This does make intuitive sense. In the TE model, poor spatial resolution only impacts the evaluation of Eq. (15) but does not change TE decay rates. Equation (14) has exact, analytical solutions even for long time steps, limiting accumulated time step error. The FDTD model, on the other hand, solves (12) directly. It uses a discrete scheme that is only second-order accurate in both time and space, introducing error that grows rapidly with decreased resolution [22]. Over the whole simulation parameter sweep, some regions show superior TE convergence, but others show slightly superior FDTD convergence. Most practical applications should use simulation parameters that give high accuracy. Within the relatively accurate region, with $\Delta t<8 \mu \text {s}$ and $\Delta x = \Delta y < 4 \mu \text {m}$, the TE simulations have slightly lower convergence error, with the TE RMSE being on average $\sim$0.0215${\% }$ lower than the FDTD RMSE. This shows a small but nonzero accuracy advantage when using the TE method.

As a whole, any convergence advantage from using the TE method is quite small, but Fig. 6 does show that the TE method does not lose accuracy relative to the FDTD method. Since the two models being compared are simulating the same physics, it is expected that they would have fairly similar accuracy results. The primary advantage of the TE method is not its slightly superior convergence but its significantly reduced runtime. Measured runtime data is shown in Fig. 7(a-b). Each runtime data point is averaged over five simulation runs to partially account for runtime variance. Note that the x and y axes of these figures are flipped relative to Figs. 6(a-b) for visual clarity.

 

Fig. 7. Averaged simulation runtimes in minutes.

Download Full Size | PDF

Observing Figs. 7(a-b) it is clear that the TE method can give fast time-domain TMI evolution. For all of the parameter sweep data points that are defined in Fig. 7(a), which have runtimes up to 86 hours, the TE method gives runtimes below 6 hours. All the FDTD data points omitted from Figs. 6(a) and 7(a), had runtimes greater than four days and were considered impractical for most applications of TMI simulation.

Averaging over the same "relatively accurate" region used above ($\Delta t <8 \mu \text {s}$ and $\Delta x < 4 \mu \text {m}$), the TE method runs on average $\sim$13.9 times faster than the thermal FDTD method. Because most prior literature is limited to fairly low spatial and temporal resolution, it is worth comparing a few specific points in the parameter sweep. Using the TE method for a relatively high-resolution simulation with $\Delta t =1 \mu \text {s}$ and $\Delta x =1 \mu \text {m}$ gives a factor of $\sim$10.2 runtime improvement, for $\Delta t =2 \mu \text {s}$ and $\Delta x =2 \mu \text {m}$ TEs reduce runtime by a factor $\sim$7.1, and for a low-resolution simulation with $\Delta t =8 \mu \text {s}$ and $\Delta x =6 \mu \text {m}$ TEs reduce runtime by a factor of $\sim$4.0.

It is also worth noting that for much of the simulation parameter sweep in Fig. 7(b) the TE simulation method can calculate 10ms of thermo-optic fiber amplifier evolution in significantly less than an hour. Obviously, runtime will depend on many factors (simulation resolution, length of time being simulated, computer hardware, and fiber design), but, for many practical applications, this simulation method can calculate time-domain TMI evolution in less than 30 minutes. For example, the simulation runs with $\Delta t =2 \mu \text {s}$ and $\Delta x =2 \mu \text {m}$, which used computational resources roughly analogous to those found in the average consumer PC, took a little under 12 minutes on average. This type of rapid simulation can easily fit into an engineering workflow with iterative design adjustment followed by computational validation.

5. Conclusion

This work implemented a TMI simulation model for $\text {Yb}^{3+}$ doped fiber amplifiers that uses the temperature eigenfunctions of the transverse heat diffusion equation to simplify the calculation of both the time evolution and the coupled-mode propagation. Compared to prior work that used fiber thermal modes, this model drops thermal expansion terms that cannot contribute to the two-mode CMT propagation. It accurately represents the time-domain dynamics of TMI, offering comparable accuracy relative to the standard, spatially resolved FDTD thermal diffusion method. By dropping unnecessary terms in the thermal expansion, this method eliminates computation and storage of unnecessary data, giving it gave a runtime reduction by an average factor of $\sim$13.9 for simulations with fairly high accuracy. It allows TMI evolution in many realistic fiber systems to be evaluated in less than an hour on a standard, low-performance personal computer.