Simulation of Optical Fiber Amplifier Gain Using Equivalent Short Fibers

TL;DR

This paper introduces an efficient modeling approach for optical fiber amplifiers using equivalent short fibers, enabling faster simulations while maintaining essential characteristics, validated through numerical studies on doped fibers.

Contribution

The paper proposes a novel scale model called the equivalent short fiber that simplifies Maxwell-based simulations of optical amplifiers, improving computational efficiency.

Findings

The equivalent short fiber model accelerates computations by a factor proportional to the length reduction.

Numerical validation shows the model's effectiveness for specific fiber types.

The model's applicability depends on fiber properties and desired accuracy.

Abstract

Electromagnetic wave propagation in optical fiber amplifiers obeys Maxwell equations. Using coupled mode theory, the full Maxwell system within an optical fiber amplifier is reduced to a simpler model. The simpler model is made more efficient through a new scale model, referred to as an equivalent short fiber, which captures some of the essential characteristics of a longer fiber. The equivalent short fiber can be viewed as a fiber made using artificial (unphysical) material properties that in some sense compensates for its reduced length. The computations can be accelerated by a factor approximately equal to the ratio of the original length to the reduced length of the equivalent fiber. Computations using models of two commercially available fibers -- one doped with ytterbium, and the other with thulium -- show the practical utility of the concept. Extensive numerical studies are conducted to assess when the equivalent short fiber model is useful and when it is not.

Tables

Table 1. Table 1. Parameters used in Yb-doped fiber simulation

Parameter	Value	Units	Parameter	Value	Units
${\lambda }_p=2\pi c/{\omega }_p$	$976\text{×}{10}^{-9}$	m	${\lambda }_s=2\pi c/{\omega }_s$	$1064\text{×}{10}^{-9}$	m
${\sigma }^{\text{abs}}\left({\omega }_p\right)$	$1.429\text{×}{10}^{-24}$	${\text{m}}^2$	${\sigma }^{\text{ems}}\left({\omega }_p\right)$	$1.776\text{×}{10}^{-24}$	${\text{m}}^2$
${\sigma }^{\text{abs}}\left({\omega }_s\right)$	$6\text{×}{10}^{-27}$	${\text{m}}^2$	${\sigma }^{\text{ems}}\left({\omega }_s\right)$	$3.58\text{×}{10}^{-25}$	${\text{m}}^2$
$N_{\text{total}}$	$3\text{×}{10}^{26}$	${\text{ions/m}}^3$	$\tau$	$8\text{×}{10}^{-4}$	s
$n_{\text{core}}$	$\mathrm{1.450\hspace{0.17em}971}$	–	NA	0.06	–
$r_{\text{core}}$	$1.25\text{×}{10}^{-5}$	m	$r_{\text{clad}}$	$2\text{×}{10}^{-4}$	m
$P_p^0$	1000	W	$P_s^0$	25	W

Table 2. Table 2. Parameters used in Tm-doped fiber simulation

Parameter	Value	Units	Parameter	Value	Units
${\lambda }_p=2\pi c/{\omega }_p$	$793\text{×}{10}^{-9}$	m	${\lambda }_s=2\pi c/{\omega }_s$	$2110\text{×}{10}^{-9}$	m
${\sigma }^{\text{abs}}\left({\omega }_p\right)$	$4.4686\text{×}{10}^{-25}$	${\text{m}}^2$	${\sigma }^{\text{ems}}\left({\omega }_p\right)$	0	${\text{m}}^2$
${\sigma }^{\text{abs}}\left({\omega }_s\right)$	$1.7423\text{×}{10}^{-27}$	${\text{m}}^2$	${\sigma }^{\text{ems}}\left({\omega }_s\right)$	$\mathrm{1.173\hspace{0.17em}97}\text{×}{10}^{-25}$	${\text{m}}^2$
${\tau }_{10}$	$6.2232\text{×}{10}^{-03}$	s	${\tau }_{20}$	$5.5179\text{×}{10}^{-03}$	s
${\tau }_{21}$	$2.5707\text{×}{10}^{-01}$	s	${\tau }_{30}$	$1.3949\text{×}{10}^{-03}$	s
${\tau }_{31}$	$1.7033\text{×}{10}^{-02}$	s	${\tau }_{32}$	$6.8446\text{×}{10}^{-02}$	s
${\Gamma }_1$	$\mathrm{2.592\hspace{0.17em}88}\text{×}{10}^{03}$	Hz	${\Gamma }_2$	$\mathrm{2.927\hspace{0.17em}55}\text{×}{10}^{07}$	Hz
${\Gamma }_3$	$\mathrm{8.059\hspace{0.17em}43}\text{×}{10}^{04}$	Hz	–	–	–
$N_{\text{total}}$	$3\text{×}{10}^{26}$	${\text{ions/m}}^3$	${\kappa }_R$	$1.17\text{×}{10}^{-21}$	${\text{m}}^3$
$n_{\text{core}}$	$\mathrm{1.439\hspace{0.17em}994}$	–	NA	0.1	–
$r_{\text{core}}$	$1.25\text{×}{10}^{-5}$	m	$r_{\text{clad}}$	$2\text{×}{10}^{-4}$	m
$P_p^0$	1100	W	$P_s^0$	30	W

Equations

$${\vec{\mathcal{E}}}_{\ell }(x,y,z,t)$$

$$\begin{array}{rl}\text{curl}{\vec{E}}_{\ell }-\hat{ı}{\omega }_{\ell }{\mu }_0{\vec{H}}_{\ell }& =0,\\ \text{curl}{\vec{H}}_{\ell }+\hat{ı}{\omega }_{\ell }{ϵ}_0{\vec{E}}_{\ell }& =-\hat{ı}{\omega }_{\ell }{\vec{P}}_{\ell },\end{array}\ell \in \{s,p\},$$

$${\vec{P}}_{\ell }^{\text{bg}}({\vec{E}}_{\ell })$$

$$\text{curl}\text{curl}{\vec{E}}_{\ell }-{\omega }_{\ell }^2{ϵ}_0{\mu }_0{\vec{E}}_{\ell }={\omega }_{\ell }^2{\mu }_0{\vec{P}}_{\ell }$$

$${\vec{P}}_{\ell }={\vec{P}}_{\ell }^{\text{bg}}+{\vec{P}}_{\ell }^{\text{ag}}={ϵ}_0(n^2-1){\vec{E}}_{\ell }-\frac{\hat{ı}{ϵ}_{0}cn}{{\omega }_{\ell }}{g}_{\ell }{\vec{E}}_{\ell }$$

$$\text{curl}\text{curl}{\vec{E}}_{\ell }-k_{\ell }^2{n}^2{\vec{E}}_{\ell }+\hat{ı}{k}_{\ell }n{g}_{\ell }{\vec{E}}_{\ell }=0,$$

$${\vec{E}}_{\ell }(x,y,z)=U_{\ell }(x,y,z){\hat{e}}_x,$$

$$-\Delta U_{\ell }-k_{\ell }^2{n}^2{U}_{\ell }+\hat{ı}{k}_{\ell }n{g}_{\ell }{U}_{\ell }=0.$$

$$({\Delta }_{xy}+k_{\ell }^2{n}^2){\phi }_m={\beta }_m^2{\phi }_m,m=1,\dots ,M_{\ell },$$

