EDFA的原理以及仿真

对于EDFA的设计，合理的EDFA仿真是必不可少的。本文通过讲解EDFA的原理，再延伸到对应的仿真模型建立，仿真模型代码，和大家一起摸索和分享。主要参考的书籍Springer出版社出版，P.C.Becker et.al.,的《Erbium-Doped Fiber Amplifiers：Funcamentals and Technology》一书。并且添加了对某些公式的推导，以及对书中一些阐述的个人更直白些的理解，学习过程中也发现书籍和其参考论文描述不一致的地方，做了一些内容上的合并。

Erbium

回顾一下元素电子分布原则，（参考一个YouTube非常优秀的视频 电子分布规则理论解释）：

1. Aufbau准则：电子应该先填充于低能级Orbital
2. 泡利不相容原则：每个Orbital只能有2个电子，而且这两个电子的自旋正好相反。
3. Hünd原则：电子应该先依次分配给同一个能级的尚没有分配电子的Orbital，然后再配对。

根据上述原则铒的元素周期表原子序数为68。电子分布为
$1s^2-2s^2-2p^6-3s^2-3p^6-4s^2-3d^{10}-4p^6-5s^2-4d^{10}-5p^6-6s^2-4f^{12}$
作为镧系稀土元素的一个特征，铒的4f亚层的轨道离原子核的半径要小于5s和5p，由于被充满电子的外层的5s和5p所形成的法拉第笼所保护，所以4f层电子不容易受到环境的影响（譬如如原子碰撞等）。同时4f电子不能参与成键，一般镧系稀土元素是通过丢失6s层的2个电子以及5p层的1个电子从而呈现3+的化学特征，所以镧系稀土元素的化学性质接近，但是4f的电子个数不同，其物理特性很不一样。

考虑到轨道半径，铒元素也可以缩写为 $(Xe)(4f{)}^{12}(6s{)}^2$其中$Xe$为元素氙（Xian，原子序号54）。这个标记形式很好的体现了轨道半径的特殊性。

4f层的电子由于原子力以及晶体场的影响会被拉扯出非常多的能级。首先原子力会把每个4f层可能的电子orbital分裂成$2S+1_{L_J}$能级。然后这些个能级又因为晶体场的影响分裂成多个子能级称为Stark level。正是因为这些能级的数量，导致镧系稀土元素存在有很多离子迁移的能级差，从而有丰富的光子释放吸收的波长范围。

离子从激发态自发迁移到基态的过程会释放出光子，这个过程称为自发辐射，且可以通过Einstein 参数A来进行描述($d{N}_2/dt=-A_{21}{N}_2$，自发辐射会导致激发态能级的population有一个指数级的衰减$N_2(t)=N_2(0)e^{-A_{21}t}$)。在实验中可以观测到这些光子的频率并不是一致的而是会在某一个中心频率周边有一定的概率分布。这样呈现出的光谱，就会有一定的展宽。对于展宽的物理参数就称之为线宽。这个展宽对于铒离子而言包含了不同的物理机制。主要可分为同类展宽（homogeneous broadening）以及非同类展宽（inhomogeneous broadening）。同类展宽意味着不管铒离子在晶体中的位置在哪里，其展宽的概率分布乃至中心频率都是一致的。非同类展宽则和铒离子所在的晶体环境有一定的关系。

典型的属于同类展宽的物理机制是lifetime broadening。根据海森堡测不准原理，一个激发态的lifetime $\Delta t$和该激发态的能级不可能同时确定。由于lifetime是有限的，所以根据如下关系：
$$\Delta E\Delta t>\frac{\hslash }{2}=>\Delta \omega \Delta t>\frac{1}{2}$$
势必会有一定的频率展宽而且该展宽大于$1/2\Delta t$。lifetime对于特定能级来说是一致的，和铒离子在晶体中的位置没有直接关系。所以这个展宽是属于同类展宽。

典型的非同类展宽的物理机制有crystal broadening。对于镧系稀土元素，虽然4f层电子不容易受到碰撞等影响，但是在晶体状态下，会收到晶体应变的影响。对于掺铒光纤而言，铒离子是掺杂在$Si{O}_2$的晶体之中，其晶体的内部原子的排列紧密或者稀疏程度是不同的。这个直接导致不同位置的铒离子会受到不同的晶体应变影响（对铒离子的电子云产生不同的影响）。晶体应变会导致能级的分裂从而生成Stark level。这个现象就会导致晶体中不同位置被激发的铒离子在回到基态上所释放的光子会有不同的频率以及线宽。其整体效果如下图所示：

如图a所示，位于不同晶体位置的铒离子集合释放出来的光子中心频率和线宽会不同，有的中心频率会有一些偏移。所有铒离子的频谱叠加在一起后就可以得到一个很宽的非同类展宽的谱型。这样一个现象也就意味着如果有一束超窄线宽激光射入到掺铒光纤中，不是所有的铒离子对于这个激光都会产生反应。只有部分和这个激光频率相同的铒离子会有相互作用。这个过程暗示了掺饵光纤对于不同频率的入射光在增益以及噪声系数上可能会有一些本地效应。可以预见这部分效应的仿真会使得模型进一步复杂。

三能级速率方程系统

一个维度的情况

Nomenclature

