Optical amplifier fundamentals
                                                H. Q. Le (c) Copyrighted

Summary

Key concepts in optical amplification process:

1- The signal is amplified with gain as in the following equation:
                                    (( d I[z ])/(d z)) =g[z] I[z]          
but gain g[z] can be saturated:
                                  g[z]=  ((Subscript[g, 0])/(1+(( I(z) )/(Subscript[I, sat])))) 
where (Subscript[g, 0]) is a characteristic value, and (Subscript[I, sat]) , the saturation intensity is: 
                                    (Subscript[I, sat]) = ((Subscript[γ, spont])/(2 (Subscript[γ, stim]))) h ν 
where (Subscript[γ, spont]) and (Subscript[γ, stim]) are the amplifier stimulated emission and spontaneous emission coefficients, respectively.
The main features are:
              1.a power (or intensity) gain is exponential at small signal
              1.b but becomes linear when the power (or intensity) is large 

2- Just like a laser has certain spectral range, gain can occur only over certain range of wavelength that depends on the medium. The variation of gain vs. wavelength is called gain spectrum,  g[λ]

3- The light coming out of an OA is not just the amplified input signal (what if we don't have any input?) but also includes Amplified Spontaneous Emission (ASE): very important effect on the OA noise characteristics, the ASE power is:
                             β ((N/2) A (Subscript[γ, spont])h ν+(((Subscript[g, 0]) (Subscript[P, s ]))/(1+((P[x])/(Subscript[P, s])))))
4- There are 2 key noise terms of an OA:
      - Signal--spontaneous beat noise
      - Spontaneous--spontaneous beat noise
Both are critical to the performance of OA.

Introduction

Discussion

Optical amplifiers amplify an optical signal.
We can describe with the simple model:

Gain

Differential gain:
g[λ,x] = (Underscript[Lim, Δx->0])(((P[x+Δx])/(P[x] ))-1)(1/(Δx)) = (1/(P[x;λ])) (( d P[x;λ])/(d x)) = ((d Log[P[x;λ] ])/(d x))
Log[((Subscript[P, out])/(Subscript[P, in]))] =Log[((P[(Subscript[x, L])])/(P[(Subscript[x, 0])]))]  = (∫_0^L)g[λ,x]dx
Net gain:
((P[(Subscript[x, L])])/(P[(Subscript[x, 0])])) = (e^((∫_0^L)g[λ,x]dx))
When gain is very small,
((Subscript[P, out])/(Subscript[P, in])) =(e^((∫_0^L)g[λ,x]dx)) ≈ 1+(∫_0^L)g[λ,x]dx = 1+ G[λ]  
EDFA simple model of differential gain:
g[λ,x]=(SuperStar[g])[λ] (Subscript[N, 2])[z] -α[λ] (Subscript[N, 1])[z] 
where 
(SuperStar[g])[λ] = (Subscript[Γ, s]) (Subscript[n, Er]) (Subscript[σ, a])[λ]

Gain saturation

Gain saturation occurs when the differential gain coefficient g[λ,x] is also a function of optical power: g[λ,x] ->g[λ,x,P]  As P increase, g decreases. A case example is: 
g =((Subscript[g, 0])/(1+(I/(Subscript[I, s])))) where (Subscript[I, s]) is called saturation intensity.

DSolve[(1/(P[z])) (∂_z)P[z] == ((Subscript[g, 0])/(1+((P[z])/Ps))),P[z],z]

The solution of this differential equation is call ProductLog (it is a function): 
P(z)->(Subscript[P, s]) ProductLog(((e^(((Subscript[c, 1])/Ps)+(Subscript[g, 0]) z))/(Subscript[P, s])))

Plot[ProductLog[x], {x,0,10}]

ProductLog(((E^(((Subscript[c, 1])/Ps)+(Subscript[g, 0]) z))/(Subscript[P, s]))) (Subscript[P, s]) /. z->0

Solve[(Subscript[P, s]) W((X/(Subscript[P, s])))==Pin,X]

Hence:   P(z)=(Subscript[P, s]) ProductLog(((e^(((Subscript[c, 1])/Ps)+(Subscript[g, 0]) z))/(Subscript[P, s])))=(Subscript[P, s]) ProductLog(((Subscript[P, in])/(Subscript[P, s])) (e^((Subscript[P, in])/(Subscript[P, s])))(e^((Subscript[g, 0]) z)))

Plot[{  ProductLog((((E^((Pin/(Subscript[P, s]))+z (Subscript[g, 0]))) Pin)/(Subscript[P, s]))) (Subscript[P, s]) /. (Subscript[P, s])->1 /. Pin->0.1
	                                                          /. (Subscript[g, 0])->1,
		         ((5 *(Subscript[g, 0]))/(1+ProductLog((((E^((Pin/(Subscript[P, s]))+z (Subscript[g, 0]))) Pin)/(Subscript[P, s]))))) /. (Subscript[P, s])->1  /. Pin->0.1
	                                                           /. (Subscript[g, 0])->1}
	 ,{z,0,10}
	     , PlotStyle->{{RGBColor[1.,0.,0]},{RGBColor[0.,0.,1]}}
, Frame-> True, GridLines->Automatic
,FrameLabel-> {"Distance (a. u)","Power (red) - Gain (blue)"} 
		                  ]

Work note: for detailed derivation

(∂_z) P[z]((1/(P[z]))+(1/(Subscript[P, s])))==(Subscript[g, 0 ])

((1/(Subscript[P, s]))+(1/(P(z)))) P'(z)==∫(Subscript[g, 0])\[DifferentialD]z

∫((1/(Subscript[P, s]))+(1/(P(z)))) P'(z)\[DifferentialD]z==∫(Subscript[g, 0])\[DifferentialD]z

log(((P(z))/(Subscript[P, s])))+((P(z))/(Subscript[P, s]))==z (Subscript[g, 0]) +c /. z->0

log(((P(z))/(Subscript[P, s])))+((P(z))/(Subscript[P, s]))==z (Subscript[g, 0]) +c /. c->log(((P(0))/(Subscript[P, s])))+((P(0))/(Subscript[P, s]))

FullSimplify[log(((P(z))/(Subscript[P, s])))+((P(z))/(Subscript[P, s]))-log(((P(0))/(Subscript[P, s])))-z (Subscript[g, 0])-((P(0))/(Subscript[P, s]))]

log(((P(z))/(Subscript[P, in])))+((P(z)-(Subscript[P, in]))/(Subscript[P, s]))==z (Subscript[g, 0])

log(((P(z))/(Subscript[P, in])))+(((P(z))/(Subscript[P, in]))-1)((Subscript[P, in])/(Subscript[P, s]))==z (Subscript[g, 0])  /. ((Subscript[P, in])/(Subscript[P, s]))->ρ

Continued

As ((P(z))/(Subscript[P, in])) >> 1 ; the equation becomes: 
ρ (((P(z))/(Subscript[P, in]))-1)==z (Subscript[g, 0]) or ((P(z))/(Subscript[P, in]))~ ((Subscript[g, 0])/ρ)z+1: the gain becomes just linear, not exponential in z.
Net gain is determined by the equation:
log(G)+ρ (G-1)==L (Subscript[g, 0])
where ρ = ((Subscript[P, in])/(Subscript[P, s ])). The larger the input, the larger ρ is and the net gain is thus lower. But the power output is the relevant quantity.

NetGain[Pin_,Psat_, gain_,L_]:=
	G /. FindRoot[log(G)+(Pin/Psat)* (G-1)== L gain,
			                              {G, 1, Exp[gain*L]}]

Plot for increasing gain

Plot[ Table[ 
	        10*Log[10.,NetGain[10^(Pin/10.),1, g,5] ]
			                                             , {g,1,2,0.2}] 
			 ,{Pin,-50,0}, PlotRange->{{-50,0},All}
, Frame-> True, GridLines->Automatic
,FrameLabel-> {"Power in (dB)","Net Power Gain (dB)"} ]

log(G)+ρ (G-1)==L (Subscript[g, 0]) Let ρ -> 0:
G = (e^(L (Subscript[g, 0]))) 
large input power: 
((log(G))/ρ)+(G-1)==((L (Subscript[g, 0]))/(ρ )) becomes:
(G-1)==((L (Subscript[g, 0]))/(ρ )) or G=1+((L (Subscript[g, 0]))/(ρ ))= 1+((L (Subscript[g, 0]))/(Subscript[P, in ]))(Subscript[P, s]) 
 (Subscript[P, out]) = (Subscript[P, in]) + (Subscript[g, 0]) L (Subscript[P, s ]) This is additive gain: power out can only be additive

Gain saturation

What is gain saturation? We use to think the growth of a signal in an OA as: (e^(g L)) where g is constant. Is this correct all the time? It simply can't be. No power or energy can grow exponentially to infinity: energy conservation stipulates that the maximum energy one can get out of a finite segment L can not exceed the energy put into it. If the exponential term were correct, the power extracted would be:
P ((e^(g L)) - 1). If P is very large already and P ((e^(g L))-1)>> Pump energy, noway that this expression is correct. We have to re-examine the gain concept.

Consider a simple 2-level system.
The power out is:  
P[x+Δx] = P[x](1-α) + ((Subscript[N, 2])-(Subscript[N, 1])) A Δx (Subscript[γ, stim]) ((P[x])/(A  h ν)) h ν 
			         =P[x](1-α) + ((Subscript[N, 2])-(Subscript[N, 1])) Δx (Subscript[γ, stim]) P[x] 
where -α P[x] is the power loss due to other mechanism other than reabsorption 
Thus:  
(1/(P[x])) ((P[x+Δx] - P[x])/(Δx)) = ((Subscript[N, 2])-(Subscript[N, 1]))  (Subscript[γ, stim]) -α 
(1/(P[x]))d (P/(d x)) = ((Subscript[N, 2])-(Subscript[N, 1]))  (Subscript[γ, stim]) -α= g (gain) -α  
The 2-level system gain is:
g = ((Subscript[N, 2])-(Subscript[N, 1]))  (Subscript[γ, stim ])
If we solve for a constant gain, the solution (e^((g-α) x)) is not sustainable: no power can grow forever. As P[x] gets larger, P[x+Δx] - P[x]= ((Subscript[N, 2])-(Subscript[N, 1])) Δx (Subscript[γ, stim]) P[x] also gets larger, but if the energy rate (power) pumped into that volume is finite, the power gain can not exceed the input power! thus, this gain can not be forever. This leads to saturation. Saturation is a NOT a fundamental effect at the optical transition levels. It is a nonlinear effect, BUT NOT a fundamental nonlinear effect. It is a microscopic effect associated with population and pump dynamic. We can expect that as the power gets higher, the rate of population going from level 1 to level 2 is almost as much as the other way. If so, the population becomes nearly equal, i.e. (Subscript[N, 2])≈(Subscript[N, 1]), then the gain becomes very small.
Thus, the issue is: for a given pump rate, what is the population distribution? 
((d (Subscript[N, 2]))/(d t))=  (Subscript[N, pump]) - (Subscript[N, 2]) (Subscript[γ, spont.]) -(( ((Subscript[N, 2])-(Subscript[N, 1])) Δx (Subscript[γ, stim]) P[x])/(A Δx  h ν )) 
At steady state:  
((d (Subscript[N, 2]))/(d t))=0=  (Subscript[N, pump]) - (Subscript[N, 2]) (Subscript[γ, spont.]) - ((Subscript[N, 2])-(Subscript[N, 1]))  (Subscript[γ, stim]) ((P[x])/(A  h ν)) 
                =  (Subscript[N, pump]) - (Subscript[N, 2]) (Subscript[γ, spont.]) - ((Subscript[N, 2])-(Subscript[N, 1]))  (Subscript[γ, stim]) ((I[x])/(h ν))  
Also: (Subscript[N, 2])+(Subscript[N, 1]) = N  
So, here is the situation: Total (Subscript[N, 2])+(Subscript[N, 1 ]) is finite: there is ONLY that many gain elements in a given volume. One can only pump that much power into such a volume. As the intensity is very large, equation  0=  (Subscript[N, pump]) - (Subscript[N, 2]) (Subscript[γ, spont.]) - ((Subscript[N, 2])-(Subscript[N, 1]))  (Subscript[γ, stim]) ((I[x])/(h ν)) tells us that ((Subscript[N, 2])-(Subscript[N, 1])) has to becomes very small, since (Subscript[N, pump]) is finite. Thus, the gain ((Subscript[N, 2])-(Subscript[N, 1]))  (Subscript[γ, stim]) has to be also small.

Work note: by solving eqs., we get:

Solve[{0==  (Subscript[N, pump]) - (Subscript[N, 2]) (Subscript[γ, spont]) - ((Subscript[N, 2])-(Subscript[N, 1]))  (Subscript[γ, stim]) ((IP[x])/(h ν)),
		        (Subscript[N, 2])+(Subscript[N, 1]) == N},{(Subscript[N, 2]),(Subscript[N, 1]) }]

(Subscript[N, 2])-(Subscript[N, 1])  /. (Subscript[N, 2])->N+((N ((Subscript[γ, spont])+((IP(x) (Subscript[γ, stim]))/(h ν)))-(Subscript[N, pump]))/(-(Subscript[γ, spont])-((2 IP(x) (Subscript[γ, stim]))/(h ν))))
                   /. (Subscript[N, 1])->-((N ((Subscript[γ, spont])+((IP(x) (Subscript[γ, stim]))/(h ν)))-(Subscript[N, pump]))/(-(Subscript[γ, spont])-((2 IP(x) (Subscript[γ, stim]))/(h ν))))

FullSimplify[N+((2 (N ((Subscript[γ, spont])+((IP(x) (Subscript[γ, stim]))/(h ν)))-(Subscript[N, pump])))/(-(Subscript[γ, spont])-((2 IP(x) (Subscript[γ, stim]))/(h ν))))]

gain == ((2(Subscript[γ, stim])(((Subscript[N, pump])/( (Subscript[γ, spont])))-(N/2)))/(1+((2 IP(x) (Subscript[γ, stim]))/(h ν (Subscript[γ, spont])))))

Thus, gain is a function of optical intensity:
g== ((2(Subscript[γ, stim])(((Subscript[N, pump])/( (Subscript[γ, spont])))-(N/2)))/(1+((2 IP(x) (Subscript[γ, stim]))/(h ν (Subscript[γ, spont])))))= ((Subscript[g, 0])/(1+(( I(x) )/(Subscript[I, sat])))) 
where (Subscript[g, 0]) = 2(Subscript[γ, stim])(((Subscript[N, pump])/( (Subscript[γ, spont])))-(N/2))  and (Subscript[I, sat]) = ((Subscript[γ, spont])/(2 (Subscript[γ, stim]))) h ν 
Thus, the pump rate must be sufficiently large that the level 2 population ((Subscript[N, pump])/( (Subscript[γ, spont]))) > (N/2) for gain to occur: this is obvious. The larger the pump rate, the higher is the unsaturated gain (Subscript[g, 0]): this also makes sense. 
But now we see that as the intensity increase, it equalizes the population: there is almost as much absorption (level 1 to level 2 transitions) as stimulated radiation: this cause the net gain to drop. The saturation intensity ((Subscript[γ, spont])/(2 (Subscript[γ, stim]))) h ν is a fundamental quantity of the two-level system. In this model, it has nothing to do with pump rate. The higher the spontaneous emission rate, the more difficult it is to saturate the system, but it is also more difficult to get gain: these two aspects go together. The lower the stimulated transition probability is, the higher is the saturation intensity. This is because low (Subscript[γ, stim]) means it takes a high intensity to equalize the population, hence, a higher saturation intensity.

Thus, to have high gain, one can always pump harder. But saturation effect will always set in at the same intensity. This brings up the aspect of pump distribution. In EDFA for example, one may want higher gain when the signal grows large and lower gain when the signal is small: this brings the configuration where the pump beam enters at the other end of the fiber. The gain distribution is such that the pump is low (low gain) as the signal enters the OA, but the signal grows rapidly. As the signal gets larger, it also encounters higher pump region, and continue to see high gain.

DSolve[(1/(P[z]))(∂_z)P[z]==(((Subscript[g, 0])[z])/(1+((P[z])/Ps))),{P[z]},z]

Amplified spontaneous emission

In the above, we ignore light that comes from spontaneous emission (SE). SE is a very important feature of amplifiers. The SE can be amplified, giving rise to what is known as ASE.
P[x+Δx] = P[x] + ((Subscript[N, 2])-(Subscript[N, 1])) Δx (Subscript[γ, stim]) P[x] +β (Subscript[N, 2]) Δx A  (Subscript[γ, spont]) h ν   
where  ((Subscript[N, 2])-(Subscript[N, 1])) Δx (Subscript[γ, stim]) P[x]  is the stimulated emission, and β (Subscript[N, 2]) Δx A  (Subscript[γ, spont]) h ν is the spontaneous emission, β is the factor that indicates how much of the total (Subscript[N, 2]) Δx A (Subscript[γ, spont]) h ν that is coupled into the beam output.
((P[x+Δx]-P[x])/(Δx)) =   ((Subscript[N, 2])-(Subscript[N, 1]))  (Subscript[γ, stim]) P[x] +β (Subscript[N, 2])  A  (Subscript[γ, spont]) h ν   
d (P/(d x)) =   ((Subscript[N, 2])-(Subscript[N, 1]))  (Subscript[γ, stim]) P[x]+β (Subscript[N, 2]) A  (Subscript[γ, spont]) h ν   
Substitute from above:
d (P/(d x)) =   ((Subscript[g, 0])/(1+((P[x])/(Subscript[P, s]))))P[x] +β (Subscript[N, 2])  A (Subscript[γ, spont]) h ν  
We need an expression for (Subscript[N, 2]):
(Subscript[N, 2])=N+((N ((Subscript[γ, spont])+((IP(x) (Subscript[γ, stim]))/(h ν)))-(Subscript[N, pump]))/(-(Subscript[γ, spont])-((2 IP(x) (Subscript[γ, stim]))/(h ν)))) = N+((((Subscript[N, pump])/(Subscript[γ, spont])) -N (1+I(x) ((Subscript[γ, stim])/(h ν (Subscript[γ, spont])))))/(1+((P[x])/(Subscript[P, s]))))
     =  N+((((Subscript[N, pump])/(Subscript[γ, spont])) -N -N ((P[x])/(2 (Subscript[P, s]))))/(1+((P[x])/(Subscript[P, s])))) = N+((((Subscript[N, pump])/(Subscript[γ, spont])) -(N/2)-(N/2)(1+  ((P[x])/( (Subscript[P, s])))))/(1+((P[x])/(Subscript[P, s]))))         
    = N+(((Subscript[g, 0])/(2 (Subscript[γ, stim]) ))/(1+((P[x])/(Subscript[P, s]))))-(N/2) = (N/2)+(((Subscript[g, 0])/(2 (Subscript[γ, stim  ]) ))/(1+((P[x])/(Subscript[P, s])) )) ; since   (Subscript[g, 0]) = 2(Subscript[γ, stim])(((Subscript[N, pump])/( (Subscript[γ, spont])))-(N/2))
Which gives us:
((d P)/(d x)) =  ((Subscript[g, 0])/(1+((P[x])/(Subscript[P, s]))))P[x]  +β ((N/2)+(((Subscript[g, 0])/(2 (Subscript[γ, stim  ]) ))/(1+((P[x])/(Subscript[P, s]))))) A (Subscript[γ, spont]) h ν  
       =  ((Subscript[g, 0])/(1+((P[x])/(Subscript[P, s]))))P[x]  +β ((N/2) A (Subscript[γ, spont])h ν+((A (((Subscript[g, 0]) (Subscript[γ, spont ])h ν)/(2 (Subscript[γ, stim ]) )))/(1+((P[x])/(Subscript[P, s]))))) 
       =  ((Subscript[g, 0])/(1+((P[x])/(Subscript[P, s]))))P[x]  +β ((N/2) A (Subscript[γ, spont])h ν+(((Subscript[g, 0]) (Subscript[P, s ]))/(1+((P[x])/(Subscript[P, s])))))  ;  since: A (((Subscript[γ, spont])h ν)/(2 (Subscript[γ, stim ])))= A (Subscript[I, s]) =(Subscript[P, s]) 
       =  ((Subscript[g, 0])/(1+((P[x])/(Subscript[P, s]))))P[x]  +β ((N/2) A (Subscript[γ, spont])h ν+(((Subscript[g, 0]) (Subscript[P, s ]))/(1+((P[x])/(Subscript[P, s])))))

Summary

Key results: with ASE included, the power gain equation is:
  ((d P)/(d x)) =  ((Subscript[g, 0])/(1+((P[x])/(Subscript[P, s]))))P[x]  +β ((N/2) A (Subscript[γ, spont])h ν+(((Subscript[g, 0]) (Subscript[P, s ]))/(1+((P[x])/(Subscript[P, s])))))

The first term ((Subscript[g, 0])/(1+((P[x])/(Subscript[P, s]))))P[x] is the power gain term as we see before. The second term, 
β ((N/2) A (Subscript[γ, spont])h ν+(((Subscript[g, 0]) (Subscript[P, s ]))/(1+((P[x])/(Subscript[P, s]))))) is the ASE power term: note that it does depend on P[x]: how much ASE there is at a location depends on how much upper state population at that location, which depends on the power. Let x=0 and P[x]=0 (no input). then: 
(Subscript[(((d P)/(d x))), x=0]) =  β ((N/2) A (Subscript[γ, spont])h ν+(Subscript[g, 0]) (Subscript[P, s]) ): it is NOT zero, we can consider this as the "seeding" power of the ASE.

No saturation region

((d P)/(d x)) = P[x]  (Subscript[g, 0]) +β ((N/2) A (Subscript[γ, spont])h ν+(Subscript[g, 0]) (Subscript[P, s]) ) . Or:
 (1/(P[x]   +β (((N A (Subscript[γ, spont])h ν)/(2 (Subscript[g, 0])))+ (Subscript[P, s]) )))d (P/(d x))= (Subscript[g, 0 ])

DSolve[ (1/(P[x]   +β (((N A (Subscript[γ, spont])h ν)/(2 (Subscript[g, 0])))+ (Subscript[P, s]) ) ))(∂_x)P[x]==(Subscript[g, 0]),{P[x]},x]

P(x)->((-2 β (Subscript[g, 0]) (Subscript[P, s])+(E^((Subscript[c, 1])+x (Subscript[g, 0])))-A h N β ν (Subscript[γ, spont]))/(2 (Subscript[g, 0]))) /. x->0

FullSimplify[((-2 β (Subscript[g, 0]) (Subscript[P, s])+X (E^(x (Subscript[g, 0])))-A h N β ν (Subscript[γ, spont]))/(2 (Subscript[g, 0]))) 
			/.X->2 β (Subscript[g, 0]) (Subscript[P, s])+A h N β ν (Subscript[γ, spont])]

P(x)-> β ((Subscript[P, s])+((A N (Subscript[γ, spont]) h ν)/(2 (Subscript[g, 0]))))  ((e^(x (Subscript[g, 0])))-1)
Since (e^(x (Subscript[g, 0]))) = G: the net gain (unsaturated), we can write:
ASE power = P(x=L)= β (1+((A N (Subscript[γ, spont]) h ν)/(2 (Subscript[g, 0]) (Subscript[P, s]))))  (Subscript[P, s]) (G-1)
Notice the interesting term: ((A N (Subscript[γ, spont]) h ν)/(2 (Subscript[g, 0]) (Subscript[P, s]))), which is: 
((A N (Subscript[γ, spont]) h ν)/(2 (Subscript[g, 0]) (Subscript[P, s])))= ((A N (Subscript[γ, spont]) h ν)/( 4(Subscript[γ, stim])(((Subscript[N, pump])/( (Subscript[γ, spont])))-(N/2))  (Subscript[P, s]))) =((A N (Subscript[γ, spont]) h ν)/( 4(Subscript[γ, stim])(((Subscript[N, pump])/( (Subscript[γ, spont])))-(N/2))  (Subscript[P, s]))) = (( N (Subscript[γ, spont ]))/( 2(((2 (Subscript[γ, stim]) (Subscript[N, pump]))/(h v (Subscript[γ, spont])))-((N (Subscript[γ, stim]))/(h v)))  (Subscript[I, s]))) 
((=  N (Subscript[γ, spont ]))/( 2((( (Subscript[N, pump]))/( (Subscript[I, s])))-((N (Subscript[γ, stim]))/(h v)))  (Subscript[I, s]))) = (( N (Subscript[γ, spont ]))/( 2( (Subscript[N, pump]) - N (Subscript[γ, stim]) ((Subscript[I, s])/(h ν)))))=(( N (Subscript[γ, spont ]))/( 2( (Subscript[N, pump]) - (N/2) (Subscript[γ, spont])))) 
Thus we can rewrite:
P(L)= β (1+(( N (Subscript[γ, spont ]))/( 2( (Subscript[N, pump]) - (N/2) (Subscript[γ, spont])))))  (Subscript[P, s]) (G-1)
P(L)= β (((2 (Subscript[N, pump   ]))/( 2(Subscript[N, pump]) - N (Subscript[γ, spont])))))  (Subscript[P, s]) (G-1)
Clearly, for significant ASE, (Subscript[N, pump]) > (N/2) (Subscript[γ, spont]): the pump rate must be more than the decay rate of 1/2 of the population density. Consider the heaviest pumping condition, that's when (Subscript[N, pump]) ~ N (Subscript[γ, spont]): in this case, most of the population is inverted (in upper state), then P(x)= 2 β (Subscript[P, s]) (G-1).
When we have weak pump (very low gain), the expression becomes: 
P(x)=  β (1+((A N (Subscript[γ, spont]) h ν)/(2 (Subscript[g, 0]) (Subscript[P, s]))))  (Subscript[P, s]) (G-1) = β  (1+((A N (Subscript[γ, spont]) h ν)/(2 (Subscript[g, 0]) (Subscript[P, s])))) (Subscript[P, s])  ((e^(x (Subscript[g, 0])))-1)
					         ~ β  (1+((A N (Subscript[γ, spont]) h ν)/(2 (Subscript[g, 0]) (Subscript[P, s])))) (Subscript[P, s]) (Subscript[g, 0]) L = β  ((Subscript[g, 0])+((A N (Subscript[γ, spont]) h ν)/(2  (Subscript[P, s])))) (Subscript[P, s]) L
                            = β  (Subscript[g, 0]) (Subscript[P, s]) L+ β  ((A N (Subscript[γ, spont]) h ν)/(2  )) L 