$$U_{\ell }(x,y,z)=\sum _{m=1}^{M_{\ell }}{A}_m^{(\ell )}(z){\phi }_m(x,y)e^{\hat{ı}{\beta }_{m}z},$$

$$\sum _{m=1}^{M_{\ell }}\frac{d{A}_m^{(\ell )}}{dz}{\phi }_m{\beta }_m{e}^{\hat{ı}{\beta }_{m}z}$$

$$\frac{d{A}_l^{(\ell )}}{dz}$$

$$K_{lm}^{(\ell )}(I_{\ell },I_{{\ell }^c})=\frac{k_{\ell }}{2{\beta }_l}{\int }_{{\Omega }_z}{g}_{\ell }(I_{\ell },I_{{\ell }^c})n(x,y){\phi }_m(x,y)\overline{{\phi }_l(x,y)}dxdy,$$

$$\frac{d{A}_1^{(p)}}{dz}=K_{11}^p(I_p,I_s)A_1^{(p)},$$

$$K_{11}^p(I_p,I_s)=\frac{1}{2\pi r_{\text{clad}}^2{n}_{\text{clad}}}{\int }_{{\Omega }_z}{g}_p(I_p,I_s)n_{\text{clad}}dxdy=\frac{1}{2}⟨g_p⟩.$$

$$I_{\ell }(x,y,z)=\frac{n}{{\mu }_{0}c}{\left|\sum _{m=1}^{M_{\ell }}{A}_m^{(\ell )}(z)e^{\hat{ı}{\beta }_{m}z}{\phi }_m(x,y)\right|}^2.$$

$$I_p(z)=\frac{n}{{\mu }_{0}c\pi r_{\text{clad}}^2}{\left|A_1^{(p)}(z)\right|}^2.$$

$$\frac{d{I}_p}{dz}=⟨g_p⟩I_p,$$

$$I_s\equiv I_s(x,y,z,A)=\frac{n}{{\mu }_{0}c}{\left|\sum _{m=1}^M{A}_m(z)e^{\hat{ı}{\beta }_{m}z}{\phi }_m(x,y)\right|}^2.$$

$$Y(z)=[I_p(z),A_1(z),A_2(z),\dots ,A_M(z){]}^t,0<z<L,$$

$$\frac{dY}{dz}$$

$$Y(0)$$

$$N_{\text{total}}=N_0(x,y,z,t)+N_1(x,y,z,t)+N_2(x,y,z,t)+N_3(x,y,z,t)$$

$${\nu }_{\ell }=\frac{I_{\ell }}{\hslash {\omega }_{\ell }},\ell \in \{s,p\}$$

$${\partial }_t{N}_3$$

$${\partial }_t{N}_2$$

$${\partial }_t{N}_1$$

$$N_{\text{total}}$$

$${\psi }_{\ell }^{\text{abs}}={\sigma }^{\text{abs}}({\omega }_{\ell }){\nu }_{\ell },{\psi }_{\ell }^{\text{ems}}={\sigma }^{\text{ems}}({\omega }_{\ell }){\nu }_{\ell },\ell \in \{s,p\}.$$

$${\delta }_i=\sum _{j=0}^{i-1}{\tau }_{ij}+{\Gamma }_i,{\gamma }_0=\frac{1}{{\psi }_p^{\text{ems}}+{\delta }_3}{\gamma }_1={\psi }_p^{\text{abs}}{\gamma }_0,{\gamma }_2=\frac{{\tau }_{32}^{-1}+{\Gamma }_3}{{\delta }_2},$$

$${\gamma }_3=\frac{{\tau }_{31}^{-1}+{\gamma }_2({\tau }_{21}^{-1}+{\Gamma }_2)+{\gamma }_1^{-1}{\psi }_s^{\text{abs}}}{{\psi }_s^{\text{ems}}+{\delta }_1},{\gamma }_4=\frac{2{\psi }_p^{\text{abs}}+{\psi }_s^{\text{abs}}}{{\psi }_p^{\text{abs}}({\psi }_s^{\text{ems}}+{\delta }_1)}.$$

$$N_0$$

$$+\frac{\sqrt{(1-{\gamma }_0{\kappa }_R{N}_{\text{total}}+{\gamma }_1(1+{\gamma }_2+{\gamma }_3){)}^2+4({\gamma }_0+{\gamma }_1{\gamma }_4){\kappa }_R{N}_{\text{total}}}}{2{\kappa }_R({\gamma }_0+{\gamma }_1{\gamma }_4)},$$

Full text

This work was supported in part by AFOSR grant FA9550-17-1-0090.

1. Introduction

“Scale models” are ubiquitous in fields such as fluid dynamics. They are physical or numerical models that preserve some of the important properties of an object being modeled while not preserving the original dimensions of the object. The main goal of this paper is to formulate and study a miniature scale model of an optical fiber laser amplifier. Our scale model reduces fiber length to increase computational efficiency. While unable to preserve all properties of the original electromagnetic solution, our numerical scale model is able to approximately replicate the original fiber’s power distribution, as we shall see in later sections. After this introductory section, we will begin by describing a simplified model of beam propagation in fibers. This model will then be used to derive, justify, and verify the scale model.

The importance of fiber amplifiers in enabling our current world of long-distance fiber optics and submarine telecommunications cannot be overlooked [4]. High power fiber amplifiers also have many other uses, for example, as defensive speed-of-light weapons. High output powers have been achieved by solid-state optical fiber laser amplifiers [9]. Numerical modeling of these optical devices has also been effectively used by many [10, 12, 17, 19]. Yet, simulation of full length fibers remains cumbersome and far from being routine. This is because of the long simulation times and the large computational resources required. Simulations using the full Maxwell system are too expensive since there are millions of wavelengths within any realistically long fiber. An an example, consider the full Maxwell simulation of Raman gain attempted in [13]: more than five million degrees of freedom was needed to simulate an extremely short fiber containing 80 wavelengths (less than 0.0001 m). Although a full Maxwell model of a realistically long (10 m) fiber can be written out (and we shall do so in Subsection 2.1), its numerical solution is beyond the reach of today’s simulation capabilities.

Therefore, simplified models form the current state of the art. It is somewhat surprising how unreasonably effective these models have proved to be, despite the drastic simplifications used in their derivation. The state of the art in fiber amplifier simulation consists of beam propagation methods using coupled mode theory (CMT). We shall introduce the reader to the CMT model in Subsection 2.2 as a simplification of the full Maxwell model. To facilitate cross-disciplinary readership, we make an effort to enunciate the assumptions behind such simplifications. Even though it is not common in the optics literature to view CMT in the backdrop of emerging developments in reduced-order models, one may view it as essentially an example of a physics-based reduced-order model. Indeed, in CMT, the electromagnetic solution is expressed using a “reduced basis” consisting of transverse guided modes of the fiber that encapsulates the energy-transport mechanism in fibers.

