# Steady-periodic method for modeling mode instability in fiber amplifiers

Open Access

## Abstract

We present a detailed description of the methods used in our model of mode instability in high-power, rare earth-doped, large-mode-area fiber amplifiers. Our model assumes steady-periodic behavior, so it is appropriate to operation after turn on transients have dissipated. It can be adapted to transient cases as well. We describe our algorithm, which includes propagation of the signal field by fast-Fourier transforms, steady-state solutions of the laser gain equations, and two methods of solving the time-dependent heat equation: alternating-direction-implicit integration, and the Green’s function method for steady-periodic heating.

© 2013 Optical Society of America

## 1. Introduction

In earlier papers we described a physical processes that can lead to modal instability [1] in high power fiber amplifiers. We also presented simulation results that illustrated how periodic modulation of the pump or signal seed can dramatically reduce the instability threshold [2]. Although we described our model and methods in those papers we did not present the mathematical details. Here we do so. The individual parts of the model are all known, but the way they are combined to form a steady-periodic model of modal instability is original.

First, we offer a brief reprise of the physics included in our fiber amplifier model. Mode instability is the degradation of output beam quality above a sharp power threshold [3]. Below threshold the light emerges in the fundamental mode, above threshold some is in higher order modes, usually LP11. Assuming most of the signal seed light is injected into the fundamental mode LP01 with a small amount populating a higher order mode, usually mode LP11, these two modes interfere along the length of the fiber. Because they have different propagation constants their interference creates a signal irradiance pattern that oscillates along the length of the fiber. Pump light is preferentially absorbed in regions of higher signal irradiance, and because a fraction of the absorbed pump is converted to heat, this creates a heating pattern that resembles the irradiance pattern. The heat pattern is converted to a similar temperature pattern, and the temperature pattern creates a refractive index pattern via the thermo optic effect. If the interference pattern is stationary, there is little or no phase shift between the thermally-induced index pattern and the irradiance pattern, so there is almost zero net transfer of power between modes. However, if the light in the higher order mode is slightly detuned in frequency from that in LP01, the irradiance pattern moves along the fiber - down stream for a red detuning and up stream for a blue detuning. The temperature pattern also moves, but it lags the irradiance pattern, and this lag produces the phase shift necessary for power transfer between modes. Red detuning of LP11 leads to power transfer from LP01 to LP11. The detuning that maximizes the mode coupling is set largely by the thermal diffusion time across the fiber core. This is approximately 1 ms, implying an optimum frequency detuning of approximately 1 kHz.

When light in LP01 is transferred to LP11 by deflection from the moving grating it experiences a frequency shift equal to the frequency offset, due to a Doppler effect [1], so the transferred light adds coherently to the light already in LP11. This mechanism can cause a substantial transfer of the power from LP01 to LP11 if the pump power exceeds a sharp mode instability threshold. This gain process can be categorized as near-forward stimulated thermal Rayleigh scattering.

Our modeling approach is to develop as general a model as feasible, and our model includes all the physical effects just described. We use diffractive beam propagation of a steady-periodic, time-dependent signal field. This field incorporates all fiber modes (including lossy modes) and all offset frequencies that are harmonics of the primary frequency offset. The time-dependent thermal profile is computed using either an alternating-direction-implicit (ADI) integration method or a steady-periodic Green’s function method. Mode coupling occurs through inclusion of the thermally-induced refractive index change in the index profile that is used to compute beam propagation. No analytic or semi-analytic expressions of mode coupling are needed, and coupling between all modes is included. We use this approach because it permits the most accurate modeling, and because it makes adding various physical effects relatively straightforward. The cost of such a general numerical model is long run-times, but using the methods described here one can model several meters of fiber per hour on a consumer grade desktop computer.

Alternative mode coupling models based on similar physical effects have been published by other authors [46] but we will not review them here.

## 2. The scalar, paraxial beam propagation equation

Equation (1) is the scalar, paraxial beam propagation equation for an isotropic medium. It is derived in numerous papers on nonlinear fiber optics, for example [7,8]. The variable Es is a complex envelope function that represents all spatial and temporal modulation of a monochromatic, plane carrier wave. The direction of propagation is z. This equation is appropriate for small index contrasts typical of step index, large mode area fiber amplifiers.

$∂Es(x,y,z,t)∂z=i2kc∇⊥2Es(x,y,z,t)+i[k2(x,y,z,t)−kc2]2kcEs(x,y,z,t)$
where kc = ωcnclad/c is the wave number of the carrier wave in the cladding, and k(x, y, z, t) = ωcn(x, y, z, t)/c. Here n(x, y, z, t) is the guiding index including the thermo-optic contribution,
$n(x,y,z,t)=ncore(x,y,z,t)+dndTT(x,y,z,t),$
where dn/dT is the thermo optic coefficient given in Table 1. The first term on the right hand side of Eq. (1) describes diffraction in a homogeneous medium with a refractive index equal to the cladding index nc. Here $∇⊥2$ is the Laplacian operator in the transverse dimensions x and y. The second term adds a correction to account for the core guidance, including the usual core refractive index step plus the thermally induced index change. In general k(x, y, z, t) is imaginary to account for gain or loss, but we will use k(x, y, z, t) to represent only the real part, and write the much smaller imaginary part as a separate term with the coefficient g(x, y, z, t),
$∂E(x,y,z,t)∂z=i2kc∇⊥2E(x,y,z,t)+i[k2(x,y,z,t)−kc2]2kcE(x,y,z,t)+g(x,y,z,t)E(x,y,z,t).$

Table 1. Model input parameters

## 3. Integrating the propagation equation

Equation (3) can be integrated by various methods - for example, using the split step method described below, or using finite-differences. Whichever method is used, the central issue is how to compute the thermally induced part of k(x, y, z, t) for use in the z integration.

One approach is to integrate the field along the full length of the fiber at t = 0 based on an initial temperature profile T(x, y, z, 0). The temperature profile is then updated based on the heat deposited during the time increment Δt, plus the previous temperature profile and the thermal boundary conditions. The new temperature profile is used to integrate the field at time Δt along the full length of the fiber. This iteration of full z integration of Es alternating with stepping T by Δt is repeated for successive time steps to find the temperature and field histories, both over the domain (x, y, z, t).

Alternatively, a time sequence of fields at z = 0 can be used to compute T(x, y, t) at the input and that time varying temperature can be used to propagate the time varying field by Δz. This is repeated for successive z steps. One limitation is that the temperature is solved in two dimensions rather than three so longitudinal heat flow is not accounted for.

Both of the methods just described can treat transient or steady-periodic cases. However, our goal is to compute behavior only under steady-periodic conditions. We are not interested in starting from an initial transient and following the amplifier evolution to a steady state since that may require integration over a time longer than the thermal diffusion time from the core to the outer diameter of the fiber. The required settling time can be greater than 100 ms, perhaps much greater, depending on the fiber design.

Because transient effects are of secondary interest, and also in order to achieve shorter run times, we base our model on the assumption that all transients have decayed and the only time dependence is periodic modulation of powers, mode content, temperature, etc. With this steady-periodic assumption we need integrate only over a single modulation period of approximately 1 ms. We first perform the time integration on the temperature equation to find T(x, y, t) for a single period at one z position, then we step the signal and pump fields by Δz.