The power is linear in L.

The expression: P(L)= β (((2 (Subscript[N, pump   ]))/( 2(Subscript[N, pump]) - N (Subscript[γ, spont])))))  (Subscript[P, s]) (G-1)=  h ν β A ((Subscript[γ, spont])/(2 (Subscript[γ, stim])))(((2 (Subscript[N, pump   ]))/( 2(Subscript[N, pump]) - N (Subscript[γ, spont]))))) (G-1)
has a very general form. The expression β A ((Subscript[γ, spont])/(2 (Subscript[γ, stim])))(((2 (Subscript[N, pump   ]))/( 2(Subscript[N, pump]) - N (Subscript[γ, spont]))))) has the unit of  1/sec

Noise from ASE

Discussion

An OA emits certain spontaneous emission power (Subscript[P, spont]). This optical power is converted to electrical signal in direct detection (DD), resulting in electrical noise.
Optical signal power: (Subscript[P, sig]) 
Net optical power out from an amplifier: (Subscript[P, sig])  G (Subscript[η, in]) (Subscript[η, out]) where (Subscript[η, in]), (Subscript[η, out]) are coupling efficiencies with the amplifier and G is the net gain.
The electrical signal power of the amplified signal is:  
S=((((Subscript[P, sig]) G (Subscript[η, in]) (Subscript[η, out])L (Subscript[R, det])))^2) =  (((G (Subscript[η, in]) (Subscript[η, out]) L (Subscript[I, sig])))^2) where (Subscript[R, det ]) is the detector responsivity, L is the optical loss between amplifier and detector.
The shot noise current of this signal is:
2 e G (Subscript[η, in]) (Subscript[η, out]) L (Subscript[I, sig]) (Subscript[BW, E]) where (Subscript[BW, E]) is the electrical bandwidth.
From the amplifier, there is ASE which also has a current:
(Subscript[I, ASE]) =(Subscript[P, spont]) L (Subscript[η, out])(Subscript[R, det])
This current has a shot noise:
2 e (Subscript[I, ASE])  (Subscript[BW, E]).
Notice that if the input coupling (Subscript[η, in]) is poor, the ASE shot noise is still at a high level while the amplified signal is low: this way, OA is less beneficial. This is a way to show the criticality of efficient coupling of initial signal into the OA.