Even the simplified CMT model is computationally too demanding to ably assist with the important open issues in the subject today. One of these issues is what is currently recognized to be the main roadblock to power scaling of beam combinable fiber amplifiers, namely the nonlinear transverse mode instability (TMI). TMI can be described as a sudden breakdown in beam quality at high power operation, first observed experimentally [5]. As pointed out in the review [9], when attempting to design highly coherent lasers capable of sustained high (average) powers, a practically uncrossable limit was encountered due to the TMI. After intensive speculations on the cause of TMI, the prevailing theory seems to be that the cause is a temperature-induced grating. We believe that numerical modeling is essential for investigating the TMI, and other nonlinearities that arise inside fiber amplifiers, since experimental evidence is mostly limited to examining the amplifier output, rather than the onset of physical effects that occur inside of the glass fiber along its length. The current difficulty in using numerical models is the excessive simulation times: indeed any numerical technique used must be able to solve for the electromagnetic field within a long fiber a vast number of times. Given the great computational burden of capturing length scales as small as 10 $\mu$m, and time scales as small as 10 $\mu$sec (for the thermal problem), techniques that further accelerate the numerical simulations have the potential to significantly enhance the ability for computer modeling to inform experimental designs and configurations in a timely manner. It is our intent to contribute such an acceleration technique by developing the above-mentioned scale model (in Sections 4–5). Studies of its application to TMI investigations are postponed to the future.

The models are tested using two commercially available examples of doped step-index fibers, one with ytterbium (Yb) doping in the fiber core, and another with thulium (Tm) doping. Both are examples of large mode area (LMA) fibers which support more than one guided transverse core mode. LMA fibers are of great interest since they permit greater light amplification per unit length and help mitigate the onset of other detrimental optical nonlinearities. Unfortunately, they are also more susceptible to the TMI, and hence stand most to benefit from advances in numerical simulation. Our active gain model for these fibers utilizes the population dynamics of Yb and Tm ions. Active gain in fiber amplifiers appears as a nonlinear coupling term between the Maxwell systems for the (1) less coherent “pump” light that supplies energy for amplification, and the (2) highly coherent “signal” (laser) light. The gain mechanism involves exciting the outer most electron of the dopant (Yb or Tm) by absorbing the pump light, and producing more coherent signal light via stimulated emission that allows the excited electrons to return to their ground state. We have included a simplified, yet very typical, mathematical formulation for the dopant ion population dynamics in Section 3. A few of the initial results obtained in this work for the simpler Yb-doped case were announced earlier in the conference proceedings [6]. We begin by deriving the CMT model next.

2. The CMT model

Physics-based reduced-order models are now being used successfully in various simulation techniques [18]. In this section, we introduce such a model for an optical fiber amplifier starting from Maxwell equations. We will start from the Maxwell system and describe the assumptions that lead us to the simplified CMT model consisting of a system of ordinary differential equations (ODE).

2.1. The full Maxwell model

Suppose the optical fiber amplifier to be modeled is aligned so that it is longitudinally centered along the $z$-axis; the transverse coordinates will be denoted by $x$ and $y$ in the Cartesian coordinate system. The core region of the fiber, $\{(x,y,z):x^{2}+y^{2}<r_{\text{core}}^{2}\}$, is enveloped by a cladding that extends to radius $r_{\text{clad}}$. The fiber is a step-index fiber, i.e., its refractive index is a piecewise constant function that takes the value $n_{\text{core}}$ in the core region and $n_{\text{clad}}$ in the cladding region. There is usually another polymer coating that surrounds this inner cladding (composed of fused silica); however, this second cladding/coating can readily be neglected for this analysis since the laser light is mostly guided in the fiber core region. We want to model a continuous wave, weakly guided ($n_{\text{core}}-n_{\text{clad}}\ll 1$), polarization maintaining, large mode area (LMA) fiber. There are various arrangements in which this fiber amplifier could be seeded and pumped. We consider the co-pumped/clad-pumped configuration, wherein a highly coherent laser light – which we shall refer to as the signal– is injected into the fiber core area at the beginning of the fiber ($z=0$). The pump light is injected into the fiber at $z=0,$ and unlike the signal, it enters both core and cladding.

Let $\vec{\mathcal{E}}_{s},\vec{\mathcal{H}}_{s}$ and $\vec{\mathcal{E}}_{p},\vec{\mathcal{H}}_{p}$ denote the electric and magnetic fields of the signal and pump light, respectively. The signal and pump fields are assumed to be time harmonic of frequencies $\omega_{s}$ and $\omega_{p}$ respectively, i.e.,

[TABLE]

Here and throughout, we use the subscript $\ell\in\{s,p\}$ to distinguish between signal and pump fields. Note that the $\vec{\mathcal{E}}_{\ell}$ and $\vec{\mathcal{H}}_{\ell}$ are real valued while $\vec{{E}}_{\ell}$ and $\vec{{H}}_{\ell}$ are complex valued. The signal field $\vec{{E}}_{s},\vec{{H}}_{s}$ and the pump field $\vec{{E}}_{p},\vec{{H}}_{p},$ are assumed to independently satisfy Maxwell equations, but are coupled through the electric polarization terms of the form $\vec{P}_{\ell}\equiv\vec{P}_{\ell}(\vec{{E}}_{s},\vec{{E}}_{p})$, $\ell\in\{s,p\}$, which appear in the following time-harmonic Maxwell system,

[TABLE]

where $\epsilon_{0}$ is the electric permittivity and $\mu_{0}$ is the vacuum magnetic permeability.

All interactions between the propagation medium and the electromagnetic field are modeled through electric polarization terms. The traditional polarization model includes linear susceptibility, namely the background material interaction $\vec{P}^{\text{bg}}_{\ell}$ given as a function of the index of refraction of the medium that the light propagates through. Other examples of polarization terms include those that account for linear loss, active laser gain ($\vec{P}^{\text{ag}}_{\ell}$), thermal effects, and optical nonlinearities such as Brillouin scattering, Raman scattering, and Kerr effects. Here we focus on active gain polarization and the linear background polarization, namely,

[TABLE]

where c is the speed of light and $g_{\ell}$ is the *active gain term that depends on $\vec{{E}}_{\ell}$ in some nonlinear * fashion. Examples of $g_{\ell}$ are given in Section 3. Typical optical operating frequencies imply that within a fiber of realistic length there are several millions of wavelengths. Even if a mesh fine enough to capture the wave oscillations is used, the pollution effect [3] in wave propagation simulations destroys the accuracy of finite element solutions at the end of millions of wavelengths. Hence, without further simplifications, the above-described full Maxwell model is not a feasible simulation tool for realistic fiber lengths. We proceed to develop a physics-based reduced model.

2.2. Coupled mode theory

Experiments indicate that the vast majority of the laser signal is contained within the guided core modes of the fiber, and, likewise, most of the pump light is within the guided cladding modes. This is the basis of an electric field ansatz that CMT uses. Before giving the ansatz, let us eliminate $\vec{{H}}_{\ell}$ from (2.1), to obtain the second order equation

[TABLE]

solely for the electric field. Substituting

[TABLE]

into (2.2), using $c=1/\sqrt{\epsilon_{0}\mu_{0}}$ and simplifying we get,

[TABLE]

where $k_{\ell}=\omega_{\ell}/c$ is the wavenumber corresponding to the frequency $\omega_{\ell}$.

Next, we assume that the electric field $\vec{{E}}_{\ell}$ can be expressed as

[TABLE]

i.e., it is linearly polarized in a fixed transverse direction, which is taken above to be the $x$-direction (where $\hat{e}_{x}$ denotes the unit vector in the $x$-direction). Furthermore, since $\vec{{E}}_{\ell}$ has high frequency oscillations along the $z$-direction, its variations along the transverse directions may be considered negligible. It is therefore standard in optics to neglect $\mathop{{\text{grad}}}\mathop{\text{missing}}{div}\vec{{E}}_{\ell}.$ These assumptions imply that the vector equation (2.4) becomes the following scalar Helmholtz equation for $U_{\ell}$,