The steady-periodic condition appears to correspond well with reported behavior of amplifiers operating near and slightly above their mode instability thresholds. Far above threshold the mode coupling may be more complicated [3]. The relevant time scale for the steady-periodic condition is the time for heat diffusion over the core rather than over the entire fiber. The temperature outside the core is determined by the time averaged heating in the core, but it cannot respond rapidly enough to follow the periodic heating variations within the core.

#### 3.2. Transverse heat flow approximation

Our method does not allow longitudinal flow of heat. We justify this by noting the large difference in length scale for the transverse and longitudinal temperature variations. The important longitudinal length scale is the modal beat length. This varies with fiber design but it typically falls in the range 3–30 mm in large mode fibers. The transverse length scale is the core diameter which typically lies in the range 20–80 μm. The ratio of longitudinal to transverse length is of order 100:1.

#### 3.3. Narrow line width approximation

Our model relies on a frequency offset between the light in different transverse modes. The frequency offset is usually in the range 300–3000 Hz. The assumption of periodic behavior implies the set of all allowed frequencies must be separated by multiples of the inverse modulation period. The light in any mode must consist of only these frequencies. The width of each frequency component is assumed to be much smaller than the frequency offset. This narrow line width is obviously not entirely realistic. Nevertheless, as we will discuss below in Sec. 17, this method can still be a highly accurate way of treating the problem of mode coupling in relatively narrow band amplifiers.

#### 3.4. Split-step integration in t and z

Because it is only necessary to treat one oscillation period under the steady-periodic assumption, we perform the integration in z on a set of time samples of the field that are evenly spaced over one period. At the fiber input we specify Es(x, y, 0, t) and the pump power. Using this information we compute the temperature profiles for each sampling time over the full period. We describe two methods for the integration of the heat equation below. The temperature profile for each time sample is then used to propagate the corresponding time sample of the signal field by Δz. This is repeated until the end of fiber is reached. This method is used because it makes efficient use of limited computer resources, and because it can be made to run fast.

We use a split-step, fast-Fourier-transform, beam propagation method (FFT BPM) to perform the z integrations. Split-step methods have been employed in numerical simulations of CW optical diffraction, including guiding in optical fibers, at least since the 1970s [7, 8]. One period of the field is discretized into Nt time steps. We use the split-step method by applying it individually to each of the field time slices. Related FFT methods for dispersive propagation are also widely used to model propagation of short light pulses propagating in a single fiber mode [9]. Although we are also propagating pulses consisting of one beat time, it is best to think of the time slices as widely separated samples that are unaffected by dispersion. We will say more about this in Sec. 17.

In the split-step method, Eq. (3) is used to advance each field time sample by Δz in three steps. In the first step, linear diffractive propagation is applied to advance the field by Δz/2 in a homogeneous medium with a refractive index equal to the cladding index. In this step, the beam propagation equation is used, keeping only the diffractive term on the right-hand-side

$∂Es(x,y,z,t)∂z=i2kc∇⊥2Es(x,y,z,t).$
In the second step, the field is advanced by Δz without diffraction, keeping only the phases induced by the guiding and thermally induced index, along with the laser gain and linear loss contained in g(x, y, z, t). The propagation equation for this step becomes
$∂Es(x,y,z,t)∂z=i[k2(x,y,z,t)−kc2]2kcEs(x,y,z,t)+g(x,y,z,t)Es(x,y,z,t).$
The gain term will be described in more detail in the next section. In the third step, linear diffractive propagation is again applied to advance the field by Δz/2. This method produces errors of order 𝒪(Δz3).

## 4. Algorithm

These methods are incorporated in our algorithm which is diagrammed in Fig. 1. In the following sections, we describe each step of the algorithm in more detail.

Fig. 1 Flow diagram for split-step FFT propagation using alternating direction implicit integration (ADI) or steady-periodic Greens function method for temperature calculation.

## 5. Read problem parameters and fiber parameters

The first step in the algorithm is to read in the problem parameters. Table 1 lists them individually. Some of them are standard silica parameters. For example: thermal conductivity K = 1.38 W/m-K; specific heat C = 703 J/kg-K; density ρ = 2201 kg/m3; thermo-optic coefficient dn/dT = 1.2 × 10−5 K−1. Usually the pump wavelength is λp = 976 nm, and the upper state ion lifetime is τ = 901 μs.

We define our problem at each z-location on an (x, y, t) grid. The grid is equally spaced in (x, y, t), with number of points typically (Nx = 64, Ny = 64, Nt = 64). The spatial grid spans a domain of size Lx × Ly with the core at its center, where Lx and Ly are typically ≈3 × dcore. We choose as the thermal boundary condition a fixed temperature on the domain boundary. We assume there is a fixed frequency offset between modes equal to Δω. The period ϒ is defined as one cycle of the beat frequency, ϒ = 2πω. The grid step sizes are then

$Δx=Lx/(Nx−1)$
$Δy=Ly/(Ny−1)$
$Δt=ϒ/Nt.$

The parameters defining the doping profile NYb(x, y), the refractive index profile ncore(x, y) and the linear loss α(x, y) vary with (x, y) position. We usually use super-Gaussian profiles for these, with super-Gaussian coefficient of order 40. The FWHM values for the super-Gaussian diameters are dcore and ddopant. The pump is confined to the pump cladding of diameter dclad. The absorption and emission cross sections for pump and signal are $σpa$, $σpe$, $σsa$, $σse$.

The forward and backward time-averaged input pump powers, $Pfp$ and $Pbp$, may be periodically modulated by specifying a pump modulation function Mp(t). Similarly, each signal mode LPm,n with time averaged power $Pm,ns$ can be modulated by specifying its modulation function $Mm,ns(t)$. The γm,n are mode specific losses (non heating).

The fiber length is L, and the fiber is bent to a radius of Rbend. The integration step size Δz is typically set to a few microns. We store information about the field every Δsample integration steps. Typically the stored information includes the local time sequences for the mode content, total signal power, effective area of the signal, and pump power.

## 6. Calculate modes

The next step of the algorithm is to compute the signal mode profiles. Here, we present a general method capable of handling bent fiber and arbitrary refractive index profiles of low contrast. This procedure is from Marcuse [10]. It is not a vector method so it is not appropriate for air holes or other guiding structure with high index contrast. Vector adaptations that permit high index contrast are available, but we do not discuss them here.

On the rectangular grid (0 < x < Lx) by (0 < y < Ly), the functions