Beat noise

Unfortunately, that is not all the noise from the OA. The (Subscript[P, spont])[t]  is not a constant-average function with a random fluctuation distributed with photon Poisson statistics. It has noises that come from the beat terms. To treat this, let's consider the optical electric field:
(Subscript[E, out])[t] = (Subscript[A, sig]) (e^(i (Subscript[ω, 0]) t )) and
(Subscript[E, spont])[t] = (∫_(-(Subscript[BW, O])/2)^((Subscript[BW, O])/2))A[ν] (e^(i (((Subscript[ω, 0]) +2 π ν) t + Φ[ν]))) dν 
where (Subscript[ω, 0]) is the central frequency component, A[ν]  and  Φ[ν] are the amplitude and phase, which is random for the latter. A minor problem with the representation of (Subscript[E, spont])[t] = (∫_(-(Subscript[BW, O])/2)^((Subscript[BW, O])/2))A[ν] (e^(i (((Subscript[ω, 0]) +2 π ν) t + Φ[ν]))) dν as an integral is that it does not represent the signal as a power signal, but more like an energy signal. So, if we let A[ν]= (∑_(m=-N)^N)(Subscript[A, m]) δ[ν-(Subscript[ν, m])]

then:
(Subscript[E, spont])[t] = (∑_(m=-N)^N)(Subscript[A, m]) (e^(i (((Subscript[ω, 0]) +2 π (Subscript[ν, m])) t + Φ[(Subscript[ν, m])]) )) with the understanding that the phase  Φ[(Subscript[ν, m])] is random, and N is an arbitrary number that does to ∞ in the limit.
These electrical fields can beat with each other. This gives rise to signal-spontaneous beat and spontaneous-spontaneous beat.

