Pound-Drever-Hall Error Signal Derivation

Derivation of PDH error signal, following March 18, 2026 tutorial here:
https://ccahilla.github.io/lasers-and-optomechanics/tutorials/

0. Fabry Perot cavity reflection as a function of phase

refl[ϕ_]:=((r1-r2 Exp[-I 2ϕ])/(1-r1 r2 Exp[-I 2ϕ]));

1. What is the total electric field incident on the front of the cavity?

Carrier input

Ein=E0 Exp[I ω0 t];

Carrier + Upper Sideband + Lower Sideband

Eintotal=Ein(1+I (Γ/2)Exp[I Ω t]+I (Γ/2)Exp[-I Ω t])

2. Write down three single-pass phases ϕ0 , ϕusb , ϕlsb that the carrier, upper sideband, and

lower sideband, will experience as they traverse the cavity, in terms of the cavity length L,

and carrier and RF frequencies ω0 and Ω.

Where L is the nominal cavity length:

params={ϕ0->((ω0 L)/c),ϕusb->(((ω0+Ω)L)/c),ϕlsb->(((ω0-Ω)L)/c)};

I have make this a list of parameters, so we can leave things in terms of ϕ0,ϕusb,ϕlsb,
Notice that we can also do another phase breakdown where 

ϕusb = ϕ0 + ϕrf and ϕlsb = ϕ0 - ϕrf
where
ϕrf = ((Ω L)/c).

3. Write down the cavity total reflected field Erefl in terms of the reflection function r(ϕ) and

ϕ0 , ϕusb , ϕlsb

r0 is the reflected carrier refl[ϕ0],
rpΩ is the reflected upper sideband refl[ϕusb],
rmΩ is the reflected lower sideband refl[ϕlsb].

reflparams={r0->refl[ϕ0],rpΩ->refl[ϕusb],rmΩ->refl[ϕlsb]};

Erefltotal=Ein(r0+I (Γ/2)rpΩ Exp[I Ω t]+I rmΩ(Γ/2)Exp[-I Ω t]);

4. Write down the cavity total reflected power (Subscript[P, refl]) (t) in terms of the reflection function r(ϕ)

and ϕ0 , ϕusb , ϕlsb .

Use ComplexExpand, which assumes all variables are real variables, with Conjugate to quickly calculate our power.
However, in Mathematica we need to be careful because r0, rusb, and rlsb are NOT real variables, they are complex.
So when we take the complex conjugate we need to keep track of them.

I keep track of this manually below, by replacing r0, rusb, and rlsb with r0c, rusbc, and rlsbc, where the “c” suffix indicates complex conjugate.  
Not elegant but it works.

TrigToExp expands our trig expressions we get from FullSimplify into exponentials, which is easier to reach given our starting point

reflcparams={r0c->refl[-ϕ0],rpΩc->refl[-ϕusb],rmΩc->refl[-ϕlsb]};

Erefltotalc=TrigToExp[ComplexExpand[Conjugate[Erefltotal]]]/.{r0->r0c,rpΩ->rpΩc,rmΩ->rmΩc}

Prefltotal=TrigToExp[FullSimplify[ComplexExpand[Erefltotal Erefltotalc]]]

Prefltotal=Prefltotal/.{(Γ^2)->0}

5. Set the coefficient for (e^(-i Ω t)) equal to some complex number z, and the coefficient for (e^(i Ω t)) equal to y.  What is z and y?  
Is there any connection between z and y? 
What is it?

First, Collect our (Subscript[P, refl])(t) by using (e^(-i Ω t)) and (e^(i Ω t)) coefficients.

Collect[Prefltotal,{(E^(I t Ω)),(E^(-I t Ω))}]

The first term above is the DC term, the second term associated with (e^(i t Ω)) is the y term, and the third term associated with (e^(-i t Ω)) is the z term:

z=Simplify[((1/2) I (E0^2) r0c rmΩ Γ-(1/2) I (E0^2) r0 rpΩc Γ)]
y=Simplify[(-(1/2) I (E0^2) r0 rmΩc Γ+(1/2) I (E0^2) r0c rpΩ Γ)]

Analyzing z and y we see that y is the complex conjugate of z.  
We can write y = (SuperStar[z]).

6. Rewrite (Subscript[P, refl]) in terms of z and y.

Prefl=(E0^2) r0 r0c+(SuperStar[z]) (e^(i t Ω)) +z (e^(-i t Ω))

7. (Subscript[P, I]) + i (Subscript[P, Q]) = (P^Ω) . Show we can rewrite the demodulation in this way.