$ϕμν=4LxLysin(πLxμx)sin(πLyνy)$
with μ and ν being integers, form a complete, orthonormal set, normalized so their areal integral is unity. Any transverse mode ψ(x, y, z) can be written
$ψ(x,y,z)=e−iβz∑μ=1∞∑ν=1∞Cμνϕμν(x,y)$
where β is the propagation constant for that mode, Cμν is the mode’s eigenvector (coefficients of the basis set ϕμν from Eq. (9)) and ψ is a solution of the Helmholtz equation
$∂2ψ∂x2+∂2ψ∂y2+∂2ψ∂z2+n2(x,y)k∘2ψ=0$
where k = ωc/c. We set an upper limit on the number of sin-sin terms to include in the expansion of ψ(Mx and My), and substituting this expansion into the Helmholtz equation yields
$∑μ=1Mx∑ν=1MyCμν[n2(x,y)k∘2−β2−π2(μ2Lx2+ν2Ly2)]ϕμν(x,y)=0.$
Next multiply this equation by ϕμν and integrate it over x and y. We define neff = β/k and we divide the equation by $k∘2$ and add and subtract $nclad2$ under the summation sign to obtain
$∑μ=1Mx∑ν=1MyAμ′ν′,μνCμν=(neff2−nclad2)Cμ′ν′.$
The A matrix is
$Aμ′ν′,μν=∫0Lx∫oLydxdy{[n2(x,y)−nclad2]ϕμ′ν′(x,y)ϕμν(x,y)}−π2k∘2(μ2Lx2+ν2Ly2)δμ,μ′δν,ν′.$
Equation (13) is an eigenvalue equation of the form
$AC=n¯2C$
to be solved for eigenvalues n̄ and eigenvectors C. The solution eigenvectors are substituted in Eq. (10) to compute the eigenmodes. Then the effective refractive index neff of mode LPm,n is related to its eigenvalue n̄m,n by
$neff(m,n)=n¯m,n2+nclad2.$
The beat length between LP01 and LP11 can be found using their values for neff.

For simple guiding index shapes such as a step index in unbent fiber, the integrals in Eq. (14) can be integrated analytically, otherwise they can be evaluated numerically. We use numerical integration because of its generality.

The modes computed this way are the low power modes that do not take into account a nonuniform temperature or any other irradiance-dependent index changes. These low power basis modes will be used to construct input fields and to analyze propagating fields into modes.

#### 6.1. Refractive index profile

The method described is general enough to consider arbitrarily-shaped refractive index profiles. For simplicity, we typically use a top hat-like two-dimensional super-Gaussian of high order to construct our refractive index and doping profiles. The super-Gaussian index profile is computed using

$n(x,y)=nclad+(ncore−nclad)×exp(−ln2{[2(x−x0)/dcore]2+[2(y−y0)/dcore]2}S),$
where nclad and ncore are the refractive indices of the cladding and core, (x0, y0) are the coordinates of the core center, dcore is the core diameter, and S is the super-Gaussian coefficient (we typically use S = 40).

We approximate the effects of bending the fiber by adding a refractive index ramp nbend to n(x, y) in the plane of the bend so the index is higher on the outside of the bend. If bending is in the x̂ plane, this added ramp can be written

$nbend(x,y)=n(x,y)(x−x0)/Rbend.$
Bend losses are calculated using the method of Marcuse [11]. In the numerical model bend and other losses from the core are enforced by included an absorbing layer near the grid boundary.

#### 6.2. Normalization

In order to easily relate our calculated modes to power in watts, we introduce a normalization factor Nm,n for each mode constructed from Eq. (10), where Nm,n is the value which satisfies

$cε02∫dxdy n(x,y)[Nm,nψm,n(x,y)]2=1W$
The normalized mode um,n(x, y) is defined by
$um,n(x,y)=Nm,nψm,n(x,y).$

## 7. Set up temperature solver

Next the temperature profile must be solved for each time point in the period ϒ using the periodic heat source Q. In an isotropic homogeneous medium, the heat equation is

$ρC∂T∂t=Q+K(∂2T∂x2+∂2T∂y2)$
where T = T(x, y, t) is the temperature, Q = Q(x, y, t) is the heating source, K is the scalar thermal conductivity, ρ is the density and C is the specific heat capacity. No 2/∂z2 term is included in Eq. (21) because we don’t allow longitudinal heat flow. We consider only the thermal boundary condition of a fixed temperature on the (x, y) grid boundary. Alternative boundary conditions can be dealt with by modifying our procedure, but we do not discuss them here.

We will describe two methods of solving Eq. (21) for T. The first method, the Green’s function method for steady-periodic heating, involves computing the temperature over the (x, y) grid as a sum of independent contributions, one from each heated grid point. The second method, the alternating direction implicit (ADI) integration, involves stepping in time via successive matrix multiplications involving all pixels contributing together.

#### 7.1. Calculate Green’s functions for steady-periodic heating

To use the Green’s function method we must first compute the steady-periodic Green’s function which gives the temperature contributions over the entire (x, y) grid due to the periodically heated point (x′, y′). These functions describe the temperature rise due to a steady-periodic heat source at a single modulation frequency. Formulations of the Green’s functions for different boundary conditions are detailed in [12]. For temperatures clamped at all four sides of the grid, the complex valued Green’s function is

$G(x,y,ω|x′,y′)=∑n=0∞Fn(y,y′)Pn(x,x′,ω)$
$Fn(y,y′)=1WKsin(nπWy)sin(nπWy′)$
$Pn(x,x′,ω)={exp[−σn(2H−|x−x′|)]−exp[−σn(2H−x−x′)]+exp[−σn|x−x′|]−exp[−σn(x+x′)]}/σn(1−exp[−2σnH])$
where
$σn2=(nπW)2+iωρCK,$
ω is the frequency of the heat, (0 ≤ xH), and (0 ≤ yW).

We compute these functions for several harmonics, making sure the time grid is sufficient to resolve them. That is, we find G(x, y, ωm|x′, y′) for ωm = (0, 1, 2, 3,...) × Δω where Δω is the specified frequency offset. The use of this set of Green’s functions to compute T(x, y, t) is described in Sec. 13.1.

An alternative, more general method of solving the heat equation is the alternating direction implicit (ADI) method [13]. This method does not require periodicity. When it is applied to a steady-periodic problem, it requires many iterations to match the steady-periodic requirement, with each iteration operating on the previous iteration’s data. This sequential process means a simple parallel computing scheme is difficult to implement for ADI integration.

The derivation of the equations for this method begins by using finite difference definitions of the partial derivatives in Eq. (21). The numerical integration of the heat equation by time step Δt is split into two halves. We first integrate explicitly in y and implicitly in x for one half time step Δt/2, then integrate explicitly in x and implicitly in y for another half time step Δt/2.

$Tx,yt+Δt/2−Tx,yt=(Tx+Δx,yt+Δt/2−2Tx,yt+Δt/2+Tx−Δx,yt+Δt/2)λx2+(Tx,y+Δyt−2Tx,yt+Tx,y−Δyt)λy2+Δt2ρCQt+Δt/4$
or, rewritten in matrix form,
$Tt+Δt/2Ax=ByTt+Δt2ρCQt+Δt/4$
$Tt+Δt/2=ByTtAx−1+Δt2ρCQt+Δt/4Ax−1,$
where Ax is the matrix for implicit integration in the x-direction, and By is the matrix for explicit integration in the y-direction. These matrices are tridiagonal, and Ax has size (Nx − 2) × (Nx − 2) while By has size (Ny − 2) × (Ny − 2). The matrices are
$Ax=[(1+λx)−λx/20⋯000−λx/2(1+λx)−λx/2⋯000⋮⋮⋮⋱⋮⋮⋮000⋯−λx/2(1+λx)−λx/2000⋯0−λx/2(1+λx)],$
$By=[(1−λy)λy/20⋯000λy/2(1−λy)λy/2⋯000⋮⋮⋮⋱⋮⋮⋮000⋯λy/2(1−λy)λy/2000⋯0λy/2(1−λy)].$
Matrix Ax has diagonal elements (1 + λx) and off-diagonal elements −λx/2, while matrix By has diagonal elements (1 − λy) and off-diagonal elements λy/2, where
$λx=KΔtρCΔx2$
$λy=KΔtρCΔy2.$