[TABLE]

Due to the high wave number $k_{\ell}$, even this simplified scalar field problem is computationally intensive. We now proceed to further reduce this scalar model using CMT.

CMT is usually useful in the analysis of the interaction between several near-resonance guided modes. For step-index fiber waveguides these modes are called linearly polarized transverse guided core modes [2], often referred to simply as LP modes. Mathematically speaking, these modes are finitely many non-trivial functions $\varphi_{m}(x,y)$, $m=1,2,\ldots,M_{\ell}$, that decay exponentially at the edge of the cladding region and satisfy

[TABLE]

where $\beta_{m}$ is the corresponding propagation constant and $\Delta_{xy}=\partial_{xx}+\partial_{yy}$ denotes the transverse Laplacian operator. The CMT approach to solve (2.5) expresses the solution using the ansatz

[TABLE]

where $A^{({\ell})}_{m}(z)$ denotes the complex field amplitude of mode $m$. Therefore, the wavenumber ($k_{\ell}n$) for the entire electric field envelop ($U_{\ell}$) is now decomposed into individual propagation constants ($\beta_{m}$) corresponding to each guided mode, and the field envelop is now decomposed into parts of amplitudes $A^{({\ell})}_{m}$ having transverse profiles described by $\varphi_{m}$.

Knowledge of the form of the solution is thus incorporated a priori into the ansatz. In particular, the physical intuition that the $\varphi_{m}$-component should oscillate longitudinally at an approximate frequency of $\beta_{m}$ is built in. This justifies the next assumption that $A^{({\ell})}_{m}(z)$ is a slowly varying function of $z$ (having built the fast variations in $z$ into the $e^{\hat{\imath}\beta_{m}z}$ term). Accordingly, for each $A^{({\ell})}_{m}$, we neglect the second-order derivative $d^{2}A^{({\ell})}_{m}/dz^{2}$ for all $m=1,\ldots,M_{\ell}$. Doing so after substituting (2.7) into (2.5) and using (2.6) we obtain

[TABLE]

The next step is to multiply both sides of (2.8) by the complex conjugate of $\varphi_{l}$, namely $\overline{\varphi}_{l}$, and integrate. We integrate over $\varOmega_{z},$ which represents the fiber cross section having the constant longitudinal coordinate value of $z$. Then, simplifying using the $L^{2}(\varOmega_{z})$-orthogonality of the modes,

[TABLE]

for $l=1,\ldots,M_{\ell}$, where $K^{({\ell})}_{lm}$ is the mode coupling coefficient, given by

[TABLE]

$\ell^{c}\in\{s,p\}\setminus\{\ell\}$, and $I_{s},I_{p}$ denote the signal and pump irradiance, respectively, which are formulated later in this subsection.

For the pump light, the number of guided cladding modes is exceedingly large: $M_{p}>10^{5}$. Rather than modeling each of these modes, it is sufficient to approximate the pump field as a plane wave, which effectively acts as the composition of all of the pump guided modes [12, 17]. Accordingly, we set $M_{p}=1$ and the normalized mode $\varphi_{1}^{p}=(\sqrt{\pi}r_{\text{clad}})^{-1}$ (without a transverse dependence). Since the cladding region is many times larger than the core, the corresponding propagation constant is estimated as if this mode travels in a uniform medium of refractive index $n_{\text{clad}}$, i.e., we set $\beta_{1}=k_{p}n_{\text{clad}}=\omega_{p}n_{\text{clad}}/c$. Then (2.9) yields

[TABLE]

for $0<z<L$, where

[TABLE]

Here $\langle{{g_{p}}}\rangle=|\varOmega_{z}|^{-1}\int_{\varOmega_{z}}g_{p}\;dx\,dy$ denotes the mean value of $g_{p}$ taken over $\varOmega_{z},$ the area of the fiber cross section out to $r=r_{\text{clad}}$.

The irradiance is proportional to the square of the field envelop magnitude, $I_{\ell}=n|U_{\ell}|^{2}/(\mu_{0}c)$. Using (2.7),

[TABLE]

For the pump plane wave, this reduces to

[TABLE]

Using the equation (2.11) and its complex conjugate, elementary simplifications lead to the following governing ODE for the pump irradiance:

[TABLE]

In view of (2.14), instead of $A^{({p})}_{1}$, we shall use $I_{p}(z)$ as our pump unknown. There is no need for the amplitude $A^{({p})}_{1}$ in the remainder of the model. Hence from now on, we write $A_{m}$ for $A^{({s})}_{m}$ dropping the superscript. We shall also simply write $M$ for $M_{s}$ and $K_{lm}$ for $K^{({s})}_{lm}$.

Next, consider the signal irradiance, namely the $\ell=s$ case in (2.13). To highlight the dependence of $I_{s}$ on $A_{m}\equiv A^{({s})}_{m}$, we use $A\equiv[A_{1}(z),\ldots,A_{M}(z)]^{t}$ to collectively denote the set of all signal mode amplitudes and write

[TABLE]

Note that the modes $\varphi_{l}(x,y)$ and the propagation constants $\beta_{l}$ may be precomputed (and the cost of this precomputation corresponds to the “off-line” computational cost in this reduced-order model).

In order to complete the CMT model (assuming we have expressions for $g_{\ell}$), we need to provide initial conditions at $z=0$, the beginning of the fiber. What is usually known is the power contained in the pump and signal light. The initial pump irradiance $I_{p}^{0}=I_{p}(0)$ can be calculated from the initial pump power $P_{p}^{0}$ provided at the inlet in a co-pumped configuration, by $I_{p}^{0}=|\varOmega_{0}|^{-1}P_{p}^{0}$. We assume that we also know how the signal light is split into various modes at the inlet, i.e., we may set $A(0)$ to some given $A^{0}=[A_{1}^{0},\ldots,A_{M}^{0}]^{t}.$ In practice, most of the signal power is usually carried in the first fundamental mode.

To summarize, the CMT model computes

$Y(z)=[I_{p}(z),A_{1}(z),A_{2}(z),\ldots,A_{M}(z)]^{t},\qquad 0<z<L,$

where each $A_{m}(z)$ is a signal mode amplitude in the fiber core, and $Y(z)$ satisfies the ODE system

$\displaystyle\frac{dY}{dz}$ $\displaystyle=\begin{bmatrix}\langle g_{p}(Y)\rangle&0\\ 0&\phi(z)\cdot K(Y)\end{bmatrix}Y,$ $\displaystyle 0<z<L,$

(2.16a)

$\displaystyle Y(0)$ $\displaystyle=[I_{p}^{0},A^{0}]^{t}$ $\displaystyle z=0,$

(2.16b)

where $\phi(z)$ is an $M\times M$ matrix defined by $\phi_{lm}(z)=e^{\hat{\imath}(\beta_{m}-\beta_{l})z},$ $K(Y)$ is a matrix of the same size whose $(l,m)$th entry is $K^{({s})}_{lm}({I_{s}},I_{p})$ defined in (2.10), and $\phi(z)\cdot K(Y)$ denotes the Hadamard product of $\phi$ and $K$, i.e., $[\phi\cdot K]_{lm}=\phi_{lm}K_{lm}$.