Reviewing old notes (see noise)

Beat terms. If we have:
(Subscript[F, 1]) (e^(i (Subscript[Φ, 1]))) +(Subscript[F, 2]) (e^(i (Subscript[Φ, 2]) )), where (Subscript[F, 1 ]) and (Subscript[F, 2 ]) are pure amplitude (real, no phase);
The intensity of the mix is: 
 ((Re[(Subscript[F, 1]) (e^(i (Subscript[Φ, 1]))) +(Subscript[F, 2]) (e^(i (Subscript[Φ, 2])))])^2)  =    ((((Subscript[F, 1]) Cos[(Subscript[Φ, 1])]+(Subscript[F, 2]) Cos[(Subscript[Φ, 2])]))^2)

Expand[((((Subscript[F, 1]) Cos[(Subscript[Φ, 1])]+(Subscript[F, 2]) Cos[(Subscript[Φ, 2])]))^2)]

TrigReduce[%]

((Re[(Subscript[F, 1]) (e^(i (Subscript[Φ, 1]))) +(Subscript[F, 2]) (e^(i (Subscript[Φ, 2])))])^2) =(1/2) (cos(2 (Subscript[Φ, 1])) (Subscript[F, 1]^2)+(Subscript[F, 1]^2)+2 cos((Subscript[Φ, 1])-(Subscript[Φ, 2])) (Subscript[F, 2]) (Subscript[F, 1])
					                                      +2 cos((Subscript[Φ, 1])+(Subscript[Φ, 2])) (Subscript[F, 2]) (Subscript[F, 1])+cos(2 (Subscript[Φ, 2])) (Subscript[F, 2]^2)+(Subscript[F, 2]^2)) 