In the second half time step, implicit/explicit integration order is reversed

$Tx,yt+Δt−Tx,yt+Δt/2=(Tx+Δx,yt+Δt/2−2Tx,yt+Δt/2+Tx−Δx,yt+Δt/2)λx2+(Tx,y+Δyt+Δt−2Tx,yt+Δt+Tx,y−Δyt+Δt)λy2+Δt2ρCQt+3Δt/4.$
Or, rewritten in matrix form,
$AyTt+Δt=Tt+Δt/2Bx+Δt2ρCQt+3Δt/4Tt+Δt=Ay−1Tt+Δt/2Bx+Δt2ρCAy−1Qt+3Δt/4,$
where Ay and Bx follow a definition similar to Eqs. 28 and 29, with λx swapped with λy and redimensioned as appropriate.

Combining the two half time step integrations, using the ADI to integrate by one full time step Δt becomes

$Tt+Δt=Ay−1ByTtAx−1Bx+Δt2ρCAy−1(Qt+Δt/4Ax−1Bx+Qt+3Δt/4),$
where, because we only have Q defined on integral numbers of time steps, we use linear interpolation to find Qtt/4 and Qt+3Δt/4 from the heat at t and t + Δt. This method is converged to 𝒪(Δt2).

## 8. Construct propagation phase array

Linear diffractive propagation that comprises the first and third steps in the split-step propagation is best evaluated in (kx, ky, z,t) space. Equation (4) can be transformed to this space by inserting the following form of Es(x, y, z, t)

$Es(x,y,z,t)=12π∫−∞∞∫−∞∞Es(kx,ky,z,t)eikxxeikyydkxdky.$
The transformed equation is
$∂Es(kx,ky,z,t)∂z=−ikx22kcEs(kx,ky,z,t)−iky22kcEs(kx,ky,z,t).$
Advancing the field E(kx, ky, z,t) by Δz/2 consists of shifting the phase of each (kx, ky) plane wave component by
$ϕ(kx,ky)=−Δz2(kx2+ky2)2kc.$
The inverse two-dimensional Fourier transform is used to convert Es(kx, ky, z,t) back to Es(x, y, z, t). In practice, fast Fourier transforms (FFTs) are used to efficiently convert fields between (x, y, z, t) and (kx, ky, z, t) spaces.

## 9. Construct input signal field and pump powers

#### 9.1. Signal

We use the set of normalized modes um,n(x, y) defined in Eq. (20) to construct the input signal field

$Es(x,y,z,t)|z=0=∑modesPm,nsum,n(x,y)Mm,ns(t)$
where $Pm,ns$ is the power in the (m, n)th mode and $Mm,ns(t)$ the periodic modulation function for that mode. For example, to model an amplifier with LP01 populated by unshifted light, and LP11 populated by light detuned by Δω, we set $M0,1s(t)=1$, $M1,1s(t)=exp(−iΔωt)$. If instead we wish to amplitude modulate the signal in both modes with a simple sinusoidal envelope, we make each modulation function $Mm,ns(t)=1+0.25acos(Δωt)$ where a is the small peak-to-peak power modulation. We can also use more complicated modulations, as long as they are periodic. For more complicated modulation we must include an adequate number of harmonics in the Green’s function temperature solver.

#### 9.2. Pump

Defining a co-propagating pump input is simple, but since we start the integration at the signal input end, specifying the counter-propagating pump that gives the target pump at the opposite end can be more difficult, especially if amplitude modulations are involved. To keep the two pumps separate we define two quantities $Pfp(z,t)|z=0$ and $Pbp(z,t)|z=0$ for forward- and backward-propagating pumps.

To generate modulated input pump powers, we follow a technique similar to the signal modulation considered above for each pump. The model is capable of treating an arbitrary periodic pump modulation, provided we include an adequate number of harmonics if solving the heat equation using the Green’s function method.

## 10. Propagate signal field

To reiterate the propagation described earlier, the signal is propagated by Fourier transforming each time slice of the field Es(x, y, z, t) to Es(kx, ky, z, t) using 2D-FFTs. The phase of each plane wave component is advanced by the propagation phase ϕ(kx, ky) given in Eq. (37), and the inverse 2D-FFTs of Es(kx, ky, z, t) gives the updated field Es(x, y, z, t). If the propagation step is not the first or last half-step in the fiber, the two consecutive half-step propagations are combined by advancing the phase by 2ϕ(kx, ky).

## 11. Update the field for laser gain

In the non diffractive part of the split-step procedure we update the signal field by adding the guiding and thermally induced phases plus any laser gain and linear loss. To compute the gain we find the upper state population profile in (x, y, t) from the signal irradiance profile, pump power, and doping profile, and use it to compute the signal gain and pump loss. The (steady-state) upper state population is

$nu(x,y,t)=Ppσpa/hνpAp+Isσsa/hνsPp(σpa+σpe)/hνpAp+Is(σsa+σse)/hνs+1/τ$
where Is = Is(x, y, t) is the signal irradiance, νs and νp are the signal and pump frequencies, $σsa$ and $σse$ are the signal absorption and emission cross sections, $σpa$ and $σpe$ are the pump absorption and emission cross sections, τ is the ion upper-state lifetime, Pp is the pump power, and Ap is the area of the pump cladding. We assume that the cross sections and lifetime are independent of temperature. The effective ion lifetime at typical amplifier powers is a few micro seconds so at modulation frequencies of order 1 kHz the steady state solution is appropriate.

The change in signal field is given by

$∂Es(x,y,t)∂z=12[−σsa+(σsa+σse)nu(x,y,t)]NYb(x,y)Es(x,y,t)−12α(x,y)Es(x,y,t),$
where NYb(x, y) is the Yb3+ doping profile, α(x, y) is the linear absorption coefficient. This method is general enough to treat an arbitrarily-shaped doping profile, but we typically consider super-Gaussian doping profiles similar to the one defined in Eq. (17). The linear absorption coefficient α can be non-uniform in (x, y) to accommodate a photodarkening model. We assume that all the power absorbed due to α(x, y) is turned into heat. Referring to Eq. (5), the laser gain and loss term g(x, y, t)Es(x, y, t) is given by the right hand side of Eq. (40).

## 12. Update pump powers

We include both forward- and backward-going pumps, with the total pump power given by

$Pp=Pfp+Pbp.$
Pp is assumed to be uniformly distributed across the pump cladding. The change in the pump power is computed directly from the ion inversion, rather than from the signal increment. This allows us to include linear signal loss and fluorescence loss correctly. The total change in pump power is given by
$dPpdz=PpAp∬[(σpa+σpe)nu(x,y)−σpa]NYb(x,y)dxdy.$
The loss is apportioned between the forward and backward pumps according to
$dPfpdz=dPpdzPfpPp$
$dPbpdz=−dPpdzPbpPp.$