3. Thulium and ytterbium doped fiber amplifiers

Thulium (Tm)-doped fiber amplifiers [7, 8] can operate in eye-safe laser wavelengths (larger than 1.4 $\mu$m) and can reach an atmospheric transmission window (2.1–2.2 $\mu$m). There are efficient high-power LEDs that operate in the range of 0.79-0.793 $\mu$m, which is a peak absorption bandwidth for Tm-doped fibers. Cross-relaxations and upconversions occur in Tm-doped amplifiers. Even though Tm-doped fibers usually have better TMI suppression compared to other rare-earth ion doped fibers [17], ytterbium (Yb)-doped fiber amplifiers have also emerged as excellent candidates for high power operation due to their high-efficiencies and low amplified spontaneous emission gain. Yb-doped amplifiers are usually pumped at 976 nm and can lase around 1064 nm very efficiently. The dynamics of both these ion populations are explained below. They complete our model by giving expressions for $g_{\ell}$ to be used in (2.16).

3.1. Tm-dopant ion dynamics

The Tm ion population dynamics are schematically represented in Figure 1. The model involves four manifolds. The total number of Tm ions (per volume) is

[TABLE]

where $N_{0}$ represents the ground state (manifold 0) ion-population concentration, while $N_{1}$, $N_{2},$ and $N_{3}$ denote ion concentrations at excitation manifolds 1,2 and 3, respectively. What we have named energy manifolds 0,1,2, and 3, represent Tm energy levels usually written as ${}^{3}H_{6}$ ${}^{3}F_{4}$, ${}^{3}H_{5}$ and ${}^{3}H_{4}$, respectively.

Pump light of frequency $\omega_{p}=793$ nm excites the Tm ground state ions into higher energy manifolds, thus depleting manifold 0 at the rate $\nu_{p}\sigma^{\text{abs}}(\omega_{p})N_{0}$ while increasing the excited manifold $j$ at the rate $\nu_{p}\sigma^{\text{ems}}(\omega_{p})N_{j}$, where $\sigma^{\text{abs}}$ and $\sigma^{\text{ems}}$ represent measurable absorption and emission cross sections of Tm [1], and

[TABLE]

represents the flux of photons of frequency $\omega_{\ell}$. We must also take into account the fact that an excited ion in manifold $j$ can decay spontaneously to a lower energy manifold $k$ at the rate $1/\tau_{jk}$. An excited ion in manifold $j$ can also decay non-radiatively to the next lower energy manifold at the rate $\Gamma_{j}$. Finally, an excited Tm ion can also undergo cross-relaxation, wherein it transfers part of its energy to a ground state ion so both can end up in an intermediate energy level. Cross-relaxation is represented by the slanted arrows in Figure 1, while the other processes are represented by up/down arrows. The rate constant for the cross-relaxation is denoted by $\kappa_{R}$. Cross-relaxation, which creates two excited Tm ions for every pump photon (a two-for-one process), increase the amplifier efficiency (while upconversions, which are neglected in our model, decrease fiber efficiency). Following [11], these processes are modeled by

[TABLE]

where

[TABLE]

In our simulations, we have set $\omega_{s}$ to correspond to signal light of wavelength 2100 nm.

Next, we make the simplifying assumption that all the time derivatives $\partial_{t}$ in (3.2) may be neglected. By doing so, we are neglecting the time variations in the ion populations that occur at an extremely small time scale of around $10^{-5}$ s. Equations (3.2a)–(3.2c) after setting $\partial_{t}=0$ immediately yield $N_{1},N_{2},N_{3}$ in terms of $N_{0}$. The last equation (3.2d) then gives a quadratic equation for $N_{0}$. To express this solution, first define

[TABLE]

Then, the steady-state solution is given explicitly by

[TABLE]

Using this, we set the gain expressions by

[TABLE]

This completes the prescription of the CMT model (2.16) for Tm-doped fiber amplifier.

3.2. Yb-dopant ion dynamics

The model for population dynamics of Yb ions is simpler as it can be modeled using only two energy states, the ground state and one excited state manifold, as shown in Figure 2. Hence, instead of (3.1), we now have

[TABLE]

where $N_{\text{total}}$ denotes the total population concentration in the fiber, $N_{\text{ground}}$ represents the ground state ion-population (in $\mathrm{{}^{2}F_{7/2}}$) and $N_{\text{excited}}$ denotes the excited state ion-population (in $\mathrm{{}^{2}F_{5/2}}$). The absorption and emission that models the two-state dynamics now result in

[TABLE]

where now we must use the absorption and emission cross section values [14] of Yb for $\sigma^{\text{abs}},\sigma^{\text{ems}}$ while computing $\psi^{\text{abs}}_{\ell},\psi^{\text{ems}}_{\ell}$. The parameter $\tau$ is the upper level radiative lifetime of the excited state. As in the Tm case, we assume that the system has already reached the steady-state solution. Putting the time derivative in (3.6) to zero, a simple calculation shows that

[TABLE]

Finally, the active gain expressions are modeled in terms of the above $N_{\text{ground}}$ and $N_{\text{excited}}$ by

[TABLE]

When this is substituted into (2.16), the model for Yb-doped fiber amplifiers is complete.

3.3. Basic simulations

We report the results obtained from simulation of the CMT model for two 10 m long fibers, one doped with Yb and the other with Tm. The fiber parameters are collected from data sheets of commercially available exemplars of these fibers (specifically Nufern™ fibers – see nufern.com). All parameters used for the simulation of both the fibers are reported in Tables 1 and 2.

We solve the CMT system (2.16) using the classical $4^{th}$ order explicit Runge-Kutta method (in complex arithmetic). The phase terms $\phi_{lm}(z)=e^{\hat{\imath}(\beta_{m}-\beta_{l})z}$ in the ODE system oscillate at a wavelength not smaller than the so-called mode beat length

[TABLE]

An ODE solver applied to solve (2.16) must take sufficient number of steps per mode beat length to capture the effect of these oscillations in the solution. Prevailing theories [12] point to the potential importance of the mode beating term in thermal effects, so we must be careful to treat these oscillations with the needed accuracy if the model is to be extendable to incorporate thermal effects in the future. In all our simulations, we used 50 ODE steps per mode beat length.

Before running the ODE solver, we precompute the propagation constants $\beta_{j}$, the mode beat length, and of course, the modes. For step-index fibers, we can compute the modes $\varphi_{l}$ exactly in closed form (see [2, 15]) as quickly described next. One first computes the propagation constants by solving the characteristic equation of the fiber as follows. Let $\mathcal{J}_{i}$ and $\mathcal{K}_{i}$ denote, respectively, the standard Bessel function and the modified Bessel function of second kind of order $i$. Then we solve for $X$ satisfying the so-called “characteristic equation” of the fiber, namely setting the fiber’s “numerical aperture” NA$=\sqrt{n_{\text{core}}^{2}-n_{\text{clad}}^{2}},$ we solve $X\mathcal{J}_{i-1}(X)\mathcal{K}_{i}(\sqrt{\text{NA}^{2}-X^{2}})+\sqrt{\text{NA}^{2}-X^{2}}\,\mathcal{J}_{i}(X)\mathcal{K}_{i}(\sqrt{\text{NA}^{2}-X^{2}})=0$ by a bisection-based root-finding method. This equation arises from the matching conditions at the core-cladding interface. For each $i$, enumerating the roots of the characteristic equation as $X_{ij}$, $j=0,1,\ldots$, the propagation constants are given by