Symbol	Description	Dimension
${\varphi }_x$	入射光强通量	$Photons/s/c{m}^2$
${\sigma }_x$	Effective cross-sectional area	$c{m}^2$
$I$	光强	$W/c{m}^2$
$N$	Total population	$ions$

假设前提为pump和signal的强度以及铒离子分布在光纤截面上是常数（先不考虑在圆柱坐标的r方向上的分布变化），从而来简化问题。铒离子的能级状态跃迁主要有如下几个途径：

1. 由于吸收了入射光场的光子
2. 自发辐射以及受激辐射
3. 其他的一些能级跃迁途径（遗留问题1）

结论一：population的变化率是和当前能级的population以及单位面积通量和截面积乘积成正比。

对于此处的截面积需要做一个额外的解释。此截面积（cross-section）是表征了离子在两个能级之间迁移的概率。常用单位为$c{m}^2$。可见其也是一个面积的量纲。直观的理解可以认为当存在入射光时候，该截面积会捕获或者释放光子来影响入射光的通量(注意不是单位面积通量）

能级1和能级3能级差对应的频率直接决定了pump光源的频率。同时定义pump光源的单位面积上的入射光子通量（通量定义为能量场在截面dS的法向量投影与截面的乘积，然后再对dS做积分，但是注意此处为单位面积。）为${\varphi }_p$。通量意味着能级跃迁和pump光源的光强是有关系的。因为一维场景假设，如果能级1到能级3的光子吸收截面积可以定义为${\sigma }_p$，则pump光源的光子通量为${\sigma }_p{\varphi }_p$（量纲为：photons/unit time也就是单位时间内有多少个光子通过某个截面。如果乘以单个光子的能量hv，则得到单位时间内光子的总能量）

类似的能级1和能级2间的能级差对应的频率直接决定了可以被放大的信号的频率范围。其单位面积上的光子通量可以定义为${\varphi }_s$，同样由于一维假设，能级2和能级1的光子吸收、激发截面积可以定义为${\sigma }_s$. 信号光源的光子通量为${\sigma }_s{\varphi }_s$.

必须要强调一下，上面的描述基于一个假设也就是吸收和辐射的cross section是对称的。如果不对称则根据对应的是吸收还是辐射过程乘以相应的吸收cross section ${\sigma }_x^{(a)}$ 或者辐射cross section ${\sigma }_x^{(e)}$，这个在之后阐述的泛化二能级速率系统有所体现

现在考虑三能级速率方程的构建。第三能级的population变化率取决于如下场景，同时根据结论一我们可以得到：

1. 铒离子从第三能级通过辐射和非辐射方式迁移到第二能级的population。第三能级变化率为负：$-{\Gamma }_{32}{N}_3$
2. 铒离子由于吸收了pump光子而从第一能级跃迁到了第三能级。第三能级的变化率是正的: $N_1{\varphi }_p{\sigma }_p$
3. 铒离子由于pump入射光导致受激辐射，失去能量后从第三能级跃迁到了第一能级。第三能级变化率为负: $-N_3{\varphi }_p{\sigma }_p$
   从而对于第三级的population速率方程可以构建为$$\frac{d{N}_3}{dt}=-{\Gamma }_{32}{N}_3+(N_1-N_3){\varphi }_p{\sigma }_p$$

第二能级的population变化率结合结论一可以得到：

1. 铒离子从第二能级通过辐射方式迁移到第一能级的population。第二能级变化率为负：$-{\Gamma }_{21}{N}_2$
2. 铒离子从第三能级通过辐射和非辐射方式迁移到第二能级的population。第二能级变化率为正：${\Gamma }_{32}{N}_3$
3. 铒离子从第一能级由于吸收了信号光功率而跃迁到第二能级。第二能级变化率为正：$N_1{\varphi }_s{\sigma }_s$
4. 铒离子从第二能级由于信号光功率而受激辐射光子，从而失去能量跃迁到第一能级。第二能级变化率为负：$-N_2{\varphi }_s{\sigma }_s$
   从而对于第二级的population速率方程可以构建为$$\frac{d{N}_2}{dt}=-{\Gamma }_{21}{N}_2+{\Gamma }_{32}{N}_3+(N_1-N_2){\varphi }_s{\sigma }_s$$

第一能级的population变化率结合结论一可以得到：

1. 铒离子从第二能级通过辐射方式迁移到第一能级的population。第一能级变化率为正：${\Gamma }_{21}{N}_2$
2. 铒离子从第一能级由于吸收了信号光功率而跃迁到第二能级。第一能级变化率为负：$-N_1{\varphi }_s{\sigma }_s$
3. 铒离子从第二能级由于信号光功率而受激辐射光子，从而失去能量跃迁到第一能级。第一能级变化率为正：$N_2{\varphi }_s{\sigma }_s$
4. 铒离子由于pump入射光导致受激辐射，失去能量后从第三能级跃迁到了第一能级。第一能级变化率为正: $N_3{\varphi }_p{\sigma }_p$
5. 铒离子由于吸收了pump光子而从第一能级跃迁到了第三能级。第一能级的变化率是负的: $-N_1{\varphi }_p{\sigma }_p$
   从而对于第一级的population速率方程可以构建为$$\frac{d{N}_1}{dt}={\Gamma }_{21}{N}_2+(N_2-N_1){\varphi }_s{\sigma }_s+(N_3-N_1){\varphi }_p{\sigma }_p$$

当对于速率方程的求解，能级的稳态是描述正常工作的很重要的一个场景。稳态顾名思义即各能级变化率为0。通过求解如下的方程组
$$\begin{array}{rl}\frac{d{N}_1}{dt}& =\frac{d{N}_2}{dt}=\frac{d{N}_3}{dt}=0\\ N& =N_1+N_2+N_3\end{array}$$
可以得到以下关系：
$$\begin{array}{rl}{N}_2-N_1& =\frac{{\sigma }_p{\varphi }_p-{\Gamma }_{21}}{{\Gamma }_{21}+2{\sigma }_s{\varphi }_s+{\varphi }_p{\sigma }_p}N\\ {\varphi }_{th}& =\frac{{\Gamma }_{21}}{{\sigma }_p}=\frac{1}{{\tau }_2{\sigma }_p}\\ I_{th}& =\frac{h{\omega }_p{\Gamma }_{21}}{{\sigma }_p}=\frac{h{\omega }_p}{{\sigma }_p{\tau }_2}\end{array}$$
为了可以实现光的放大，我们所需要的一些关键信息有：

1. 第三能级population的存在时间可以短些，这样主要的population都在第一能级和第二能级上
2. 需要确定实现population inversion也就是$N_2-N_1$为正值的条件。当population inversion实现的时候，更多的铒离子会从第二能级通过辐射迁移到第一能级，从而实现信号放大。反之则会因为吸收信号光能量从而导致信号的衰减。
3. 使用如何的pump光功率才可以实现信号放大增益。

上述问题根据pump阈值光强$I_{th}$的表达式，在pump频率确定的情况下，可以得到如下结论：

1. 为了让所需的pump光强尽可能小，针对pump的effective cross section ${\sigma }_p$要大
2. 亚稳态能级（第二能级）的population维持存在时间要长

第二点对于在光纤中的铒离子而言是正好契合。其在亚稳态能级的存在时间可以达到10ms级。从而所需要的实现population inversion的pump阈值功率非常的小(0.5mW)级别。
同时，实验证明第三能级离子的lifetime是很短的，在10us级别。

小信号增益

Nomenclature

Symbol	Description	Dimension
${\varphi }_x$	光子通量	$Photons/s/c{m}^2$
$N$	Population密度	$ions/c{m}^3$
$I$	光强	$W/c{m}^2$

入光功率的变化是由于铒离子与入射光子间的相互作用而导致的。铒离子对于光子的吸收会导致入光功率的衰减。铒离子由于入射光子的激发从亚稳态跃迁到基态时候释放出的额外一个相干光子，这个受激辐射对于入光功率而言起到放大作用。该过程可以由如下公式描述。其中$dV$为体积元等于${\sigma }_{s}dz$。
$$\begin{array}{rl}dP& =P_{em}-P_{abs}\\ dP& =N_{2}dV{\sigma }_s{I}_s-N_{1}dV{\sigma }_s{I}_s\\ dP& =N_2{\sigma }_{s}dz{\sigma }_s{I}_s-N_1{\sigma }_{s}dz{\sigma }_s{I}_s\\ \frac{d\frac{P}{{\sigma }_s}}{dz}& =N_2{\sigma }_s{I}_s-N_1{\sigma }_s{I}_s\\ \frac{d{I}_s}{dz}& =(N_2-N_1){\sigma }_s{I}_s\end{array}$$
同时已知${\varphi }_s=I_s/h\omega$，从而也可以得到如下关系。
$$\frac{d{\varphi }_s}{dz}=(N_2-N_1){\sigma }_s{\varphi }_s$$
类似的，对于pump光源：
$$\frac{d{\varphi }_p}{dz}=(N_3-N_1){\sigma }_p{\varphi }_p$$
考虑稳态时得到的$N_2-N_1$，可以得到信号以及pump光强随着传播距离变化的微分方程：
$$\begin{array}{rl}\frac{d{I}_s}{dz}& =\frac{\frac{I_p{\varphi }_p}{h{\omega }_p}-{\Gamma }_{21}}{{\Gamma }_{21}+2\frac{I_s{\sigma }_s}{h{\omega }_s}+\frac{I_p{\sigma }_p}{h{\omega }_p}}N{\sigma }_s{I}_s\\ & =\frac{1-{\Gamma }_{21}^{\prime }}{{\Gamma }_{21}^{\prime }+2\frac{I_s{\sigma }_s{\omega }_p}{I_p{\sigma }_p{\omega }_s}+1}N{\sigma }_s{I}_s\end{array}$$
考虑小信号模型，则$I_s<<I_p$, 所以分母第2项趋向于0，从而得到
$$\frac{d{I}_s}{dz}=\frac{1-{\Gamma }_{21}^{\prime }}{1+{\Gamma }_{21}^{\prime }}N{\sigma }_s{I}_s$$
其中${\Gamma }_{21}^{\prime }=\frac{h{\omega }_p}{I_p{\varphi }_p{\tau }_{21}}=\frac{I_{th}}{I_p}=\frac{1}{I_p^{\prime }}$，可以得到
$$\frac{d{I}_s}{dz}=\frac{I_p^{\prime }-1}{I_p^{\prime }+1}N{\sigma }_s{I}_s$$
此方程求解可得
$$I_s(z)=I_s(0)e^{{\alpha }_{p}z}$$其中
$${\alpha }_p=\frac{I_p^{\prime }-1}{I_p^{\prime }+1}N{\sigma }_s$$
针对上述公式，可以看到

1. pump功率大于阈值功率才能出现增益
2. 如果pump功率远远大于阈值功率，则${\alpha }_p=N{\sigma }_s$。信号增益只和铒离子密度和信号的等效cross section相关。

饱和区间

当信号光功率过大的时候，可能出现如下情况：

1. 铒离子高能级population数即使全部翻转但是通过受激辐射生成的光子不足以提供相对于信号而言足够的功率增益。
2. 提高pump功率，是可以提高最大饱和输出功率的。

最优EDF长度

最优EDF长度的确定是需要和最后优化的目标相关的。

1. 最大化小信号增益：EDF最优长度位于pump功率较低到功率阈值时候的长度。从而保证pump光都用于提供了正向增益。
2. 增加pump功率变化下，增益的稳定度：可以考虑让EDF工作在饱和区域，从而pump的变化不影响增益。
3. 优化NF：TBC
4. 优化最大饱和输出功率：TBC

光纤横截面光功率与EDF的重叠区域等效

对于光放大器而言，在光纤中传播的模式和铒离子的重叠区域是非常关键的参数。该参数可以将pump，信号或者ASE噪声的功率和用于激发铒离子的有效光强以及掺铒光纤的铒离子分布更好的联系起来。

假设单模掺铒光纤而言，传输的光纤模式为$L{P}_{01}$基模，其能量在单模光纤界面上的分布可以近似为高斯分布，同时掺铒光纤的横截面上铒离子浓度是和横截面的直径以及方位角有关即$n(r,\varphi )$。进一步的模型简化可以考虑认为离子浓度只和横截面的直径有关系而和方位角无关$n(r,\varphi )->n(r)$，这个意味着离光纤中心任意半径的环状区域具有相同的铒离子掺杂浓度。

为了简化分析，可以考虑一个与实际掺杂浓度等效的模型，也就是认为在某一个半径$R$之内铒离子浓度相同，而且在整个R内的区域的铒离子数和实际掺杂的铒离子数相同，这个可以由如下公式表达(同时考虑和方位角无关）：
$${\int }_0^{2\pi }{\int }_0^{\infty }n(r,\varphi )rdrd\varphi =2\pi {\int }_0^{\infty }n(r)rdr=N\pi R^2$$
从这个公式可以看到，$N{R}^2$可以确定，但是N和R的值还需要单独定义。R值的定义如下
$$\pi R^2=2\pi {\int }_0^{\infty }\frac{n(r)}{n(0)}rdr$$
其中$n(0)$为光纤纤芯处的铒离子浓度。此处笔者认为其公式的含义是很明确，是相对于$n(0)$做了归一化后的权重做了积分，而且应该是包含了$n(0)$是最大掺杂浓度的强假设。为何如此定义，一个出发点可能如下：在归一化后，我们需要向一个高度为1，底面面积为$\pi r^2$的圆柱体中进行填充，如果高度需要超过1，则进一步扩大r直到正好等于$R$，这样最后的离子浓度等效就可以用一个常数$N$来进行灵活缩放。

所以$R$可以如下计算得到:
$$R=(2{\int }_0^{\infty }\frac{n(r)}{n(0)}rdr{)}^{-0.5}$$
以及
$$N=\frac{2{\int }_0^{\infty }n(r)rdr}{R^2}$$
如果考虑一个离子吸收光子从而进入激发态的过程，则在光纤截面的某一个位置$(r,\varphi )$，降低的pump功率是${\sigma }_p^{(a)}{I}_p(r,\varphi )n(r,\varphi )$，考虑整个截面的话可以得到：
$$\begin{array}{rl}\Delta P_p& ={\sigma }_p^{(a)}{\int }_0^{2\pi }{\int }_0^{\infty }{I}_p(r,\varphi )n(r,\varphi )rdrd\varphi \\ & ={\sigma }_p^{(a)}{\int }_0^{2\pi }{\int }_0^{\infty }{P}_p{I}^{(n)}(r,\varphi )n(r,\varphi )rdrd\varphi \\ & =N{\sigma }_p^{(a)}{P}_p{\int }_0^{2\pi }{\int }_0^{\infty }{I}^{(n)}(r,\varphi )\frac{n(r,\varphi )}{N}rdrd\varphi \end{array}$$
$I^{(n)}(r,\varphi )$为归一化的入射pump截面上某点光强，且${\int }_0^{2\pi }{\int }_0^{\infty }{I}^{(n)}(r,\varphi )rdrd\varphi =1$。至此，我们可以得到重叠系数的一般表达式为：
$$\Gamma ={\int }_0^{2\pi }{\int }_0^{\infty }{I}^{(n)}(r,\varphi )\frac{n(r,\varphi )}{N}rdrd\varphi$$

如果考虑与方位角无关以及截面积铒离子浓度等效模型则可以得到：
$$\begin{array}{rl}\Delta P_p& =2\pi N{\sigma }_p^{(a)}{P}_p{\int }_0^{\infty }{I}^{(n)}(r)\frac{n(r)}{N}rdr\\ & =N{\sigma }_p^{(a)}{P}_{p}2\pi {\int }_0^R{I}^{(n)}(r)rdr\end{array}$$
然后重叠系数可以简化为$$\Gamma =2\pi {\int }_0^R{I}^{(n)}(r)rdr$$其中$2\pi {\int }_0^{\infty }{I}^{(n)}(r)rdr=1$

对于在单模光纤中可以传播的$L{P}_{01}$基模，其$$I^{(n)}(r)=\frac{e^{-r^2/A_{spot}^2}}{\pi A_{spot}^2}$$
$A_{spot}$为spot大小，其数值可以由如下公式得到
$$A_{spot}=\frac{1}{\sqrt{2}}(0.65+\frac{1.619}{V^{1.5}}+\frac{2.879}{V^6})$$
所以对于阶跃光纤的基模而言，其重叠系数为$${\Gamma }^{(01)}=1-e^{-R^2/A_{spot}^2}$$
在定义了重叠系数的表达式之后，无论那种入射光其会和掺杂的铒离子产生互相作用的光强就可以简洁的通过下述表达式得到：
$$I(z)=\frac{P(z)\Gamma }{\pi R^2}=\frac{P(z)\Gamma }{A_{eff}}$$

三能级速率系统简化为二能级速率系统

由于铒离子在第三能级的lifetime非常的短，所以实际上三能级速率系统可以考虑简化为而能级速率。也就是第三能级的population为0，而本来跃迁到第三能级的铒离子数可以全部归于第二亚稳态上的population。这个假设适合于pump波长为980nm的情况。基于上述讨论，原三能级速率系统考虑$N_3=0$，而且也不存在pump photon emission ${\sigma }_p^{(e)}=0$，然后$N_1{\sigma }_p^{(a)}{\varphi }_p$可以直接归算于$N_2$ population的增加。

另外一种情况是当pump波长为1480nm的情况。在这个pump波长下，铒离子跃迁到的第三能级和第二亚稳态实际上是有重叠的。更准确的描述是，第三能级和第二能级由于热分布的原因，会有不同的子能级（sublevel）。子能级的集合叫multiplet。在1480nm下基态的铒离子会吸收pump光子从而跃迁到所谓的第三能级。而这个第三能级属于第二亚稳态的一个较高能量的子能级。如下图所示：

所以1480nm pump下的铒离子跃迁，也可以类比于980nm下明显分离的三能级系统。但是其实1480nm pump下，其本身可以通过二能级系统来表征，即可以认为第三能级的population即为$N_2$的一部分。对于1480nm pump，是会通过受激辐射产生pump频率的光子，也就是${\sigma }_p^{(e)}\ne 0$。

泛化二能级速率系统

基于上述讨论，能级模型可以统一为如下图示。

population的速率方程可以和推导三能级速率方程一样的方法得到：
$$\begin{array}{rl}\frac{d{N}_2}{dt}& =-{\Gamma }_{21}{N}_2+N_1{\varphi }_p{\sigma }_p^{(a)}-N_2{\varphi }_p{\sigma }_p^{(e)}+N_1{\varphi }_s{\sigma }_s^{(a)}-N_2{\varphi }_s{\sigma }_s^{(e)}\\ \frac{d{N}_1}{dt}& ={\Gamma }_{21}{N}_2-N_1{\varphi }_p{\sigma }_p^{(a)}+N_2{\varphi }_p{\sigma }_p^{(e)}-N_1{\varphi }_s{\sigma }_s^{(a)}+N_2{\varphi }_s{\sigma }_s^{(e)}\\ N& =N_1+N_2\end{array}$$
需要强调的是$N_x$依然为population density （$ions/c{m}^3$）。显而易见，这两个微分方程不是相互独立的。因为$$\frac{d{N}_1}{dt}=-\frac{d{N}_2}{dt}$$
考虑稳态则$\frac{d{N}_2}{dt}=0$可得
$$N_2(z)=\frac{\frac{{\tau }_{21}{\sigma }_p^{(a)}}{h{\omega }_p}{I}_p(z)+\frac{{\tau }_{21}{\sigma }_s^{(a)}}{h{\omega }_s}{I}_s(z)}{1+({\sigma }_p^{(a)}+{\sigma }_p^{(e)})\frac{{\tau }_{21}}{h{\omega }_p}{I}_p(z)+({\sigma }_s^{(a)}+{\sigma }_s^{(e)})\frac{{\tau }_{21}}{h{\omega }_s}{I}_s(z)}$$
也可以用同样的方式得到信号和pump光强的微分方程：
$$\begin{array}{rl}\frac{d{I}_p(z)}{dz}& =(N_2{\sigma }_p^{(e)}-N_1{\sigma }_p^{(a)})I_p(z)\\ \frac{d{I}_s(z)}{dz}& =(N_2{\sigma }_s^{(e)}-N_1{\sigma }_s^{(a)})I_s(z)\end{array}$$
我们更为关心的是能够产生信号增益的pump光强，所以令$\frac{d{I}_s(z)}{dz}>0$以及将上述$N_2(z)$代入，可以得到符合该目标的pump的阈值光强为
$$I_{th}=\frac{h{\omega }_p}{{\tau }_{21}}\frac{1}{{\sigma }_p^{(a)}\left(\frac{{\sigma }_s^{(e)}}{{\sigma }_s^{(a)}}\right)-{\sigma }_p^{(e)}}$$

对于多pump，多信道，上述$N_2(z)$的结论可以很容易地进行如下扩展：
$$N_2(z)=\frac{{\sum }_{pi}\frac{{\tau }_{21}{\sigma }_{pi}^{(a)}}{h{\omega }_{pi}}{I}_{pi}(z)+{\sum }_{si}\frac{{\tau }_{21}{\sigma }_{si}^{(a)}}{h{\omega }_{si}}{I}_{si}(z)}{1+{\sum }_{pi}({\sigma }_{pi}^{(a)}+{\sigma }_{pi}^{(e)})\frac{{\tau }_{21}}{h{\omega }_{pi}}{I}_{pi}(z)+{\sum }_{si}({\sigma }_{si}^{(a)}+{\sigma }_{si}^{(e)})\frac{{\tau }_{21}}{h{\omega }_{si}}{I}_{si}(z)}$$

放大自发辐射（ASE）

ASE的物理过程还是很有必要回顾一下的。受激发的离子会有一定概率自发从高能级迁移到基态，并且释放一个和信号不相干的光子。这个不相干的光子激发更多的高能级离子释放和其相干的光子。所以ASE其实也有受激辐射的过程在其中，只是生成的光子和信号的光子并不coherent。

和之前信号或者pump的速率方程是有地方需要区别考虑的。对于ASE，在距离z的地方，除了之前ASE信号功率，还有本地新的ASE噪声的生成。而且这个ASE噪声是和频率相关的，可能会产生新的频率的光子。所以对应的速率方程有如下几个部分：
$$\frac{d{P}_{ASE}(\omega )}{dz}=(N_2{\sigma }^{(e)}(\omega )-N_1{\sigma }^{(a)}(\omega ))P_{ASE}(\omega )+N_2{\sigma }^{(e)}(\omega )P_{ASE}^0(\omega )$$
其中$P_{ASE}^0(\omega )=2h\omega \Delta \omega$

考虑ASE的情况下对应的population equation 需要扩展为如下：
$$N_2(z)=\frac{{\sum }_{pi}\frac{{\tau }_{21}{\sigma }_{pi}^{(a)}}{h{\omega }_{pi}}{I}_{pi}(z)+{\sum }_{si}\frac{{\tau }_{21}{\sigma }_{si}^{(a)}}{h{\omega }_{si}}{I}_{si}(z)+{\sum }_{ni}\frac{{\tau }_{21}{\sigma }_{ni}^{(a)}}{h{\omega }_{ni}}{I}_{ni}(z)}{1+{\sum }_{pi}({\sigma }_{pi}^{(a)}+{\sigma }_{pi}^{(e)})\frac{{\tau }_{21}}{h{\omega }_{pi}}{I}_{pi}(z)+{\sum }_{si}({\sigma }_{si}^{(a)}+{\sigma }_{si}^{(e)})\frac{{\tau }_{21}}{h{\omega }_{si}}{I}_{si}(z)+{\sum }_{ni}({\sigma }_{ni}^{(a)}+{\sigma }_{ni}^{(e)})\frac{{\tau }_{21}}{h{\omega }_{ni}}{I}_{ni}(z)}$$
需要强调是的ni和si的集合也是和z相关的。考虑到可能在z有新的波长信道加入以及自发辐射引入的新噪声频率。

阶段性总结一

至此可以在假设了阶跃单模光纤的情况下，得到如下三组光场传播微分方程：
$$\begin{array}{rl}\frac{d{I}_p(z)}{dz}& =(N_2{\sigma }_p^{(e)}-N_1{\sigma }_p^{(a)})I_p(z)\\ \frac{d{I}_s(z)}{dz}& =(N_2{\sigma }_s^{(e)}-N_1{\sigma }_s^{(a)})I_s(z)\\ \frac{d{P}_{ASE}(\omega )}{dz}& =(N_2{\sigma }^{(e)}(\omega )-N_1{\sigma }^{(a)}(\omega ))P_{ASE}(\omega )+N_2{\sigma }^{(e)}(\omega )P_{ASE}^0(\omega )\end{array}$$
通过重叠系数，可以将光强和入射光的功率进一步相联系起来。
这一组微分方程将是分析pump，信号以及ASE噪声沿着z方向传播起到重要作用。
但是需要指出的是其

1. 假设所有的铒离子都是同类的，并不考虑掺铒光纤中的晶体应变导致的对于入射光子不同频率的不同反应。
2. 没有区分前向和反向ASE噪声
3. 没有考虑沿着z方向或者-z方向的光纤衰减

EDFA 模型公式

所述模型的前提条件

1. 所有离子都是同类展宽（homogeneous broadening）
2. 考虑等效截面。铒离子的分布在光纤截面中心到R的圆形范围内均匀分布。而且R小于铒纤的纤芯半径（铒离子没有掺杂到光纤包层，称为confined EDF）。

考虑光纤本身的衰耗，可以得到如下的一组微分方程，这组微分方程的求解即为EDFA模型建立的基础。
$$\begin{array}{rl}\frac{d{P}_p(z)}{dz}& =(N_2{\sigma }_p^{(e)}-N_1{\sigma }_p^{(a)}){\Gamma }_p{P}_p(z)-{\alpha }_p{P}_p(z)\\ \frac{d{P}_s(z)}{dz}& =(N_2{\sigma }_s^{(e)}-N_1{\sigma }_s^{(a)}){\Gamma }_s{P}_s(z)-{\alpha }_s{P}_s(z)\\ \frac{d{P}_{ASE}^{+}(v_j)}{dz}& =(N_2{\sigma }_{v_j}^{(e)}-N_1{\sigma }_{v_j}^{(a)})P_{ASE}^{+}(v_j)+N_2{\sigma }_{v_j}^{(e)}{\Gamma }_{s}h{v}_j\Delta v_j-{\alpha }_{v_j}{P}_A^{+}(v_j)\\ \frac{d{P}_{ASE}^{-}(v_j)}{dz}& =-(N_2{\sigma }_{v_j}^{(e)}-N_1{\sigma }_{v_j}^{(a)})P_{ASE}^{-}(v_j)-N_2{\sigma }_{v_j}^{(e)}{\Gamma }_{s}h{v}_j\Delta v_j+{\alpha }_{v_j}{P}_A^{-}(v_j)\end{array}$$

需要指出的是，由于ASE噪声也靠近1550nm，所以${\sigma }_{v_j}={\sigma }_s(v_j)$。前向和反向传播ASE的边界条件是非常清晰的：$P_A^{+}(z=0)=0$以及$P_A^{-}(z=L)=0$。当$\Delta v_j$取值较小的时候，可以认为${\sigma }_{v_j}$是个常数。

上述耦合微分方程的求解牵涉到dual boundary problem，可以使用Runge-Kutta method求数值解。以下篇章将用于详细阐述数值求解的过程。

Runge-Kutta Method
这里尽量给出读者一个直白的解释，然后附加经过笔者仔细审阅可以帮助读者尽快掌握的资源连接。
如果可以科学上网，那如下的资源是很值得推荐的：
港科大 Prof. Jeffery Chasnov: Numerical Method for Engineer 49~53
从简单的单个ODE数值求解出发，首先明确需要解决的问题如下面所述：

$\frac{dx}{dt}=f(t,x)$ with $x(t_0)=x_0$，需要知道的是$x(t_N)$时候的值。
基于数值求解，那就不需$x$的具体解析解。很多情况下，解析解也是很难求出的。我们需要找出$x(t_{n+1})=x(t_n)+Expression$ 这么一个递归形式的表达式， 从而可以从边界条件出发递归到所需要时间点时候x的值即$x(t_0)->x(t_N)$。为了简化标记$x_{n+1}=x(t_{n+1})$。

以下分别介绍几种数值解ODE的方法：

1. Euler Method
   如下图所示Euler方法是最为直观简单的一个方法，即用当前已知点$(t_n,x_n)$处的斜率来推算$x_{n+1}$，斜率自然就是$ f(t_n, x_n)$。所以
   $$x_{n+1}=x_n+\Delta t\cdot f(t_n,x_n)$$
   

由欧拉方法得到的$x_{n+1}$即为图中的红点。可以看到，当两个点斜率相差很大的时候，其估计值和实际值可能会出现很大的偏差，由$ϵ$表示。通过减小步长是可能提高数值解的准确度的，但是运算量肯定会提高。

2. Modified Euler Method
   改进欧拉方法的改进点在于既然直接用$(t_n,x_n)$处的斜率来推算会导致较大偏差，那不如同时考虑$(t_n,x_n)$和$(t_{n+1},x_{n+1})$处的斜率，且就用两个点的平均斜率来做运算，则递归表达式如下：
   $$x_{n+1}=x_n+\Delta t\cdot \frac{f(t_n,x_n)+f(t_{n+1},x_{n+1})}{2}$$
   但是很明显$f(t_{n+1},x_{n+1})$ 中$x_{n+1}$并不知道。为了解决这个问题，$x_{n+1}$是用和欧拉方法一样的办法通过$(t_n,x_n)$处的斜率来得到这个$x_{n+1}$预测值，这里记为${\hat{x}}_{n+1}$。这个方法被广泛称之为"predictor and corrector method"，所以迭代公式变为：
   $$x_{n+1}=x_n+\frac{\Delta t}{2}\cdot (f(t_n,x_n)+f(t_{n+1},x_n+\Delta t\cdot f(t_n,x_n)))$$
   
   现在将迭代公式做一下等效替换，令：
   $k_1=\Delta t\cdot f(t_n,x_n)$
   $k_2=\Delta t\cdot f(t_n+\Delta t,x_n+k_1)$
   $x_{n+1}=x_n+\frac{1}{2}\cdot (k_1+k_2)$
   综上，改进型欧拉方法本质上是考虑了$t_n$和$t_n+\Delta t$，两个时间点处的斜率，且对斜率做了加权平均后来对$x_{n+1}$做预估算。这样就引来两个问题：
   a. 是否可以考虑更多中间的时间点譬如$t_n+\frac{\Delta t}{2}$, $t_n+\frac{\Delta t}{3}$
   b. 如何对这些中间点得到的值在$t_n+\Delta t$进行额外的运算从而使得$x_{n+1}$的预测值比较准确。

3. Runge-Kutta Method
   为了解释上面的问题：在符合一定的限制条件情况下，可以使用不同的时间点和对应的不同的估计值来估计$x_{n+1}$。
   如果考虑使用2个点来做预估，则问题可以描述如下
   $k_1=\Delta t\cdot f(t_n,x_n)$
   $k_2=\Delta t\cdot f(t_n+\alpha \Delta t,x_n+\beta k_1)$
   $x_{n+1}=x_n+a\cdot k_1+b\cdot k_2$
   则
   $$x_{n+1}=x_n+a\Delta t\cdot f(t_n,x_n)+b\Delta t\cdot f(t_n+\alpha \Delta t,x_n+\beta k_1)$$
   这样只要选择合适的$\alpha ,\beta ,a,b$就可以进行ODE的数值求解。
   通过对$f(t_n+\alpha \Delta t,x_n+\beta \Delta t\cdot f(t_n,x_n))$ 做在$t_n$的一阶泰勒级数展开（需要用到全导数）。
   $$f(t_n+\alpha \Delta t,x_n+\beta k_1)=f(t_n,x_n)+\alpha \Delta t{f}_t(t_n,x_n)+\beta \Delta t{f}_x(t_n,x_n)f(t_n,x_n)$$
   从而得到：
   $$\begin{array}{r}{x}_{n+1}=x_n+(a+b)\Delta t\cdot f(t_n,x_n)+\alpha b\Delta t^2{f}_t(t_n,x_n)+\beta b\Delta t^2{f}_x(t_n,x_n)f(t_n,x_n)\end{array}$$
   现在我们再思考一个问题，记$x(t)$为ODE的解析解，已知$x_n$然后需要求得$x_{n+1}$。我们可以将$x(t)$在$t=t_n$处做泰勒展开（只到2阶），从而得到：
   $$\begin{array}{rl}{x}_{n+1}& =x_n+\frac{(t_{n+1}-t_n)}{1!}\frac{dx}{dt}+\frac{(t_{n+1}-t_n{)}^2}{2!}\frac{d^{2}x}{d{t}^2}\\ & =x_n+\Delta t\frac{dx}{dt}{|}_{t_n}+\Delta t^2\frac{d^{2}x}{d{t}^2}{|}_{t_n}\end{array}$$
   其中$\frac{dx}{dt}$ 就是$f(t,x(t))$，后面2阶通过全导数链式法则可以得到：
   $$\frac{d^{2}x}{d{t}^2}=\frac{df}{dt}=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial x}\cdot \frac{dx}{dt}=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial x}\cdot f(t,x)$$
   所以将两个不同路径得到的$x_{n+1}$展开式写在一起后，可以得到：
   $$\begin{array}{rl}{x}_{n+1}& =x_n+(a+b)\Delta t\cdot f(t_n,x_n)+\alpha b\Delta t^2{f}_t(t_n,x_n)+\beta b\Delta t^2{f}_x(t_n,x_n)f(t_n,x_n)\\ x_{n+1}& =x_n+\Delta t\cdot f(t_n,x_n)+\frac{\Delta t^2}{2}{f}_t(t_n,x_n)+\frac{\Delta t^2}{2}{f}_x(t_n,x_n)f(t_n,x_n)\end{array}$$
   至此相应的constraints也比较明确，也就是：
   $a+b=1$
   $\alpha b=1/2$
   $\beta b=1/2$
   只要取值符合上述的限制条件，就可以构造出2阶Runge Kutta Method族。
   可见Modified Euler Method是一个特例，也就是a=b=1/2, $\alpha =\beta =1$。
   另外一个特例为Midpoint Method, $\alpha =1/2,b=1,\beta =1/2,a=0$