## 13. Compute T(x,y,t) and thermo-optic index

For the computation of T(x, y, t) we use either the Green’s function method or the ADI method. We calculate the heat deposition rate from the absorbed pump and the quantum defect according to

$Q(x,y,t)=NYb(x,y)[νp−νsνp][σpa−(σpa+σpe)nu(x,y,t)]Pp(t)Ap+α(x,y)Is(x,y,t),$
where the upper state population nu(x, y, t) is given by Eq. (39).

#### 13.1. Green’s function method

The heat as a function of time over the full cycle at each (x′, y′) pixel is resolved into its Fourier components ωm = (0, 1, 2, 3...)×Δω by performing a temporal Fourier transform on Q(x′, y′, t)

$q(x′,y′,ωm)=ΔxΔy∑i=0Nt−1Q(x′,y′,ti)exp(−iωmti).$
Here, q(x′, y′, ωm) is a complex coefficient that includes the phase as well as the amplitude of the heat deposition. These q(x′, y′, ωm) values are used to weight the Green’s functions in computing the steady-periodic temperature over the entire transverse grid. We always include frequencies (0, 1)×Δω, and we add higher frequency terms as needed for convergence. The time grid must be fine enough to resolve the highest frequency. Higher frequency terms are usually necessary only above the mode instability threshold.

Using the Green’s functions computed as described in Sec. 7.1, the temperature T(x, y, t) is found by summing the contributions of each q(x′, y′, ωm) according to

$T(x,y,t)=Real[∑m=0max∑x′,y′q(x′,y′,ωm)G(x,y,ωm|x′,y′)exp(iωmt)].$

In Sec. 7.2, we described a method of integrating the heat equation using the ADI method. This method uses Eq. (34) to integrate one time step of Δt. For steady-periodic heating, we enforce the steady-state criterion by integrating for several periods, reusing the heat Q(x, y, t) from one period to the next. We terminate this process when the difference in temperatures between periods has reached an acceptably small residual.

The ADI method is more general than the steady-periodic Green’s function method, because it does not require periodicity. It can be used to model transient heating, for example, if the steady-periodic condition is not enforced. It is unconditionally stable, but its accuracy suffers if the time-step is too large.

## 14. Add thermo-optic and guiding index phases to the field

In the non diffraction portion of the split-step propagation we also advance the phase of the field to account for the guiding index plus the thermal index over Δz using

$ϕ(x,y,t)=Δzk2(x,y,t)−kc22kc.$
The phase can be rewritten using Eq. (2) as
$ϕ(x,y,t)=Δzωcc([n(x,y,t)−nclad]+[n(x,y,t)−nclad]22nclad).$
We write it in this form merely to show that the phase is (ωcΔzΔn/c) plus a small correction from the quadratic term.

## 15. Resolve time-dependent modal power contents

We use the normalized modes defined in Eq. (20) to resolve the content of the signal field Es(x, y, z, t) into fiber modes. We compute the inner product of the field and the normalized mode from

$Fm,n(z,t)=∬dxdyEs(x,y,z,t)um,n(x,y)∬dxdy[um,n(x,y)]2,$
where Fm,n(z, t) is a complex quantity. The mode amplitude Fm,n(z, t) is converted to power in watts using
$Pm,n(z,t)=|Fm,n(z,t)|2.$

#### 15.1. Spectral analysis

To find the frequency spectrum of mode (m, n) we Fourier transform Fm,n(z, t) from time to frequency. The allowed frequencies for the periodic function Fm,n(z, t) are the ωm.

#### 15.2. Compute effective area

We also compute the effective area of the signal field using

$Aeff(z,t)=[∬|Es(x,y,z,t)|2dxdy]2∬|Es(x,y,z,t)|4dxdy.$

## 16. An example computation

In Fig. 2 we show a sample result produced by the model. The surface shows the power in LP11 over one LP01-LP11 beat length at the input end of an amplifier versus t and z. The motion of the surface ridges reflect the change in phase of LP11 relative to LP01 due to the difference in propagation constants for the two modes. In this example the LP11 seed light is Stokes shifted by 600 Hz relative to LP01, but the details of the fiber are unimportant for this illustration since all large mode fiber amplifiers produce qualitatively similar surfaces. A small fraction of the gain is due to laser gain but most is due to thermally induced mode coupling.

Fig. 2 Power in LP11 versus time (T) and distance (Z) at the input end of the amplifier similar to the one specified in [2]. This illustrates the growth of the low power mode over one beat cycle and over one beat length. The pump is unmodulated in this example.

The offset frequency that produces the highest mode coupling gain varies along the length of the fiber due to the changing shape of the heat profile. We integrate a set of frequencies and powers over the full length of the fiber to find the lowest threshold. When operating near or below threshold the gain curve has a well defined frequency of highest gain even though the shape of the gain varies somewhat along the fiber due to changing heat profiles.

## 17. A narrow band width model for broad band light?

Our model assumes that inter modal frequency shifts of order 1 kHz are responsible for thermal mode coupling. How is it possible that a model based on such small shifts can be applicable to light that has much broader line widths? Few seed lasers have sub kHz line widths, and in fact, the seed light is often intentionally phase modulated to widths of several GHz to combat SBS. Our answer is that even broad band light can be frequency shifted as a whole by 1 kHz. Further, for phase modulated light the beat between two waves with identical, but frequency shifted spectra is just as strong as the beat between two frequency shifted monochromatic waves. Additionally, the frequency shift in our model is caused by coupling two modes via a moving thermal grating. The Doppler frequency shift induced by this process shifts the entire spectrum of the scattered light so the interference between modes has full strength. This means that both in the mode coupling process and in the mode beating process responsible for heating, light with a broad spectrum behaves the same as monochromatic light. Also, a phase modulation of several GHz has little impact on the ion population evolution. The steady state expression for the upper state population is still accurate.

This equivalence of narrow and broad band light breaks down if the two interfering fields shift out of time synchronization. This can happen for large band widths because different modes have slightly different group velocities. The difference in modal propagation times is of order 1 ps/m, so an estimate of the line width limit is 100 GHz for a typical fiber amplifier. However, this is substantially larger than typical SBS suppressing line widths.

If the signal light consists of a train of short pulses, and the pulse separation is shorter than a few micro seconds, the ion population responds primarily to the average power. Slow modulations of the average power, in the kHz range, are still effective in driving the thermal grating that leads to mode coupling. If the pulse width is more than a ps or so the pulses in different modes will stay synchronized well enough to produce strong inter modal interference.

## 18. Sources of frequency offset HOM seed light

It is likely that in most amplifiers an amplitude or spectral modulation of the pump light [2], combined with accidental injection of a small fraction of seed light into LP11, produces the initial frequency shifted light in LP11. Alternatively, amplitude modulation of the seed, plus accidental injection of a fraction into LP11 could provide the initial light. If both of these modulations can be sufficiently suppressed, the initial light would be provided by quantum noise. This unavoidable initial light would produce the highest possible instability threshold. Our model does not incorporate a true quantum noise model. Instead we estimate the starting noise power from a classical stochastic electrodynamics noise model [14] often used in estimating thresholds for lasing or for Raman gain. The noise level is set to half a photon per mode. In our case there is a single transverse mode and a gain bandwidth of approximately 1 kHz. The expression for starting power is the same as for Raman noise power,

