optimization of output coupling

define units, constants

m=1;cm=0.01;mm=0.001;um=(10^(-6));nm=(10^(-9));
J=1;mJ=(10^(-3));uJ=(10^(-6));
s=1;ps=(10^(-12));ns=(10^(-9));us=(10^(-6));ms=(10^(-3));
W=1;MW=(10^6)W;
eV=1.6 (10^(-19))J;
c=3 (10^8) m/s;
hbar=(10^(-34))J s;

define relations

the cavity modeling is very simple, assuming 4-level system, no space-dependent pumping and mode effects
aside from the ratio of λlaser/λpump, there are no pump efficiencies accounted for. 

define cavity loss parameters

γp[Γp_]:=-Log[1-Γp]; (* total passive loss, not OC *)
γoc[Γoc_]:=-Log[1-Γoc];(* passive OC loss = (Subscript[γ, 2]) in Svelto *)
γ[Γp_,Γoc_]:=γp[Γp]+γoc[Γoc]/2;(* total single-pass, not RT, passive cavity loss *)

cross-section, saturation intensity
mode area

σ:=ephot/Fsat;(* stimulated emission coeff, in terms of sat. fluence *)
Isat:=(ephot/(σ τ));
modeArea:=π (modeD^2)/4;(* area of resonator mode *)

define threshold power and output power
(Subscript[P, th])=γ ((A(Subscript[, b]))/σ)((h ν)/τ)
(Subscript[P, out])=(Subscript[A, b])(Subscript[I, s])((Subscript[γ, 2])/2)(((Subscript[P, in])/((Subscript[P, th])[γ]))-1)

Pth[Γp_,Γoc_]:=γ[Γp,Γoc]((modeArea ephot)/(σ τ)); (* threshold power *)
Pout[Pin_,Γp_,Γoc_]:=Max[modeArea Isat ((γoc[Γoc])/2)((Pin/(Pth[Γp,Γoc]))-1),0];

λphot=633nm;
ephot=hbar ((2 π c)/(λphot));
Fsat=0.1 mJ/(cm^2);
modeD=1mm;
τ=100ns;

Γp=0.05;
ΓocTst=0.03;
Print["threshold power: ",Pth[Γp,ΓocTst]];
Plot[Pout[Pin,Γp,ΓocTst],{Pin,0,3},GridLines->{{Pth[Γp,ΓocTst]},None}]

vary output coupling, see how output power changes.
note how threshold power goes up with more OC

Manipulate[Plot[Pout[Pin,Γp,Γoc],{Pin,0,3},PlotRange->{0,1.5},GridLines->{{Pth[Γp,Γoc]},None}],{Γoc,0.001,0.4}]

this shows the variation of the threshold power with OC
the y-intercept is non-zero because of the fixed internal passive losses

Plot[Pth[Γp,Γoc],{Γoc,0,0.4},AxesOrigin->{0,0}]

now vary the output coupling to see where the output power maximizes for fixed input power

PinTst=2;
Plot[Pout[PinTst,Γp,Γoc],{Γoc,0,0.4}]

other relations (not used)

Nth:=γ/(σ lenCry); (* threshold inversion density *)
RpCrit:=Nth/τ; (* critical inversion density *)
ϕ0:=modeArea (γ/σ)((τc)/τ)((Rp/RpCrit)-1); (* CW number of circulating internal photons*)
lenEff:=lenCav +(nCry-1)lenCry;(* effective length of cavity *)
τc:=lenEff/(γ c);(* photon lifetime in cavity *)
τRT:=2 lenEff/c;(* cavity round trip time *)
Va:=modeArea lenCry;(* volume of active medium *)
Vcav:=modeArea lenEff;(* volume of cavity (assumes same mode size throughout *)
Estor[Ninv_]:=Ninv ephot Va;(* stored energy from inversion density *)
Fstor[Ninv_]:=Ninv lenCry ephot;(* stored fluence *)
G0[Ninv_]:=Exp[Fstor[Ninv]/Fsat]; (* small-signal single-pass gain *)
Pabs[Rp_]:=Rp ephotPump Va;(* absorbed pump power required to maintain Rp *)
PoutFromϕ[ϕ_]:=((γoc c)/(2 lenEff))ephot ϕ;(* output power of laser, from ϕ *)
Pint[ϕ_]:=Pout[ϕ]/Γoc; (* internal circulating power based on output power, OC *)
B:=((σ c)/Vcav); (* effective Einstein B value *)