For square-law detector, time-averaged over optical frequency results in all high-frequency term vanish: cos(2 (Subscript[Φ, 1])), cos((Subscript[Φ, 1])+(Subscript[Φ, 2])), and cos(2 (Subscript[Φ, 2])); we are left with DC and low frequency terms:
(1/2) ((Subscript[F, 1]^2)+(Subscript[F, 2]^2)+2 cos((Subscript[Φ, 1])-(Subscript[Φ, 2])) (Subscript[F, 2]) (Subscript[F, 1]))

DC components:  (1/2) ((Subscript[F, 1]^2)+(Subscript[F, 2]^2)) 
Beat component: (Subscript[F, 2]) (Subscript[F, 1])cos((Subscript[Φ, 1])-(Subscript[Φ, 2])) 
Modulation depth: ((2 (Subscript[F, 2]) (Subscript[F, 1]))/((Subscript[F, 1]^2)+(Subscript[F, 2]^2)))

Useful relation

We can also use the useful relation of time-averaging of complex harmonics. Recalling this relation:  let S be any harmonic wave quantity,  
((Re[S])^2) = (1/4) (S+(SuperStar[S]))(S+(SuperStar[S])) = (1/4)((S^2) + 2 (SuperStar[S])S +((SuperStar[S])^2)). Let S be the sum of a number of waves. Then (S^2) and ((SuperStar[S])^2) consist of high frequency terms, which vanish with time average. Thus, we have the relation:
(Subscript[⟨((Re[S])^2)⟩, time average]) =  (1/2) (SuperStar[S])S =(1/2)(((|S|))^2)  
which is sometime more convenient to work with.
(Subscript[⟨((Re[(Subscript[A, sig]) (e^(i (Subscript[ω, 0]) t ))])^2)⟩, time average]) =  (1/2) (SuperStar[(Subscript[A, sig])])(Subscript[A, sig]) =(1/2)(((|(Subscript[A, sig])|))^2) = signal optical power 
So, for unit: (Subscript[A, sig]): (Sqrt[2 × power ])