$P=hνΔν=(6.6×10−34)(2.8×1014)Δν=1.8×10−19Δν$
which gives 10−16 W for a bandwidth of 0.5 kHz and a wavelength of 1060 nm. This is the minimum starting power in the higher order mode. The gain process will pick out from the sea of quantum noise the light that best matches the light in LP01 both in amplitude and phase over the beat cycle. This is the light with highest gain. There will be fluctuations in frequency and power, but the average power within the gain bandwidth will be well approximated by Eq. (53).

## 19. Desktop computer implementation

#### 19.1. Memory requirements

Storing the full, four-dimensional double-precision, complex [64 × 64 × 64 × Lz] array of the signal field would require on the order of 1 – 10 terabytes in a meter long amplifier with step sizes of a few microns. Therefore, we do not store E fields at each step. We only store computed properties such as modal content Fm,n(z, t), total signal power I(z, t), and effective area Aeff(z, t) at positions separated by (Δsample × Δz). This reduces the amount of memory required to at most a few gigabytes.

#### 19.2. Parallel computing

A substantial portion of the model’s run time is spent solving the thermal problem. This makes the Green’s function method attractive because it is generally much faster than the ADI method. One reason the ADI is relatively slow is that it is necessary to integrate through several cycles of the heat to ensure steady-periodic condition is enforced.

Equally important, the Green’s function method is easy to parallelize while the ADI method is not. The ADI method is sequential by nature. It requires updating the entire grid to advance the time by Δt. In contrast, the Green’s function involves a summation of the temperature contribution from each heated pixel, so the calculations for pixels are independent and easy to run in parallel. This parallelization is relatively simple to implement using the shared-memory multiprocessing library OpenMP.

In addition, we can parallelize other steps of the model. Propagation of the signal field array can be performed independently for each time slice of the field. Similarly, the laser gain calculation, modal content decomposition, and effective area calculation can be performed independently for each of the time slices.

#### 19.3. Execution speed considerations

We have both MATLAB and Fortran versions of our model. Both versions use the FFT library FFTW [15]. Our experience is that FFTW is substantially faster than other FFT routines in Fortran.

OpenMP, the shared-memory multiprocessing library and compiler directives, is implemented in the version of Fortran we use, GNU Fortran 4.6.3. Using a desktop computer based on an Intel Core-i7 3770 (Ivy-Bridge) processor with four physical cores, we are able to obtain an excellent speed up. Using four threads instead of one, the time required to run the same fiber setup decreases by a factor of ≈ 3.2.

Another important speedup was obtained by solving the thermal problem and laser gain equations on a coarser z-grid than the FFT propagation problem. This can dramatically reduce the model run time as well. Care in choosing this grid is important because a grid that is too coarse can reduce the computed gain and increase the instability threshold power.

The run times on our computer lie in the range of 0.25–1.5 hour/m depending mostly on the z step size and the number of harmonics included in the Green’s function. Larger cores use larger z steps, and near-threshold runs require only a single harmonic, permitting times near 0.25 hour/m.

## 20. Approximations of model

All numerical models make judicious approximations. A brief list of ours follows:

• Single λ pump with single absorption and emission cross sections
• Pump power is uniformly distributed across the pump cladding
• All signal light is identically polarized
• Steady-periodic heating is required for application of Green’s function
• Fixed period ϒ, which allows only signal frequency offsets 1/ϒ × (0, ±1, ±2,...)
• Thermal boundary condition is fixed temperature on square boundary
• Thermal boundary size approximately three times core diameter
• Heat equation solved in two dimensions (no longitudinal heat flow)
• Thermal properties assumed uniform and isotropic
• Temperature dependence of cross-sections, dn/dT, and τ not included
• No refractive index dependence on nu
• Signal bandwidth < 100 GHz
• Low contrast refractive index profile, e.g. no air-holes as in PCF
• Upper state population nu follows Is instantaneously (steady-state expression)

## 21. Attributes of model

Attributes of our model include:

• Highly numeric - general and simple to add additional physical effects
• Model a variety of refractive index, linear absorption, and doping profiles
• Steady-periodic model produces well defined thresholds
• The ADI method can be used to study transient behavior if desired
• All transverse modes are automatically included
• Thermal lensing automatically included
• Comparatively short run times and minimal memory requirements
• Variety of thermal boundary conditions possible using Green’s functions
• Green’s function method offers large speedup using multiple processors

## Acknowledgments

This work partially supported under funding from the Air Force Research Laboratory Directed Energy Directorate.

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

2. A. V. Smith and J. J. Smith, “Influence of pump and seed modulation on the mode instability thresholds of fiber amplifiers,” Opt. Express 20, 24545–24558 (2012). [CrossRef]   [PubMed]

3. H.-J. Otto, F. Stutzki, F. Jansen, T. Eidam, C. Jauregui, J. Limpert, and A. Tünnermann, “Temporal dynamics of mode instabilities in high-power fiber lasers and amplifiers,” Opt. Express 20, 15710–15722 (2012). [CrossRef]   [PubMed]

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

5. B. Ward, C. Robin, and I. Dajani, “Origin of thermal modal instabilities in large mode area fiber amplifiers,” Opt. Express 20, 11407–11422 (2012). [CrossRef]   [PubMed]

6. C. Jauregui, T. Eidam, H.-J. Otto, F. Stutzki, F. Jansen, J. Limpert, and A. Tünnermann, “Physical origin of mode instabilities in high-power fiber laser systems,” Opt. Express 20, 12912–12925 (2012). [CrossRef]   [PubMed]

7. M. D. Feit and J. A. Fleck, “Computation of mode properties in optical fiber waveguides by a propagating beam method,” Appl. Opt. 19, 1154–1164 (1980). [CrossRef]   [PubMed]

8. M. D. Feit and J. A. Fleck, “Computation of mode eigenfunctions in graded-index optical fibers by the propagating beam method,” Appl. Opt. 19, 2240–2246 (1980). [CrossRef]   [PubMed]

9. G. P. Agrawal, Nonlinear fiber optics, second ed. (Academic Press, 1995).

10. D. Marcuse, Theory of dielectric optical waveguides, 2nd ed. (Academic Press, 1991).

11. D. Marcuse, “Curvature loss formula for optical fibers,” J. Opt. Soc. Am. 66, 216–220 (1976). [CrossRef]

12. K. D. Cole, “Steady-periodic Green’s functions and thermal-measurement applications in rectangular coordinates,” J. Heat Trans. 128, 706–716 (2006); DOI: [CrossRef]  .

13. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd ed. (Cambridge Univ. Press, 1992).

14. P. W. Milonni, The quantum vacuum: an introduction to quantum electronics (Academic Press, 1994).

### References

• View by:

1. A. V. Smith and J. J. Smith, “Mode instability in high power fiber amplifiers,” Opt. Express 19, 10180–10192 (2011).
[Crossref] [PubMed]
2. A. V. Smith and J. J. Smith, “Influence of pump and seed modulation on the mode instability thresholds of fiber amplifiers,” Opt. Express 20, 24545–24558 (2012).
[Crossref] [PubMed]
3. H.-J. Otto, F. Stutzki, F. Jansen, T. Eidam, C. Jauregui, J. Limpert, and A. Tünnermann, “Temporal dynamics of mode instabilities in high-power fiber lasers and amplifiers,” Opt. Express 20, 15710–15722 (2012).
[Crossref] [PubMed]
4. K. R. Hansen, T. T. Alkeskjold, J. Broeng, and J. Laegsgaard, “Thermally induced mode coupling in rare-earth doped fiber amplifiers,” Opt. Lett. 37, 2382–2384 (2012).
[Crossref] [PubMed]
5. B. Ward, C. Robin, and I. Dajani, “Origin of thermal modal instabilities in large mode area fiber amplifiers,” Opt. Express 20, 11407–11422 (2012).
[Crossref] [PubMed]
6. C. Jauregui, T. Eidam, H.-J. Otto, F. Stutzki, F. Jansen, J. Limpert, and A. Tünnermann, “Physical origin of mode instabilities in high-power fiber laser systems,” Opt. Express 20, 12912–12925 (2012).
[Crossref] [PubMed]
7. M. D. Feit and J. A. Fleck, “Computation of mode properties in optical fiber waveguides by a propagating beam method,” Appl. Opt. 19, 1154–1164 (1980).
[Crossref] [PubMed]
8. M. D. Feit and J. A. Fleck, “Computation of mode eigenfunctions in graded-index optical fibers by the propagating beam method,” Appl. Opt. 19, 2240–2246 (1980).
[Crossref] [PubMed]
9. G. P. Agrawal, Nonlinear fiber optics, second ed. (Academic Press, 1995).
10. D. Marcuse, Theory of dielectric optical waveguides, 2nd ed. (Academic Press, 1991).
11. D. Marcuse, “Curvature loss formula for optical fibers,” J. Opt. Soc. Am. 66, 216–220 (1976).
[Crossref]
12. K. D. Cole, “Steady-periodic Green’s functions and thermal-measurement applications in rectangular coordinates,” J. Heat Trans. 128, 706–716 (2006); DOI: .
[Crossref]
13. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd ed. (Cambridge Univ. Press, 1992).
14. P. W. Milonni, The quantum vacuum: an introduction to quantum electronics (Academic Press, 1994).
15. M. Frigo and S. G. Johnson, “FFTW Home Page,” http://www.fftw.org/ .

#### 2006 (1)

K. D. Cole, “Steady-periodic Green’s functions and thermal-measurement applications in rectangular coordinates,” J. Heat Trans. 128, 706–716 (2006); DOI: .
[Crossref]

#### Agrawal, G. P.

G. P. Agrawal, Nonlinear fiber optics, second ed. (Academic Press, 1995).

#### Cole, K. D.

K. D. Cole, “Steady-periodic Green’s functions and thermal-measurement applications in rectangular coordinates,” J. Heat Trans. 128, 706–716 (2006); DOI: .
[Crossref]

#### Flannery, B. P.

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd ed. (Cambridge Univ. Press, 1992).

#### Marcuse, D.

D. Marcuse, Theory of dielectric optical waveguides, 2nd ed. (Academic Press, 1991).

#### Milonni, P. W.

P. W. Milonni, The quantum vacuum: an introduction to quantum electronics (Academic Press, 1994).

#### Press, W. H.

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd ed. (Cambridge Univ. Press, 1992).

#### Teukolsky, S. A.

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd ed. (Cambridge Univ. Press, 1992).

#### Vetterling, W. T.

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd ed. (Cambridge Univ. Press, 1992).

#### J. Heat Trans. (1)

K. D. Cole, “Steady-periodic Green’s functions and thermal-measurement applications in rectangular coordinates,” J. Heat Trans. 128, 706–716 (2006); DOI: .
[Crossref]

#### Other (5)

G. P. Agrawal, Nonlinear fiber optics, second ed. (Academic Press, 1995).

D. Marcuse, Theory of dielectric optical waveguides, 2nd ed. (Academic Press, 1991).

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd ed. (Cambridge Univ. Press, 1992).

P. W. Milonni, The quantum vacuum: an introduction to quantum electronics (Academic Press, 1994).

### Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

### Figures (2)

Fig. 1 Flow diagram for split-step FFT propagation using alternating direction implicit integration (ADI) or steady-periodic Greens function method for temperature calculation.
Fig. 2 Power in LP11 versus time (T) and distance (Z) at the input end of the amplifier similar to the one specified in [2]. This illustrates the growth of the low power mode over one beat cycle and over one beat length. The pump is unmodulated in this example.

### Tables (1)

Table 1 Model input parameters

### Equations (54)