Cos[Ωt]+i Sin[Ωt] = (e^(i Ω t)), so we can combine our two real demodulations into a single complex demodulation.  This will gives us a complex number whose real and imaginary parts represent our usual I and Q outputs.

8. Calculate the 1Ω demodulated power ((Subscript[P, refl])^Ω)(ϕ0, ϕusb , ϕlsb). 
What can we see pops out?

Looking at our Prefl expression above:
Prefl=(E0^2) r0 r0c+(SuperStar[z]) (e^(i t Ω)) +z (e^(-i t Ω)) 
it’s clear that if we multiply by (e^(i Ω t)) and integrate Prefl by one cycle in Ωt, the DC term and (e^(i t Ω)) will go to (e^(i t Ω)) and (e^(2i t Ω)), and will integrate to zero over one cycle.

BUT the (e^(-i t Ω)) term goes to stationary in time, so we should get exactly z out of our demodulation:
((Subscript[P, refl])^Ω) = z.

PreflΩ=z

9. Calculate the derivative of ((Subscript[P, refl])^Ω) with respect to ϕ0, then evaluate it about zero.

First, we must substitute all instances of ϕ0: ϕusb -> ϕ0+ϕrf, ϕlsb -> ϕ0 - ϕrf, as well as our reflection terms r(ϕ0).
Then we take the derivative, 
then evaluate at ϕ0 -> 0.

dPreflΩdϕ0=Simplify[D[PreflΩ /. reflparams /. reflcparams /. {ϕusb->ϕ0+ϕrf,ϕlsb->ϕ0-ϕrf},ϕ0]];

dPreflΩdϕ0at0=Simplify[dPreflΩdϕ0/.{ϕ0->0}]

10. Plot the real and imaginary parts of ((Subscript[P, refl])^Ω)((Subscript[ϕ, 0]), ϕusb, ϕlsb ) as a function of ϕ0 ∈ [−π, π].

Below we evaluate our (Subscript[(P^Ω), refl]) signal.  
The main signal appears in the real quadrature, 
with some small blip in Q-quadrature due to a slight overlap in carrier and RF sideband resonances.
The discriminant (Subscript[(dP^Ω), refl])/(Subscript[dϕ, 0])×(Subscript[ϕ, 0]) is plotted as a dashed line, and gives our optical gain about resonance.

This error signal has a very distinctive shape, and is highly non-linear outside the region of carrier resonance.
But inside the resonance, it is very discriminating, and able to provide an excellent, strong, linear error signal.
This can be detected with a radio-frequency photodetector (RFPD). 

The RF sidebands do show up in our PDH error signal sweep, as the side lobes on either side of the carrier resonance.
The side lobes represent when the Upper Sideband or Lower Sideband starts to resonant inside the cavity.
The RF sidebands have a strong response in the imaginary quadrature, 
while the carrier does not.

values={L->1,T1->0.1,T2->0.1,r1->(Sqrt[1-T1]),r2->(Sqrt[1-T2]),
E0->1,Γ->0.1,ϕusb->ϕ0+ϕrf,ϕlsb->ϕ0-ϕrf,ϕrf->(π/4)};

PreflΩ/.reflparams/.reflcparams;

plotPreflΩ=PreflΩ/.reflparams/.reflcparams//.values;

plotdPreflΩdϕ0at0=dPreflΩdϕ0at0//.values;

Plot[{Re[plotPreflΩ],Re[plotdPreflΩdϕ0at0 ϕ0]},{ϕ0,-(π/2),(π/2)},PlotRange->{-0.06,0.06},PlotStyle->{Thickness[0.005],Dashed},PlotLegends->{"PDH Error Signal (Subscript[P, refl])((Subscript[ϕ, 0]))", "PDH discriminant ((Subscript[dP, refl])/(Subscript[dϕ, 0]))\[Cross](Subscript[ϕ, 0])"},AxesLabel->{"Carrier\nPhase (Subscript[ϕ, 0])\n[rads]","Re[(Subscript[P, refl])]"},ImageSize->Large]
Plot[{Im[plotPreflΩ],Im[plotdPreflΩdϕ0at0 ϕ0]},{ϕ0,-(π/2),(π/2)},PlotRange->{-0.06,0.06},PlotStyle->{Thickness[0.005],Dashed},AxesLabel->{"Carrier\nPhase (Subscript[ϕ, 0])\n[rads]","Im[(Subscript[P, refl])]"},ImageSize->Large]