So now we have: ((Subscript[A, sig]) (e^(i (Subscript[ω, 0]) t ))+(∑_(n=-N)^N)(Subscript[A, spont])[(Subscript[ν, n])] (e^(i (((Subscript[ω, 0]) +2 π (Subscript[ν, n])) t + Φ[(Subscript[ν, n])]) )), such that (Subscript[ν, n]) spands the entire optical bandwidth (Subscript[BW, O]), and: 
(∑_(n=-N)^N)(((|(Subscript[A, spont])[(Subscript[ν, n])]|))^2) = 2 × total spontaneous power within (Subscript[BW, O]).
We can immediately drop high frequency terms, and left with 
((Subscript[A, sig])^2) + (∑_(n=-N)^N)(∑_(m=-N)^N)(SuperStar[(Subscript[A, spont])])[(Subscript[ν, m])] (Subscript[A, spont])[(Subscript[ν, n])](e^(i (2 π ((Subscript[ν, n])-(Subscript[ν, m])) t + Φ[(Subscript[ν, n])]-Φ[(Subscript[ν, m])])))  
+  (Subscript[A, sig])(∑_(n=-N)^N)(Subscript[A, spont])[(Subscript[ν, n])] (e^(i (((Subscript[ω, 0]) +2 π (Subscript[ν, n])) t + Φ[(Subscript[ν, n])]))) (e^(-i (Subscript[ω, 0]) t ))
+  (Subscript[A, sig])(∑_(m=-N)^N)(SuperStar[(Subscript[A, spont])])[(Subscript[ν, m])]  (e^(-i (((Subscript[ω, 0]) +2 π (Subscript[ν, m])) t + Φ[(Subscript[ν, m])])))(e^(i (Subscript[ω, 0]) t )) dν
 
Or:
((Subscript[A, sig])^2) + (∑_(n=-N)^N)(∑_(m=-N)^N)(SuperStar[(Subscript[A, spont])])[(Subscript[ν, m])] (Subscript[A, spont])[(Subscript[ν, n])](e^(i (2 π ((Subscript[ν, n])-(Subscript[ν, m])) t + Φ[(Subscript[ν, n])]-Φ[(Subscript[ν, m])]))) 
+  (Subscript[A, sig])(∑_(n=-N)^N)(Subscript[A, spont])[(Subscript[ν, n])] (e^(i  (2 π (Subscript[ν, n]) t + Φ[(Subscript[ν, n])]))) 
+  (Subscript[A, sig])(∑_(m=-N)^N)(SuperStar[(Subscript[A, spont])])[(Subscript[ν, m])]  (e^(-i (2 π (Subscript[ν, m]) t + Φ[(Subscript[ν, m])])))
Out of 4 terms, only ((Subscript[A, sig])^(2 ))is the signal, all others have random phase Φ[ν] which are noises. But the diagonal of (∑_(n=-N)^N)(∑_(m=-N)^N)(SuperStar[(Subscript[A, spont])])[(Subscript[ν, m])] (Subscript[A, spont])[(Subscript[ν, n])](e^(i (2 π ((Subscript[ν, n])-(Subscript[ν, m])) t + Φ[(Subscript[ν, n])]-Φ[(Subscript[ν, m])]))) for (Subscript[ν, n])=(Subscript[ν, m]) is the total spontaneous power, which is NOT the beat term beween the spontaneous power. 
It is: (1/2)(∑_(n=-N)^N) (((Subscript[A, spont])[(Subscript[ν, n])])^2) = (Subscript[P, spont])[(Subscript[ω, 0])] (Subscript[BW, 0]), where (Subscript[P, spont])[(Subscript[ω, 0])] is the power spectral density of the spont. emission in unit of W/Hz. (Note that this will be different from the optical noise PSD, which is the noise of optical power which is in unit of (W^2)/Hz).

Signal--spontaneous beat noise:

Signal-spontaneous beat noise is
  (Subscript[A, sig])(∑_(n=-N)^N)(Subscript[A, spont])[(Subscript[ν, n])] ((e^(i  (2 π (Subscript[ν, n]) t + Φ[(Subscript[ν, n])])))+(e^(-i  (2 π (Subscript[ν, n]) t + Φ[(Subscript[ν, n])])))) 
 The power spectral density of the sig-spont noise is:  
((Subscript[A, sig])^2)(∑_(n=-N)^N)(((Subscript[A, spont])[(Subscript[ν, n])])^2) (δ[f-(Subscript[ν, n])]+δ[f+(Subscript[ν, n])])= 2 ((Subscript[A, sig])^2)(Subscript[P, spont])[(Subscript[ω, 0])] where (Subscript[P, spont])[(Subscript[ω, 0])] is described as above. (note the unit: it is W·W/Hz =(W^2)/Hz).
What about the noise distribution?
We write the beat term again:
 (Subscript[A, sig])(∑_(n=-N)^N)(Subscript[A, spont])[(Subscript[ν, n])] ((e^(i  (2 π (Subscript[ν, n]) t + Φ[(Subscript[ν, n])])))+(e^(-i  (2 π (Subscript[ν, n]) t + Φ[(Subscript[ν, n])]))))
           = 2 (Subscript[A, sig])(∑_(n=-N)^N)(Subscript[A, spont])[(Subscript[ν, n])]Cos[2 π (Subscript[ν, n]) t + Φ[(Subscript[ν, n])]]
 If Φ[(Subscript[ν, m])] is a white noise random function, the central limit theorem dictate that the summation at very large N approaches that of a Gaussian. Thus, this is why we can use Gaussian noise model.

Spontaneous-spontaneous beat noise

Spontaneous-spontaneous term noise is
 (∑_(n=-N)^N)(∑_(m=-N)^N)(SuperStar[(Subscript[A, spont])])[(Subscript[ν, m])] (Subscript[A, spont])[(Subscript[ν, n])](e^(i (2 π ((Subscript[ν, n])-(Subscript[ν, m])) t + Φ[(Subscript[ν, n])]-Φ[(Subscript[ν, m])]))) 
We'll assume the optical bandwidth is divided evenly, so that: 
(Subscript[ν, n])=n Δν; Δν=((Subscript[BW, O])/(2N)). Then:
(∑_(n=-N)^N)(∑_(m=-N)^N)(SuperStar[(Subscript[A, spont])])[(Subscript[ν, m])] (Subscript[A, spont])[(Subscript[ν, n])](e^(i 2 π(n-m) Δν t ))(e^(i  ((Subscript[Φ, n])-(Subscript[Φ, m])))) 
We can rewrite, substituting  n-m = p  or  n = m + p  
(∑_(p=-2N)^(2N))(∑_(m=Max[-N,-N-p])^(Min[N,N-p]))(Subscript[A, spont])[(Subscript[ν, m])] (Subscript[A, spont])[ (Subscript[Δν, m+p])](e^(i 2 π p Δν t ))(e^(i  ((Subscript[Φ, m+p])-(Subscript[Φ, m])))) 
For  given p, m=n-p
p= -2 N; m=N, n= -N  => m=-p-N to -p-N  -  (2N+p+1)=-p-N-2N-p-1=-2p -3N -1
p=-2N+1, m=N-1, n=-N; or m=N, n=-N+1
p=0; m=-N to +N: 2 N term
p=2N-1; m=-N+1,n=N or m=-N, n=N-1
p=2N; m=-N; n=N
The PSD is:
(∑_(m=-N)^N)(∑_(p=Max[-N,-N-m])^(Min[N,N-m]))(((Subscript[A, spont])[(Subscript[ν, m])])^2) (((Subscript[A, spont])[ (Subscript[Δν, m+p])])^2)δ[f-(Subscript[ν, p])] 
It is thus a triangular function with  (((Subscript[P, spont])[(Subscript[ω, 0])])^2) (Subscript[BW, O]) (1-(f/(Subscript[BW, O])))

Use of W-K theorem

Autocorrelation signal: R[τ]=(∫_(-∞)^∞)x[t](SuperStar[x])[t+τ]dt
Let: x[t]=(∫_(-∞)^∞)X[f](e^(-i 2 π f t)) df, then:
                              R[τ]=(∫_(-∞)^∞)(∫_(-∞)^∞)X[f'](e^(-i 2 π f' t)) df' (∫_(-∞)^∞)(SuperStar[X])[f](e^(i 2 π f (t+τ))) dfdt
                              R[τ]=(∫_(-∞)^∞)(∫_(-∞)^∞)X[f'](SuperStar[X])[f](e^(i 2 π f τ)) ((∫_(-∞)^∞)(e^(i 2 π (f'-f)t)) dt )df df'
                              R[τ]=(∫_(-∞)^∞)(∫_(-∞)^∞)X[f'](SuperStar[X])[f](e^(i 2 π f τ)) δ[f'-f]df df'  
                              R[τ]=(∫_(-∞)^∞)|X[f]((|)^2) (e^(i 2 π f τ)) df  