$∂ E s ( x , y , z , t ) ∂ z = i 2 k c ∇ ⊥ 2 E s ( x , y , z , t ) + i [ k 2 ( x , y , z , t ) − k c 2 ] 2 k c E s ( x , y , z , t )$
$n ( x , y , z , t ) = n core ( x , y , z , t ) + d n d T T ( x , y , z , t ) ,$
$∂ E ( x , y , z , t ) ∂ z = i 2 k c ∇ ⊥ 2 E ( x , y , z , t ) + i [ k 2 ( x , y , z , t ) − k c 2 ] 2 k c E ( x , y , z , t ) + g ( x , y , z , t ) E ( x , y , z , t ) .$
$∂ E s ( x , y , z , t ) ∂ z = i 2 k c ∇ ⊥ 2 E s ( x , y , z , t ) .$
$∂ E s ( x , y , z , t ) ∂ z = i [ k 2 ( x , y , z , t ) − k c 2 ] 2 k c E s ( x , y , z , t ) + g ( x , y , z , t ) E s ( x , y , z , t ) .$
$Δ x = L x / ( N x − 1 )$
$Δ y = L y / ( N y − 1 )$
$Δ t = ϒ / N t .$
$ϕ μ ν = 4 L x L y sin ( π L x μ x ) sin ( π L y ν y )$
$ψ ( x , y , z ) = e − i β z ∑ μ = 1 ∞ ∑ ν = 1 ∞ C μ ν ϕ μ ν ( x , y )$
$∂ 2 ψ ∂ x 2 + ∂ 2 ψ ∂ y 2 + ∂ 2 ψ ∂ z 2 + n 2 ( x , y ) k ∘ 2 ψ = 0$
$∑ μ = 1 M x ∑ ν = 1 M y C μ ν [ n 2 ( x , y ) k ∘ 2 − β 2 − π 2 ( μ 2 L x 2 + ν 2 L y 2 ) ] ϕ μ ν ( x , y ) = 0.$
$∑ μ = 1 M x ∑ ν = 1 M y A μ ′ ν ′ , μ ν C μ ν = ( n eff 2 − n clad 2 ) C μ ′ ν ′ .$
$A μ ′ ν ′ , μ ν = ∫ 0 L x ∫ o L y d x d y { [ n 2 ( x , y ) − n clad 2 ] ϕ μ ′ ν ′ ( x , y ) ϕ μ ν ( x , y ) } − π 2 k ∘ 2 ( μ 2 L x 2 + ν 2 L y 2 ) δ μ , μ ′ δ ν , ν ′ .$
$A C = n ¯ 2 C$
$n eff ( m , n ) = n ¯ m , n 2 + n clad 2 .$
$n ( x , y ) = n clad + ( n core − n clad ) × exp ( − ln 2 { [ 2 ( x − x 0 ) / d core ] 2 + [ 2 ( y − y 0 ) / d core ] 2 } S ) ,$
$n bend ( x , y ) = n ( x , y ) ( x − x 0 ) / R bend .$
$c ε 0 2 ∫ d x d y n ( x , y ) [ N m , n ψ m , n ( x , y ) ] 2 = 1 W$
$u m , n ( x , y ) = N m , n ψ m , n ( x , y ) .$
$ρ C ∂ T ∂ t = Q + K ( ∂ 2 T ∂ x 2 + ∂ 2 T ∂ y 2 )$
$G ( x , y , ω | x ′ , y ′ ) = ∑ n = 0 ∞ F n ( y , y ′ ) P n ( x , x ′ , ω )$
$F n ( y , y ′ ) = 1 W K sin ( n π W y ) sin ( n π W y ′ )$
$P n ( x , x ′ , ω ) = { exp [ − σ n ( 2 H − | x − x ′ | ) ] − exp [ − σ n ( 2 H − x − x ′ ) ] + exp [ − σ n | x − x ′ | ] − exp [ − σ n ( x + x ′ ) ] } / σ n ( 1 − exp [ − 2 σ n H ] )$
$σ n 2 = ( n π W ) 2 + i ω ρ C K ,$
$T x , y t + Δ t / 2 − T x , y t = ( T x + Δ x , y t + Δ t / 2 − 2 T x , y t + Δ t / 2 + T x − Δ x , y t + Δ t / 2 ) λ x 2 + ( T x , y + Δ y t − 2 T x , y t + T x , y − Δ y t ) λ y 2 + Δ t 2 ρ C Q t + Δ t / 4$
$T t + Δ t / 2 A x = B y T t + Δ t 2 ρ C Q t + Δ t / 4$
$T t + Δ t / 2 = B y T t A x − 1 + Δ t 2 ρ C Q t + Δ t / 4 A x − 1 ,$
$A x = [ ( 1 + λ x ) − λ x / 2 0 ⋯ 0 0 0 − λ x / 2 ( 1 + λ x ) − λ x / 2 ⋯ 0 0 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ 0 0 0 ⋯ − λ x / 2 ( 1 + λ x ) − λ x / 2 0 0 0 ⋯ 0 − λ x / 2 ( 1 + λ x ) ] ,$
$B y = [ ( 1 − λ y ) λ y / 2 0 ⋯ 0 0 0 λ y / 2 ( 1 − λ y ) λ y / 2 ⋯ 0 0 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ 0 0 0 ⋯ λ y / 2 ( 1 − λ y ) λ y / 2 0 0 0 ⋯ 0 λ y / 2 ( 1 − λ y ) ] .$
$λ x = K Δ t ρ C Δ x 2$
$λ y = K Δ t ρ C Δ y 2 .$
$T x , y t + Δ t − T x , y t + Δ t / 2 = ( T x + Δ x , y t + Δ t / 2 − 2 T x , y t + Δ t / 2 + T x − Δ x , y t + Δ t / 2 ) λ x 2 + ( T x , y + Δ y t + Δ t − 2 T x , y t + Δ t + T x , y − Δ y t + Δ t ) λ y 2 + Δ t 2 ρ C Q t + 3 Δ t / 4 .$
$A y T t + Δ t = T t + Δ t / 2 B x + Δ t 2 ρ C Q t + 3 Δ t / 4 T t + Δ t = A y − 1 T t + Δ t / 2 B x + Δ t 2 ρ C A y − 1 Q t + 3 Δ t / 4 ,$
$T t + Δ t = A y − 1 B y T t A x − 1 B x + Δ t 2 ρ C A y − 1 ( Q t + Δ t / 4 A x − 1 B x + Q t + 3 Δ t / 4 ) ,$
$E s ( x , y , z , t ) = 1 2 π ∫ − ∞ ∞ ∫ − ∞ ∞ E s ( k x , k y , z , t ) e i k x x e i k y y d k x d k y .$
$∂ E s ( k x , k y , z , t ) ∂ z = − i k x 2 2 k c E s ( k x , k y , z , t ) − i k y 2 2 k c E s ( k x , k y , z , t ) .$
$ϕ ( k x , k y ) = − Δ z 2 ( k x 2 + k y 2 ) 2 k c .$
$E s ( x , y , z , t ) | z = 0 = ∑ modes P m , n s u m , n ( x , y ) M m , n s ( t )$
$n u ( x , y , t ) = P p σ p a / h ν p A p + I s σ s a / h ν s P p ( σ p a + σ p e ) / h ν p A p + I s ( σ s a + σ s e ) / h ν s + 1 / τ$
$∂ E s ( x , y , t ) ∂ z = 1 2 [ − σ s a + ( σ s a + σ s e ) n u ( x , y , t ) ] N Yb ( x , y ) E s ( x , y , t ) − 1 2 α ( x , y ) E s ( x , y , t ) ,$
$P p = P f p + P b p .$
$d P p d z = P p A p ∬ [ ( σ p a + σ p e ) n u ( x , y ) − σ p a ] N Yb ( x , y ) d x d y .$
$d P f p d z = d P p d z P f p P p$
$d P b p d z = − d P p d z P b p P p .$
$Q ( x , y , t ) = N Yb ( x , y ) [ ν p − ν s ν p ] [ σ p a − ( σ p a + σ p e ) n u ( x , y , t ) ] P p ( t ) A p + α ( x , y ) I s ( x , y , t ) ,$
$q ( x ′ , y ′ , ω m ) = Δ x Δ y ∑ i = 0 N t − 1 Q ( x ′ , y ′ , t i ) exp ( − i ω m t i ) .$
$T ( x , y , t ) = Real [ ∑ m = 0 max ∑ x ′ , y ′ q ( x ′ , y ′ , ω m ) G ( x , y , ω m | x ′ , y ′ ) exp ( i ω m t ) ] .$
$ϕ ( x , y , t ) = Δ z k 2 ( x , y , t ) − k c 2 2 k c .$
$ϕ ( x , y , t ) = Δ z ω c c ( [ n ( x , y , t ) − n clad ] + [ n ( x , y , t ) − n clad ] 2 2 n clad ) .$
$F m , n ( z , t ) = ∬ d x d y E s ( x , y , z , t ) u m , n ( x , y ) ∬ d x d y [ u m , n ( x , y ) ] 2 ,$
$P m , n ( z , t ) = | F m , n ( z , t ) | 2 .$
$A eff ( z , t ) = [ ∬ | E s ( x , y , z , t ) | 2 d x d y ] 2 ∬ | E s ( x , y , z , t ) | 4 d x d y .$
$P = h ν Δ ν = ( 6.6 × 10 − 34 ) ( 2.8 × 10 14 ) Δ ν = 1.8 × 10 − 19 Δ ν$

### Metrics

Select as filters

Cancel