Manipulate plot of PDH error signal ((Subscript[P, refl])/Γ)(ϕ0,ϕrf)

For fun, in the manipulate plot, I’ve added in the Demod Angle ζ, which is just a simple rotation of the detected I and Q quadrature.
The way to consider the demod phase is similar to the arbitrary phase accrued while transmitting to the photodetector.
We can explicitly consider it while demodulating by demodulating over Cos[Ωt+ζ] and Sin[Ωt+ζ].  
If we set ζ = (π/2), we’ll flip which quadrature our signal is in.

manivalues2={L->1,r1->(Sqrt[1-T1]),r2->(Sqrt[1-T2]),E0->1,ϕusb->ϕ0+ϕrf,ϕlsb->ϕ0-ϕrf}

maniPreflΩ[ϕ0_,ϕrf_,Γ_,T1_,T2_]:=Evaluate[PreflΩ/.reflparams/.reflcparams//.manivalues2]

Manipulate[
Column[{
Plot[Re[maniPreflΩ[ϕ0,ϕrf,Γ,T1,T2]Exp[I ζ]],{ϕ0,-(π/2),(π/2)},PlotRange->Full,PlotLegends->{"PDH Error Signal (Subscript[P, refl])((Subscript[ϕ, 0]))", "PDH discriminant ((Subscript[dP, refl])/(Subscript[dϕ, 0]))\[Cross](Subscript[ϕ, 0])"},AxesLabel->{"Carrier\nPhase (Subscript[ϕ, 0])\n[rads]","Re[(Subscript[(P^Ω), refl])]"},ImageSize->Large,PlotPoints->200],
Plot[Im[maniPreflΩ[ϕ0,ϕrf,Γ,T1,T2]Exp[I ζ]],{ϕ0,-(π/2),(π/2)},PlotRange->Full,AxesLabel->{"Carrier\nPhase (Subscript[ϕ, 0])\n[rads]","Im[(Subscript[(P^Ω), refl])]"},ImageSize->Large,PlotPoints->200]
}],
{{ϕrf,(π/4),"RF Phase ϕrf [rads]"},0,(π/2)},
{{Γ,0.1,"RF Modulation Depth Γ [rads]"},0,0.5},
{{T1,0.1,"Input Mirror Transmission (Subscript[T, 1])"},0,1},
{{T2,0.1,"End Mirror Transmission (Subscript[T, 2])"},0,1},
{{ζ,0,"Demod Angle ζ [rads]"},-π,π},
Paneled->False
]

11. Plot ((Subscript[P, refl])^Ω) as a phasor as a function of φ0 ∈ [− (π/2),(π/2)]

Below we evaluate our (Subscript[(P^Ω), refl]) signal as a phasor.
The carrier resonance is the line almost perfectly along the real axis.
The circles are the RF sideband sidelobes.

Those circles are representative of the reflection of the RF sidebands off the cavity, rotated by either + or -90°.  
If we increase the cavity finesse by lowering the mirror transmission T2 or T1, our signals get stronger and cleaner,
as the improved finesse of the cavity shrinks the size of the resonance, preventing mixing.

ParametricPlot[{Re[plotPreflΩ],Im[plotPreflΩ]},{ϕ0,-(π/2),(π/2)},AspectRatio->1,PlotRange->Full,AxesLabel->{"Re[(Subscript[(P^Ω), refl])]","Im[(Subscript[(P^Ω), refl])]"},ImageSize->Large]

Manipulate plot of PDH phasor

Manipulate[
ParametricPlot[{Re[maniPreflΩ[ϕ0,ϕrf,Γ,T1,T2]Exp[I ζ]],Im[maniPreflΩ[ϕ0,ϕrf,Γ,T1,T2]Exp[I ζ]]},{ϕ0,-(π/2),(π/2)},PlotRange->Full,PlotLegends->{"PDH Error Signal (Subscript[P, refl])((Subscript[ϕ, 0]))"},AxesLabel->{"Re[(Subscript[(P^Ω), refl])]","Im[(Subscript[(P^Ω), refl])]"},ImageSize->Large,PlotPoints->200],
{{ϕrf,(π/4),"RF Phase ϕrf [rads]"},0,(π/2)},
{{Γ,0.1,"RF Modulation Depth Γ [rads]"},0,0.5},
{{T1,0.1,"Input Mirror Transmission (Subscript[T, 1])"},0,1},
{{T2,0.1,"End Mirror Transmission (Subscript[T, 2])"},0,1},
{{ζ,0,"Demod Angle ζ [rads]"},-π,π},
Paneled->False
]