常用的4阶Runge-Kutta Method：
RK4是使用比较广泛的方法，证明过程可以类比，这里描述其构造式。
$k_1=\Delta tf(t_n,x_n)$
$k_2=\Delta tf(t_n+\frac{1}{2}\Delta t,x_n+\frac{1}{2}{k}_1)$
$k_3=\Delta tf(t_n+\frac{1}{2}\Delta t,x_n+\frac{1}{2}{k}_2)$
$k_4=\Delta tf(t_n+\Delta t,x_n+k_3)$
$x_{n+1}=x_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)$

Runge-Kutta-Method数值求解EDF 双边界条件耦合ODE系统

代码附在资源里，请下载支持一下作者。
Solving EDFA coupled ODE using RK4
代码的功能如下：

1. 支持C-Band DWDM全频信号仿真
2. 支持C-Band ASE功率谱仿真
3. 铒纤关键参数可配，使用者可以根据实际情况调节

代码局限在于

1. 双边界问题求解。目前可以通过类似于shooting method去猜测反向ASE在z=0出的功率值，然后得到解后观察ASE在z=L处的值是否接近于0来判断所得结果是否收敛。
   系统性求解boundary value problem for coupled ODE会在后面继续解释。
2. EDF 非线性问题会在下文继续深入研究。