[TABLE]

Set $\mathcal{R}_{ij}=X_{ij}/r_{\text{core}}$ and $\mathcal{G}_{ij}=\sqrt{\beta_{ij}^{2}-n_{\text{clad}}k_{s}^{2}}$. The exact LP modes take the following form in polar coordinates:

[TABLE]

The mode $\varphi_{ij}$ is usually called the “LP$ij$” mode.

For the particular case of the Tm parameters in Table 2, we find that the fiber only has the LP01 and LP11 modes, while for the Yb fiber with the parameters set in Table 1, we found four modes LP01, LP11, LP21 and LP02. In our simulation the fiber geometry was meshed using finite elements (with curved elements at the cladding boundary and at the core-cladding interface) and the relevant LP modes were interpolated into the degree $p$ Lagrange finite element space based on the mesh. Integration involving finite element functions is broken into a sum over integrals over all mesh elements and a sufficiently high quadrature rule is used to approximated an element integral. This is how we approximate all required integrals, such as in the computation of the coupling coefficient (2.10), as well as in power computations. Note that each step of the multi-stage ODE solver requires many such integrations.

To quantitatively describe the light amplification results of the simulation, we compute the signal and pump power, after the approximate $Y(z)=[I_{p}(z),{A(z)}^{t}]^{t}$ has been computed, as follows:

[TABLE]

The initial condition $Y(0)$ is set so that the entire signal power is fed into the LP01 mode at the inlet $z=0$. Initial pump power $P_{p}^{0}$ was set 1000 W for the Yb case and 1100 W for the Tm case. Figure 3 shows the distribution of the computed $P_{s}$ and $P_{p}$ (marked “signal” and “pump” there) for the Tm and Yb-doped fibers. The energy transfer from the pump light to the signal light is clearly evident. We used $p=5$ Lagrange elements for these plots. The use of 50 steps per mode beat length implies that the Yb case required 421014 RK4 steps, while the Tm case required 302340 steps of the ODE solver to cover the 10 m fiber.

Each of these hundreds of thousands of steps required (multiple) integrations over the fiber cross section (to compute integrals such as the one in (2.10)). As mentioned above, these integrations were performed using finite element quadratures. In unreported experiments, we have attempted to reduce the cost of these integrations by hyper-reduction techniques common in reduced-order models [16]. One such technique is to use reduced-order quadratures to approximate the cross-section integrals instead of using finite elements to perform the integration precisely. Our pilot studies into this used Gaussian quadrature rules on a disc (core) and an annulus (cladding) of order as high as 20. In cases where this resulted in substantial reductions in computational cost, we unfortunately also observed unacceptably large deviations from the results presented above. Further studies are needed to conclude if other hyper-reduced quadratures, specifically taking the modes into account, might prove more useful. In the next section, we describe a completely different line of inquiry that has yielded considerable acceleration in our simulations.

4. The equivalent short fiber concept

In this section, we present the concept of a nearly equivalent short fiber, which is an artificially short fiber with unphysical parameters that can mimic a longer physical fiber in some respects. Being shorter, the equivalent fiber can be solved using fewer steps of an ODE solver, thus providing significant reductions in computational cost.

To explain the rationale behind the equivalent short fiber approach, first consider applying an ODE solver to solve the CMT model (2.16). As mentioned in the previous section, very large number of ODE steps were needed to solve the CMT system (2.16) on a 10 m long fiber. Therefore, it would be extremely useful to reduce the fiber length (and hence the number of ODE steps) while still preserving the relevant physical processes in the fiber amplifier. We shall now show that this is possible to some extent using the computational scale model of an equivalent short fiber described below.

To begin with, one might consider shortening the $z$-domain in (2.9) using a dimensional analysis. Note that the left hand side of (2.9) has dimension $\mathrm{V/m}$ (volts per meter), and $K_{lm}$ has units of $\mathrm{m}^{-1}$. Therefore, by non-dimensionalization, one is led to believe that a shorter fiber of length $\tilde{L}\ll L$ might, in some ways, behave similarly to the original fiber of length $L$, provided its coupling coefficient is magnified by $L/\tilde{L}$. However, not all nonlinear systems admit scale models that are perfect replicas of the original. Below we shall identify what properties of such a shorter fiber can be expected to be close to the original.

We introduce the variable change

[TABLE]

A fiber of length $L$, under the variable change $\tilde{z}=\zeta^{-1}(z)=z\tilde{L}/L$ becomes one of length $\tilde{L}$. Under this variable change, (2.9) and (2.14) become

[TABLE]

for all $0<\tilde{z}<\tilde{L}$. In other words, defining $\hat{A}_{l}=A_{l}\circ\zeta$ and $\hat{I}_{p}=I_{p}\circ\zeta$, the above system may be rewritten as the following system on the shorter domain $0<\tilde{z}<\tilde{L}$ for $\hat{Y}=[\hat{A}_{l},\hat{I_{p}}]^{t}$,

[TABLE]

Supplemented with the same initial data at at $z=\tilde{z}=0$, (4.3) is exactly equivalent to (2.16), i.e.,

[TABLE]

In other words, the solution of (4.3), being the pull back of the original solution $Y$ to the shorter domain, is a perfect replica of the original solution $Y$.

Unfortunately,  (4.3) on $0<\tilde{z}<\tilde{L}$ offers no computational advantages over the original system (2.16) on $0<z<L$. This is because the mode beat length of (4.3) has been reduced by a factor of $\tilde{L}/L$ due to the variable change. So in order to solve the ODE system (4.3), keeping the same number of steps per mode beat length, the total number of steps needed to solve the system has not been reduced. This leads us to consider another mode coupling system with the same mode beat length as the original system (2.16).

Let $\tilde{Y}(\tilde{z})=[\tilde{I}_{p}(\tilde{z}),\tilde{A}_{1}(\tilde{z}),\cdots,\tilde{A}_{M}(\tilde{z})]^{t}$ solve

$\displaystyle\frac{d\tilde{Y}}{d\tilde{z}}$ $\displaystyle=\begin{bmatrix}\langle(L/\tilde{L})g_{p}(\tilde{Y})\rangle&0\\ 0&p(\tilde{z})\cdot(L/\tilde{L})K(\tilde{Y})\end{bmatrix}\tilde{Y},$ $\displaystyle 0<\tilde{z}<\tilde{L},$

(4.5a)

$\displaystyle\tilde{Y}(0)$ $\displaystyle=[I_{p}^{0},A^{0}]^{t}$ $\displaystyle\tilde{z}=0.$

(4.5b)

Clearly, (4.5) is not the same as (4.3) due to the differences in the phase factors. Therefore, unlike the solution $\hat{Y}$ of (4.3), the solution $\tilde{Y}$ of (4.5) is not a perfect replica of the original solution $Y$. Nonetheless, we shall now proceed to argue that (4.5) is a practically useful scale model of (2.16) as it approximately preserves the power distribution from the original. Power, unlike the amplitude $A$, is the quantity that can be, and actually is, experimentally measured.