Or, the autocorrelation function is the Fourier transform of |X[f]((|)^2), which is the power/energy spectral density of the signal x[t]. This is Wiener-Kintchin theorem. We can also write:
                              |X[f]((|)^2)=(∫_(-∞)^∞) R[τ] (e^(-i 2 π f τ)) df

Let x[t]=(e^(-i 2 π f t)), what is its autocorrelation function?
         R[τ]=(∫_(-∞)^∞)x[t](SuperStar[x])[t+τ]dt =(∫_(-∞)^∞)(e^(-i 2 π f t))(e^(i 2 π f (t+τ)))dt
                               =(∫_(-∞)^∞)(e^(i 2 π f τ))dt = (e^(i 2 π f τ)) (∫_(-∞)^∞)1dt
Now, we will let: x[t]=(e^(-i 2 π f t-α t)) u[t],
         R[τ]=(∫_(-∞)^∞)x[t](SuperStar[x])[t+τ]dt =(e^(i 2 π f τ-α τ ))(∫_(-∞)^∞)(e^(- 2 α t))u[t] u[t+τ] dt
                =(e^(i 2 π f τ-α τ ))(∫_(Max[0,-τ] )^∞)(e^(- 2 α t))dt

(∫_(σ )^∞)(E^(- 2 α t))\[DifferentialD]t