Let $P_{l}$ and $\tilde{P}_{l}$ be respectively the powers contained in the $l^{th}$ mode for the physical and equivalent fiber, defined by

[TABLE]

One may express these in terms of

[TABLE]

as $P_{l}(z)=|a_{l}|^{2}\Phi_{l}$, where $a_{l}(z)=A_{l}(z)e^{\hat{\imath}\beta_{l}z}$.

To obtain an equation for $P_{l}(z)$, we may start from the second equation of the block system (2.16), or equivalently from (2.9), which can be rewritten as

[TABLE]

Then using $da_{l}/dz=e^{\hat{\imath}\beta_{l}z}\partial_{z}A_{l}+\hat{\imath}\beta_{l}a_{l},$ we have

[TABLE]

Using also the complex conjugate of this equation, we have

[TABLE]

i.e.,

[TABLE]

for all $l=1,\ldots,M$, or equivalently,

[TABLE]

where

[TABLE]

for $l=1,\ldots,M$.

To the system (4.6), let us also add the pump power using the index $l=0$, i.e., let $P_{0}(z)\equiv P_{p}(z)$ as defined in (3.11). Then integrating (2.14), we obtain $dP_{0}/dz=\langle g_{p}\rangle P_{0}.$ All together, we have thus obtained an equation for $P_{l}$ for all $l=0,\ldots,M$,

[TABLE]

where $P=[P_{0},P_{1},\ldots,P_{M}]^{t}$ and $\text{diag}[\cdot]$ denotes the diagonal part of a matrix.

To understand the motivation for the remaining arguments, we now highlight an observation concerning (4.8). A scale model providing a perfect replica of the original power distribution is easy to obtain if the system (4.8) were an autonomous system: indeed, if there exists a function $F$ of $P$ alone such that $dP/dz=F(P)$, then by merely scaling $F$ by $L/\tilde{L}$, we obtain an equivalent system that provides perfect replicas of the original power distribution on the shorter fiber of length $\tilde{L}$. However (4.8) is not autonomous, in general. Yet, for practical fibers, our numerical experience suggests that (4.8) behaves almost like an autonomous system. Therefore our strategy now is to view (4.8) as a perturbation of an autonomous system.

Of particular interest is the fact that if the fiber amplifier was robustly single-mode ($M=1$ for the laser signal), then the governing system (4.8) would be autonomous. This can be achieved by not using a LMA amplifier, but one of a smaller fiber core size and/or a lower numerical aperture (NA) such that the fiber core can only support only one guided core mode, the fundamental mode (indexed by $m=1$), at the signal wavelength. However, even with a LMA fiber, if one were to account for fiber bending effects, which cause the higher-order core modes (indexed by $1<m\leq M$) to leak into the cladding region more so than for the fundamental mode, then the fiber would operate nearly as a single-mode fiber. Actual fiber amplifiers are almost always wrapped on a spool rather than stretched out straight, thus ensuring this fiber bending effect. This provides us with greater confidence of autonomous system-like behavior, even in real-world implementations of fiber laser amplifier systems.

Recall from (2.10) that $K_{lm}$ is defined using $g_{s}(I_{s},I_{p})$, where $I_{s}$ takes the form in (2.13). We define the following perturbation of $I_{s}$,

[TABLE]

It seems difficult to characterize when $I_{s}-{\mathcal{I}_{s}}$ is small a priori (as it depends, e.g., on the localization and orthogonality of the specific fiber modes) but after a CMT calculation, we may check if this difference is small a posteriori. Deferring for the moment the matter of the size of $I_{s}-{\mathcal{I}_{s}}$, let us proceed to define $\gamma_{\ell}(P)=g_{\ell}({\mathcal{I}_{s}}(P),I_{p})=g_{\ell}({\mathcal{I}_{s}}(P),P_{0}/|\varOmega_{z}|),$ for $\ell\in\{s,p\}.$ They represent the gain functions obtained by replacing $I_{s}$ by ${\mathcal{I}_{s}}$. The new gain functions in turn prompt the definition of a new mode coupling coefficient: instead of (2.10), we now consider

[TABLE]

for all $l,m=1,\ldots,M$. Additionally let

[TABLE]

and $\kappa_{0l}=\kappa_{l0}=0,$ for all $l=1,\ldots,M$. We may now view these $\kappa_{lm}$ as entries of an $(M+1)\times(M+1)$ matrix, using which (4.8) can be expressed as

[TABLE]

where $\eta\in\mathbb{R}^{M+1}$ is defined by

[TABLE]

We view $\eta$ as a function of $z$, i.e., $\eta:[0,L]\to\mathbb{R}^{M+1}$. The $z$-dependence is clear once we express the $z$-dependence of the solution $Y\equiv Y(z)$ and power $P\equiv P(z).$ Equation (4.9) shows that power is governed by a perturbation of an autonomous system whenever $\eta$ is small enough to be viewed as a perturbation.

Returning to consider (4.5), we define analogous quantities for the short fiber, namely

[TABLE]

for $l=1,\ldots,M$. Then we may repeat the above arguments starting from (4.5) to obtain the following analogue of (4.9).

[TABLE]

where $\tilde{\eta}:[0,\tilde{L}]\to\mathbb{R}^{M+1}$ is now given by

[TABLE]

Note that $\rho(\tilde{Y})$ is defined by (4.7) after replacing not only $Y$ by $\tilde{Y},$ but also $a_{l}$ (which depends on $Y$) by $\tilde{a}_{l}$ (which depends on $\tilde{Y}$).

To conclude this analysis, it now suffices to compare (4.10) and (4.9). Applying the change of variable $\zeta$ to (4.9), we get

[TABLE]

Comparing (4.10) and (4.11) we see that when $\eta$ and $\tilde{\eta}$ are negligibly small compared to the other terms, $P_{l}\circ\zeta$ and $P_{l}$ solve approximately the same equation, and consequently

[TABLE]

We summarize this discussion as follows.

The system (4.5) is an equivalent short fiber model of (2.9) in the sense that the power $P_{l}$ contained in the $l^{\text{th}}$ mode is approximately preserved from the original fiber model (2.9) through a change of variable, under the above assumptions.

5. Computational verification of equivalent fiber concept

In this section, we perform extensive numerical experiments to verify the pratical utility of the equivalent fiber concept introduced in Section 4. We shall compare the relative differences in the powers obtained from the original fiber and its equivalent short fiber for various settings to gauge the practical effectiveness of the approximation (4.12). In Subsections 5.1 and 5.2, we show a way to understand the equivalent short fiber as a fiber with artificial parameters (with values not physically realizable) for the Tm and Yb cases, respectively.

5.1. Realizing the equivalent short fiber for the Tm-doped case

The equations of the equivalent short fiber, namely (4.5), can be realized for a dopant medium if we can find a set of “artificial” parameters that would scale the original $g_{p}$ and the original $K$ by $L/\tilde{L}$. In view of (2.10), this effect is achieved by scaling the original $g_{\ell}$ by $L/\tilde{L}$ for $\ell\in\{s,p\}$. Now consider the expressions for $g_{\ell}$ for Tm-doped fiber, given in (3.4) and (3.5). Clearly, in view of these expressions, $g_{\ell}$ will scaled by $L/\tilde{L}$ if all the ion populations $N_{i}$ are so scaled.

This observation, in turn, leads us to consider the expressions for $N_{i}$ we derived in (3.3). Let