Thus:                R[τ]=(e^(i 2 π f τ-α τ ))((e^(-2 α  Max[0,-τ]))/(2 α))
If τ>0:               R[τ]=(e^(i 2 π f τ-α τ ))(1/(2 α))
if τ<0:               R[τ]=(e^(i 2 π f τ-α τ ))((e^(2 α τ))/(2 α)) = (e^(i 2 π f τ+α τ ))(1/(2 α))   
Combine both:  R[τ]=(1/(2 α))(e^(i 2 π f τ-α |τ|))

(∫_0^∞)(1/(2 α))(E^(-α τ))(E^(I 2 π (f-g) τ)) \[DifferentialD]τ

(∫_(-∞)^0)(1/(2 α))(E^(α τ))(E^(I 2 π (f-g) τ)) \[DifferentialD]τ

FullSimplify[(1/(2 (α^2)-4 I f π α+4 I g π α))+(1/(2 (α^2)+4 I f π α-4 I g π α))]

(1/(4 (π^2) (((f-g))^2)+(α^2))) /. f->(Subscript[f, 0]) /. g-> f

We see that S[f]=(1/(4 (π^2) (((f-(Subscript[f, 0])))^2)+(α^2))) is the energy spectral density of the signal.

Cross term

Example: Constant A[ν]:

((Subscript[A, spont])^4) (∑_(m=-N)^N)(∑_(p=Max[-N,-N-m])^(Min[N,N-m]))δ[f-(Subscript[ν, p])] 
For a given p the number of allow m and n is a decreasing triangular.
So: ((Subscript[A, spont])^(4 ))Abs[1-(( f)/BW)]
At f=0; the sum is: ((Subscript[A, spont])^4) * 2N = total pont power; Thus, the PSD at f=0 = (((Subscript[P, spont])^2)/(Subscript[BW, 0]))

N[0]=((((Subscript[A, 0])^4)(|(∫_(-(Subscript[BW, O])/2)^( (Subscript[BW, O])/2-0))  dν|))^2) =((Subscript[A, 0])^4) ((Subscript[BW, O])^2): this is the maximum spont-spont noise power spectral density. How does it vary vs. f? 
We write:
 (((|(∫_(-(Subscript[BW, O])/2)^( (Subscript[BW, O])/2-f))  (e^(i ( Φ[ν]-Φ[ν + f]))) dν|))^2)=
	                     (∫_(-(Subscript[BW, O])/2)^( (Subscript[BW, O])/2-f)) (∫_(-(Subscript[BW, O])/2)^( (Subscript[BW, O])/2-f))  (e^(i ( Φ[ν]-Φ[ν + f]))) (e^(-i ( Φ[ν']-Φ[ν' + f])))dν'dν
= (∫_(-(Subscript[BW, O])/2)^( (Subscript[BW, O])/2-f)) (∫_(-(Subscript[BW, O])/2)^( (Subscript[BW, O])/2-f))  (e^(i ( Φ[ν]-Φ[ν']))) (e^(-i ( Φ[ν + f]-Φ[ν' + f])))dν'dν
Since Φ[ν] is random function, any uncorrelation between Φ[ν]-Φ[ν'] would indicate that the integration must be zero for ν≠ν'. However, we must think in terms of discrete distribution of ν: Every photon must have a finite bandwidth and its phase is fixed within a certain decay time. Thus, converting the integration into discrete sum of N frequencies: 
 Power= ((Subscript[A, 0])^2)(Subscript[BW, O]); Δν = ((Subscript[BW, O])/(N  ))
(∫_(-(Subscript[BW, O])/2)^( (Subscript[BW, O])/2-f))  (e^(i ( Φ[ν]-Φ[ν']))) (e^(-i ( Φ[ν + f]-Φ[ν' + f]))) dν' 
		                        →(a^2) Δν(∑_(m=1)^(P=(((Subscript[BW, O])-f)/(Δν))))(e^(i ( Φ[(Subscript[ν, k])]-Φ[(Subscript[ν, m])']))) (e^(-i ( Φ[(Subscript[ν, k+r]) ]-Φ[(Subscript[ν, m+r])'])))
                                   =(a^2) Δν  
Thus:
((Subscript[A, 0])^4)(∫_(-(Subscript[BW, O])/2)^( (Subscript[BW, O])/2-f)) (∫_(-(Subscript[BW, O])/2)^( (Subscript[BW, O])/2-f))  (e^(i ( Φ[ν]-Φ[ν']))) (e^(-i ( Φ[ν + f]-Φ[ν' + f])))dν'dν
                              →((Subscript[A, 0])^4) ((Δν)^2)(∑_(k=1)^(P=(((Subscript[BW, O])-f)/(Δν))))1 = ((Subscript[A, 0])^4) ((Δν)^2) (((Subscript[BW, O])-f)/(Δν))=((Subscript[A, 0])^4) ((BW((Subscript[BW, O])-f))/N)
N[f]=((Subscript[A, 0])^4) (Subscript[BW, O]) ((Subscript[BW, O])-f)