[TABLE]

The value of the expression for $N_{0}$ in (3.3a) will be scaled by $L/\tilde{L}$ if we replace $\kappa_{R}$ by $\tilde{\kappa}_{R}$ and $N_{\text{total}}$ by $\tilde{N}_{\text{total}}$, i.e., (3.3a) implies

[TABLE]

Let $\tilde{N}_{0}=LN_{0}/\tilde{N_{0}},$ the left hand side above. Proceeding to analyze the expressions in (3.3b), we find that the same change in $\kappa_{R}$ and $N_{\text{total}}$, and the consequent change in $N_{0}$ to $\tilde{N}_{0}$ per (5.1), also scales all other $N_{i}$ by $L/\tilde{L},$ i.e.,

[TABLE]

Therefore, all the ion populations $N_{i}$ are scaled by $L/\tilde{L}$, and so are $g_{s}$ and $g_{p}$. We have thus arrived at our main observation of this subsection:

A short fiber of length $\tilde{L}$ is equivalent to a Tm-doped fiber of length $L$ if the fiber’s original parameters $N_{\text{total}}$ and $\kappa_{R}$ are changed to $\tilde{N}_{\text{total}}=LN_{\text{total}}/\tilde{L}$ and $\tilde{\kappa}_{R}=\tilde{L}\kappa_{R}/L$, respectively, i.e., this change realizes (4.5).

To see how this idea works in practice, we consider two scenarios, both with an equivalent short fiber of $\tilde{L}=0.1$ m representing the 10 m long Tm fiber we simulated in Figure 3. (All parameters are as in Table 2 except for $N_{\text{total}}$ and $\kappa_{R}$, which were modified for the equivalent fiber as stated above.) In the first scenario, 100% of the input signal power is carried in the LP01 mode at the inlet (the same setting as in the computation reported in Figure 3). In the left panel of Figure 4, we find that the plots of the computed powers for the equivalent short fiber and the real fiber are virtually identical. Even though the difference between them appear to be zero visually, we have quantified this difference in the bottom left plot of Figure 4: since the domains of the two power functions to be compared are different, we pull back the original powers to the shorter domain and plot $P_{l}\circ\zeta-\tilde{P}_{l}$ (for the two modes, LP01 and LP11) on the shorter domain. Clearly, from the scale of the plot, the absolute values of these differences are found to be of the order of $10^{-9}$, so indeed the differences between the two sets of power curves are negligible. The practical value of the equivalent short fiber calculation lies in the fact it gave essentially the same power curves about 100 times faster than the real-length fiber calculation of Figure 3.

In the second scenario, the total input power of 30 W is distributed equally between the LP01 and LP11 modes. From the top right panel of Figure 4, we find that LP01 mode amplifies more than the LP11 mode. Moreover, as in the left panel, the results from the real and equivalent short fiber are visually indistinguishable. However, a more careful examination of the difference $P_{l}\circ\zeta-\tilde{P}_{l}$ in the bottom right plot shows that maximal absolute power differences are about 0.3 near the inlet of the fiber. Although this is many fold larger than the first scenario, the relative power error of $3\times 10^{-4}$ is still quite small enough to make the equivalent short fiber a useful practical tool. Note that the difference $P_{l}\circ\zeta-\tilde{P}_{l}$ is now highly oscillatory, due to the interactions between the two modes.

5.2. Realizing the equivalent short fiber for the Yb-doped case

The equivalent short fiber in the Yb-doped case is more easily realizable than the Tm-case as the Yb population dynamics is simpler. The following conclusion can be arrived at easily proceeding similarly as in Subsection 5.1.

A short fiber of length $\tilde{L}$ is equivalent to a Yb-doped fiber of length $L$ if the fiber’s original parameter $N_{\text{total}}$ is changed to $\tilde{N}_{\text{total}}=LN_{\text{total}}/\tilde{L}$, i.e., this change realizes (4.5).

Figure 5 gives some indication of the practical performance of this equivalent short fiber. As in the experiments for the Tm-fiber reported in Figure 4, here we consider two scenarios, the first where all input signal power is given to the LP01 mode, and the second where the input power is distributed to the four LP modes equally (25% each). The left panel in Figure 5 shows the former, while the right panel shows the latter. The equivalent fiber is less faithful in the latter case, but the scale of the errors observed in the bottom plots in both cases are well within the acceptable error ranges in engineering practice. (Laboratory power measurement uncertainties tend to be about $\pm 5\%$.)

5.3. Increase of error with respect to some parameters

We want to understand how relative power differences between the equivalent and real fiber vary with respect to two important input parameters $P_{p}^{0}$ and the short fiber length $\tilde{L}$. We consider both the Tm and Yb fibers, holding the original fiber length $L$ fixed to 10 m.

The solutions of the original and equivalent fiber models vary as initial conditions are changed. Therefore to compare one with the other in the worst case scenario, we take the maximum of the power error measures over the set

[TABLE]

i.e., the set $\mathcal{A}$ is the set of all input distributions yielding the same initial signal power $P_{s}^{0}$, which is set for Tm and Yb fiber per Tables 2 and 1, respectively. The initial pump power $P_{p}^{0}$ is varied in the range 1000–5000 W (thus providing a corresponding range of initial values for the $I_{p}$-component in the model). We solve the full CMT model and the equivalent short fiber model, not only for this range of $P_{p}^{0}$, but also for decreasing values of the short fiber length $\tilde{L}$. The following quantity is then computed across all such solutions:

[TABLE]

Thus $\varepsilon$ represents the maximal possible power deviations between the equivalent and original models over all input signal distributions and over all mode components, as a function of initial pump power $P_{p}^{0}$ and the fictitious length $\tilde{L}$. Values of $\varepsilon$ will thus inform us of the ranges of $P_{p}^{0}$ and $\tilde{L}$ where the equivalent short fiber is more useful.

To practically compute $\varepsilon$, we replace the maximum over the infinite set $\mathcal{A}$ by a computable maximum over a finite set obtained by assigning each mode component all possible values from 0 to 100% in 10% increments (while constraining the total signal power to $P_{s}^{0}$). In the case of the 2-mode thulium fiber, this resulted in 11 input power distributions, while for the ytterbium-doped fiber having 4 modes, 286 distributions were required. The maximum over $z$ in (5.2) is replaced by the maximum over the points where ODE solver traversed. We used polynomial degree $p=5$ for the finite element approximation of modes and the 7-stage Dormand-Prince Runge Kutta method for solving the ODE system. Collecting data from hundreds of simulations, we then plot $\varepsilon$ in a two-dimensional grid of $P_{p}^{0}$ and $\tilde{L}$ values.

The resulting contour plots of the function $\varepsilon$ are given in Figure 6 for Yb and Tm fibers, for a range of $P_{p}^{0}$ and $\tilde{L}$ values. We find that relative error $\varepsilon$ varies mildly with respect to $P_{p}^{0}$ for any fixed $\tilde{L}$, indicating that the absolute error in the powers increases more or less linearly as $P_{p}^{0}$ is increased. Looking vertically at the plots of Figure 6, we find that holding $P_{p}^{0}$ fixed, there are significant variations in $\varepsilon$ with respect to $\tilde{L}$. The errors definitively increase as $\tilde{L}$ decrease. Figure 6 clearly indicates that excessively short equivalent fiber lengths are not recommendable.