NLO
NLO
NLO
S OMMERSEMESTER 2021
U LF P ESCHEL
2
1 Phenomenological representation of non-resonant
nonlinearities
We start with M AXWELL’s equations in the optical domain without external charges, currents
and no magnetization
⃗ ∂B ⃗
∇×E =− ∇·D =0 (1.1)
∂t
⃗ ∂D ⃗
∇ × B = µ0 ∇ · B = 0, (1.2)
∂t
where E ist the electric field, B the magnetic induction and D the dielectric displacement
which is defined as
D = ε0 E + P , (1.3)
consisting of a vacuum contribution ε0 E and the polarization P of the material. We note that
M AXWELL’s equations are linear in all components.
The polarization of the material P is the source of nonlinearities and is driven by the elec-
tromagnetic fields. It cannot be deduced from electrodynamics, but we need quantum me-
chanics/solid state physics. Furthermore it has to obey symmetries which will be discussed
in this chapter.
We define the polarization P (r , t ) at a fixed point in space r and time t . In principle it de-
pends on electromagnetic fields acting on all other places r ′ and other times t ′ as P (r , t ) =
P [E (r ′ , t ′ ), B (r ′ , t ′ )]. However, we make several assumptions to simplify our calculations:
• There are only local interactions: r = r ′ . We explicitly exclude the propagation of a
polarization in the material. This means the polarization of one unit of the material
is influence by light acting on just that unit only. The reason for that is that nonlocal
effects on atomic scales take place at ∼ 10−10 m while the variations of the optical field
are ∼ 10−6 m.
• Causality: The polarization cannot depend on fields from the future t ′ ≤ t .
• No external action on the polarization within the memory time of the material. The
polarization depends on time differences only P (r , t ) = P [E (r , t − τ), B (r , t − τ)] with
τ > 0.
• No magnetic effects P (r , t ) = P [E (r , t − τ)]
The above assumptions are well justified for conventional materials. Unfortunately, a further
simplification is not straightforward.
3
1.2 The Taylor expansion of the polarization NLO
The general idea ist that for non-resonant nonlinearities the interaction between photons
via the material is weak. This implies that only a very limited number of photons hits one
atom of the material within the relaxation time/memory of the polarization. We can classify
the polarization with respect to the number of interacting photons
The linear polarization P (1) (t ) is induced by a single photon and is by far the dominant po-
larization. P (2) (t ) is the nonlinear polarization of quadratic/second order and is induced by
two-photon processes.
3
ˆ∞
P i(1) (t ) = ε0 dτ1 R i(1) (τ1 )E j (t − τ1 ),
X
j
(1.5)
j =1
0
where R i(1)
j
(τ1 ) is the linear response function. It is a time dependent second rank tensor
which is decaying rapidly in time.
3
ˆ∞ ˆ∞
P i(2) (t ) = ε0 dτ2 R i(2) (τ , τ )E j (t − τ1 )E k (t − τ2 ).
X
dτ1 jk 1 2
(1.6)
j ,k=1
0 0
4
1.3 Polarization in frequency space NLO
ˆ∞ ˆ∞
3
Jˆ[P i(2) (t )] = Jˆε0 dτ2 R i(2) (τ , τ )E j (t − τ1 )E k (t − τ2 )
X
dτ1 jk 1 2
i , j =1
0 0
3
ˆ∞ ˆ∞
= ε0 dτ2 JˆR i(2) (τ , τ ) JˆE j (t − τ1 ) JˆE k (t − τ2 )
X
dτ1 jk 1 2
i , j =1
0 0
3
ˆ∞ ˆ∞
−P i(2) (t ) = ε0 dτ2 R i(2) (τ , τ )(−1)E j (t − τ1 )(−1)E k (t − τ2 )]
X
dτ1 jk 1 2
j ,k=1
0 0
−P i(2) (t ) = P i(2) (t ). (1.7)
This implies that P i(2) (t ) must be equal to zero for all times t which can only be fulfilled with
R i(2)
jk
= 0. In a similar way we can prove that all even order nonlinearities must vanish in cen-
trosymmetric and amorphous materials.
In optics we usually deal with monochromatic fields and are more interested in the spectral
amplitudes than in the fast temporal evolution. Often it is required to quantify the nonlinear
response on certain frequency components of the electric field.
a.) We start with the decomposition of the optical field into its frequency components via
Fourier transform which are defined in this lecture as follows
ˆ∞ ˆ∞
1
E (t ) = dω Ẽ (ω)e −iωt
and Ẽ (ω) = dt E (t )eiωt . (1.8)
2π
−∞ −∞
The temporal field E (t ) is thus a sum of oscillating fields each having a spectral amplitude
Ẽ (ω). Note that E (t ) is real valued which implies that Ẽ (ω) is symmetric
ˆ∞ ˆ∞
∗
1 1
Ẽ (−ω) = dt E (t )e−iωt = dt E (t )eiωt = Ẽ ∗ (ω). (1.9)
2π 2π
−∞ −∞
5
1.3 Polarization in frequency space NLO
This means that the spectral amplitudes at negative frequencies Ẽ (−|ω|) are nonzero, but
strictly correlated to their positive counterparts and thus need not be considered separately.
b.) Now we can find an expression for P̃ (1) (ω) by inserting (1.8) into (1.5)
3
ˆ∞
P i(1) (t ) = ε0 dτ1 R i(1) (τ1 )E j (t − τ1 )
X
j
j =1
0
3
ˆ∞ ˆ∞
= ε0 dτ1 R i(1) dω Ẽ j (ω)e−iω(t −τ1 )
X
j
(τ1 )
j =1
0 −∞
ˆ∞ 3
ˆ∞ ˆ∞
dω ε0 dτ1 R i(1) (τ1 )eiωτ1 Ẽ j (ω)e−iωt = dω P̃ i(1) (ω)e−iωt .
X
= j
(1.10)
j =1
−∞ 0 −∞
P̃ i(1) (ω) is the spectral amplitude of the linear polarization. It depends on fields oscillating
with the same frequency in a linear way
ˆ
∞
3 3
P̃ i(1) (ω) = ε0 dτ1 R (1) (τ1 )eiωτ1 Ẽ j (ω) = ε0 χ(1)
X X
ij ij
(−ω|ω)Ẽ j (ω). (1.11)
j =1 j =1
0
3
ˆ∞ ˆ∞
P i(2) (t ) = ε0 dτ2 R i(2) (τ , τ )E j (t − τ1 )E k (t − τ2 )
X
dτ1 jk 1 2
(1.13)
i , j =1
0 0
3
ˆ∞ ˆ∞ ˆ∞ ˆ∞
= ε0 dτ1 dτ2 R i(2) (τ , τ ) −iω1 (t −τ1 )
dω2 Ẽ k (ω2 )e−iω2 (t −τ2 )
X
jk 1 2
dω1 Ẽ j (ω1 )e
i , j =1
0 0 −∞ −∞
ˆ∞ ˆ∞ 3
ˆ ∞ ˆ∞
dω2 ε0 dτ2 R i(2) (τ , τ )ei(ω1 τ1 +ω2 τ2 ) Ẽ j (ω1 )Ẽ k (ω2 )e−i(ω1 +ω2 )t .
X
= dω1 dτ1 jk 1 2
i , j =1
−∞ −∞ 0 0
6
1.3 Polarization in frequency space NLO
ˆ∞
1
P̃ i(2) (ωp ) = dt eiωp t P i(2) (t )
2π
−∞
ˆ∞ ˆ∞ ˆ∞
1 h
= dt eiωp t dω1 dω2
2π
−∞ −∞ −∞
3
ˆ∞ ˆ∞ i
ε0 dτ2 R i(2) (τ , τ )ei(ω1 τ1 +ω2 τ2 ) Ẽ j (ω1 )Ẽ k (ω2 )e−i(ω1 +ω2 )t
X
dτ1 jk 1 2
i , j =1
0 0
ˆ∞
1
= dt ei(ωp −ω1 −ω2 )t ...
2π
−∞
ˆ∞ ˆ∞ ˆ∞ ˆ∞
3
dω2 ε0 dτ2 R i(2) (τ , τ )ei(ω1 τ1 +ω2 τ2 ) Ẽ j (ω1 )Ẽ k (ω2 ).
X
dω1 dτ1 jk 1 2
i , j =1
−∞ −∞ 0 0
(1.14)
ˆ∞ ˆ∞ 3
P̃ i(2) (ωp ) = χ(2) (−ωp |ω1 , ω2 )Ẽ j (ω1 )Ẽ k (ω2 )δ(ωp − ω1 − ω2 )
X
dω1 dω2 i jk
j ,k=1
−∞ −∞
ˆ∞ 3
dω1 ε0 χ(2) (−ωp |ω1 , ωp − ω1 )Ẽ j (ω1 )Ẽ k (ωp − ω1 )
X
= i jk
(1.16)
j ,k=1
−∞
ˆ∞ ˆ∞
χ(2)
i jk
(−ωp |ω1 , ω2 ) = dτ1 dτ2 R i(2) (τ , τ )ei(ω1 τ1 +ω2 τ2 ) .
jk 1 2
(1.17)
0 0
7
1.3 Polarization in frequency space NLO
• Because the spectral amplitudes of the electric field can also have arguments with neg-
ative frequencies, equation (1.16) incorporates difference as well as sum frequency
generation.
d.) Higher order nonlinear polarizations are introduced in the same way as to derive equa-
tion (1.16)
ˆ∞ ˆ∞ 3 3 h
dωn ε0 χ(n)
αp α1 ... αn (−ωp |ω1 , ..., ωn )Ẽ α1 (ω1 )
X X
P̃ αp (n)(ωp ) = dω1 ... ...
α1 =1 αn =1
−∞ −∞
n i
... Ẽ αn (ωn )δ ωp − ωk )
¡ X
(1.18)
k=1
• The susceptibility carries the symmetries of the respective material. The higher the
symmetry of the material the lower the number of independent coefficients.
• A K RAMERS -K RONIG transformation can be performed with respect to one of the fre-
quencies. It links real and imaginary parts of the coefficients as
ˆ∞ χ(n)
1 αp α1 ... αn (−ωp (ϖ)|ω1 , ...ϖ...ωn )
χ(n)
αp α1 ... αn (−ωp (ω0 )|ω1 , ...ω0 ...ωn ) = P dϖ
iπ ϖ − ω0
−∞
The quantity |E λ | is the peak amplitude of the electric field component with frequency ωλ .
Fourier spectrum of sum of cw-waves is equal to a sum of delta functions
1 X£
E λ δ(ω − ωλ ) + E λ∗ δ(ω + ωλ ) .
¤
Ẽ (ω) = (1.21)
2 λ
8
1.3 Polarization in frequency space NLO
The nonlinear polarization of n-th order induced by cw-waves is a sum of harmonic oscilla-
tions
1 X ³ (n) −iωp t ´
P (n) (t ) = P ωp e + c.c. . (1.22)
2 ωp
We now want to determine the amplitudes of harmonic oscillations of the nonlinear polar-
ization of n-th order induced by cw-waves P ω(n)
p
by following these steps:
1. Formal integration of equation (1.18) using the Fourier spectrum from (1.21).
2. Inserting the nonlinear coefficient χ(n)
αp α1 ... αn (−ωp |ω1 , ..., ωn ) defined in equation (1.19).
3
(∗)
P ω(n) = ε0 K (−ωp | ± ω1 , ... ± ωn ) χ(n) (∗)
αp α1 ... αn (−ωp | ± ω1 , ... ± ωn )E 1,α1 E n,αn
X X
p αp
(1.24)
{±ωλ } α1 ... αn =1
• ±ωλ : one of the frequencies of the interacting cw-waves, either positive or negative
(∗)
• E λ,α : vector component αλ of the field amplitude E λ of the continuous wave oscillat-
λ
ing with ωλ . In case of the negative frequency the complex conjugate field profile E λ,α
∗
λ
is chosen.
• {±ωλ }: all combinations of interacting ±ωλ including negative frequencies, which gen-
erate ωp , but no permutations.
The factor K (−ωp |±ω1 , ...±ωn ) incorporates the effect of permutations. It includes all factors
1
2 used to express the amplitudes of the electric fields and the polarization. We can generalize
it in the following way:
9
1.4 Symmetries and conserved quantities NLO
χαp α1 ... αk ... αn (−ωp |ω1 . . . ωk ωn ) = χαk α1 ... αp ... αn (ωk |ω1 . . . −ωp ωn ). (1.26)
This assumption is justified far away from resonance which means that nonlinear coeffi-
cients can be expressed by a Drude model
1 n 1
χ(n)
Y
αp α1 ... αn (−ωp |ω1 ...ωn ) ∼ , (1.27)
ω20 − ω2p λ=1 ω20 − ω2λ
where ω0 is the resonance frequency of the material. The frequency of the polarization does
not play an extraordinary role and can be mixed with the other components.
The d -coefficients are a very common representation of quadratic nonlinearities. They are
orientation dependent (x- and y-axis usually aligned parallel to the crystal axis). We can also
write down the K -factor for all cases possible:
• All frequencies different and nonzero: K = 1
• ω1 = ω2 or ω1 = −ω2 : K = 12
• ω1 = 0 or ω2 = 0: K = 2
10
1.5 Cubic nonlinearities in isotropic materials NLO
For non-degenerate processes |ωα | ̸= |ωβ | and arbitrary α and β including the induced po-
larization we have photon flux conservation
Wα Wβ
− = const. (1.29)
ωα ωβ
Cubic nonlinearities can usually be neglected if quadratic nonlinearities are present. There-
fore they are only relevant if quadratic coefficients vanish (e. g. in inversion symmetric ma-
terials as glass) or if quadratically nonlinear processes are suppressed (e. g. due to the lack of
phase matching). The susceptibility tensor
χ(3)
i j kl
(−(ω1 + ω2 + ω3 )|ω1 , ω2 , ω3 ) (1.30)
possesses 34 = 81 elements.
From now on will will only consider isotropic materials with the highest possible degree of
symmetry, which means they have the minimum number of independent elements.
E = (E j + E k + E l )êk .
Thus, the induced polarization must point into the same direction as E because the material
has no internal direction. Since E ||P we have
χ(3)
i kkk
= 0 for i ̸= k and χ(3) (3) (3)
xxxx = χ y y y y = χzzzz . (1.31)
The three non-vanishing elements must be equal as no direction is preferred. Thus we only
have one independent element.
11
1.5 Cubic nonlinearities in isotropic materials NLO
E = (E j + E k )êk + E l êl .
The vector components of P (3) at a fixed frequency are P i(3) ∼ χ(3) E E E . We assume that
i kkl j k l
the material is rotated around the singular axis (l ) by 180°. In the new coordinate system the
electric field is
As the material is isotropic the nonlinear coefficient has not changed by rotation χ(3)
i kkl
=
′
χ(3)
i kkl
. The vector components of P (3) in the rotated coordinate system are
′
P i(3) ∼ χ(3)
i kkl
(−E j )(−E k )E l = P i(3) . (1.32)
This is only possible if P (3) points along the axis of rotation. It follows that χ(3)
i kkl
= 0 for i ̸= l .
Thus we have 18 non-vanishing coefficients (all permutations are taken into account, for
different frequencies permutations give different coefficients)
χ(3) (3)
x y y x (−(ω1 + ω2 + ω3 )|ω1 , ω2 , ω3 ) ̸= χxx y y (−(ω1 + ω2 + ω3 )|ω1 , ω2 , ω3 ). (1.33)
For this process we consider the case j = k ̸= l = i with the nonlinear coefficient
χ(3)
i j kl
(−3ω, ω, ω, ω).
Since the arguments are equal the frequencies can be exchanged freely
(3)
χ(3) (3) (3)
xx y y = χx y x y = χx y y x = ... = χi i kk . (1.35)
12
1.5 Cubic nonlinearities in isotropic materials NLO
The relation between P (3) and E must be independent of the rotation. Therefore we can
equate (1.37) and (1.38)
3 (3) 1
χ y zz y + χ(3) = χ(3)
yyyy. (1.39)
2 2 yyyy
This means in an isotropic material THG is determined by a single independent coefficient
(3) (3)
χ(3) (3)
xx y y = χx y x y = ... = χTHG and χ(3) (3) (3)
xxxx = χ y y y y = χzzzz = 3χTHG . (1.40)
We can generalize the total third harmonic polarization by first looking at a single compo-
nent
1
P y(3) = ε0 (χ(3) 3 (3) 2 (3) 2
y y y y E y + 3χ y zz y E y E z + 3χ y xx y E y E x )
4
3 ³ ´
= ε0 χ(3) E 2
+ E 2
+ E 2
x Ey
4 THG y z
3
= ε0 χ(3) E 2E y . (1.41)
4 THG
Thus the total third harmonic polarization is
3
P (3) = ε0 χ(3) E 2E . (1.42)
4 THG
The third harmonic generated in an isotropic medium:
13
1.5 Cubic nonlinearities in isotropic materials NLO
We want to give an explanation for this phenomenon. In an isotropic medium the direction
of polarization is entirely determined by the direction of the exciting field. In case of circu-
larly polarized electric field the polarization can only rotate with the rotation speed of the
fundamental harmonic, hence one rotation per optical cycle of the fundamental harmonic.
Therefore, only polarization at fundamental harmonic frequency is generated.
χ(3) (3)
xx y y (−ω|ω, ω, −ω) ̸= χx y y x (−ω|ω, ω, −ω, ). (1.43)
For SPM we have three frequency permutations (p = 3) and therefore K = 3/4 (see equa-
tion (1.25)). We use the same arguments as for THG
1 (3)
χ(3) (3) (3)
y y y y = (χ y zz y + 2χ y y zz + χ y y y y ) ⇒ χ(3) (3) (3)
y y y y = χ y zz y + 2χ y y zz . (1.45)
2
h i
P y(3) = ε0 K 2χ(3)
y y zz (|E x | 2
+ |E y | 2
+ |E z | 2
)E y + χ (3)
y zz y |E y | 2
E y + χ (3)
y zz y (E 2
x + E 2 ∗
z )E y
h i
= ε0 K 2χ(3) 2 2 2 (3) 2 2
y y zz (|E x | + |E y | + |E z | )E y + χ y zz y (E x + E y + E z )E y .
2 ∗
(1.46)
3 ³ ´
P (3) (ω) = ε0 2χ(3) 2 (3) 2 ∗
y y zz |E | E + χ y zz y E E . (1.47)
4
We find two independent coefficients and want two discuss both nonlinear actions:
• First term: |E |2 E
depends on the intensity only and is a pure phase modulation or a nonlinear loss
14
1.5 Cubic nonlinearities in isotropic materials NLO
• Second term: E 2 E ∗
is sensitive to the state of polarization of the electric field and completely vanishes for
circularly polarized light as E 2 (see end of section 1.5.1). It bears half of the contribu-
tion as the first term for linearly polarized light.
We can further simply the result for K LEINMANN symmetry (no frequency dependence):
χ(3) (3)
y zz y = χ y y zz . Only a single coefficient is left and we obtain
3
P (3) = ε0 χ(3) 2 2 ∗
¡ ¢
y y zz 2|E | E + E E . (1.48)
4
9 9 (3)
P (3) = ε0 χ(3) 2 2
y y zz |E | E = ε0 ∆χ with ∆χ = χ y y zz |E | . (1.49)
4 4
Thep
action of the nonlinear susceptibility can be expressed by a change of the refractive index
n = 1 + χ as
q ∆χ 9 (3)
∆n = n 02 + ∆χ − n 0 ≈ = χ |E |2 , (1.50)
2n 8n 0 y y zz
where n 0 is the refractive index in the absence of a nonlinear action. As there is no gen-
erally accepted definition of the amplitude of an oscillating field one prefers to refer to the
intensity
ε0
I= cn 0 |E |2 . (1.51)
2
We can replace the electric field in (1.50) and find
(3)
9 χ y y zz
∆n = n 2 I with n 2 = . (1.52)
4 ε0 cn 02
m2
silica n 2 = 2,3 · 10−20 (1.53)
W
2
m
Al0.18 Ga0.82 As(1,5 µm) n 2 = 2,1 · 10−17 . (1.54)
W
15
2 Resonant nonlinearities - a quantum mechanical
representation
In what follows we will shortly discuss the kind of nonlinear response, which shows up close
to a resonance of the material. If the energy of the incident photons just roughly matches
the energy of a quantum mechanical transition, the optical properties of the excited mate-
rial may change dramatically as a function of the incident intensity. A Taylor expansion of
the induced polarization with respect to powers of the acting electrical field as done in equa-
tion (1.4) is still applicable, however, the treatment applied in the last section is no longer
useful. The expansion can usually not be restricted to a few orders and resulting coefficients
turn out to be highly dispersive. Hence, a treatment taking into account the material re-
sponse more seriously has to be applied.
1
Ĥ = (p̂ + e A real )2 + V (r ). (2.1)
2m
The term V (r ) is the potential, which localizes the electrons. A real is the vector potential,
which represents the optical field by
∂
E real = − A real . (2.2)
∂t
The subscript real is used to emphasize that observable fields are real valued and that they
contain components oscillating at +ω0 and −ω0 . We consider the vector potential as classical
which means that
• The strong field contains many photons
• quantum fluctuations can be neglected
• The field commutes with all operators
2
• The vector potential is weak compared to atomic fields ⇒ A real is neglected.
16
2.1 Equations of motion of a two-level-system NLO
We can now decompose the Hamiltonian into an unperturbed part Ĥ0 and a perturbation
Ŵ
1 2 e e 2 2
Ĥ = p̂ + V (r ) + p̂ A real + A . (2.3)
|2m {z } |m {z } 2m real
Ĥ0 Ŵ
The unperturbed part Ĥ0 is the Hamiltonian of the unperturbed two-level-system in absence
of any electromagnetic ¯ field.
® It defines ® energy levels E a and follows the SCHRÖDINGER
the
equation Ĥ0 ψa = E a ψa , where ψa/b is a orthonormal basis of eigenstates ψa |ψb ¯ψa |ψb =
¯ ® ¯ ¯ ®
¯ ¯ ¯
δab . The evolution is restricted to the two eigenstates, which means that the system is either
in state a or b. The wave function of the whole system can be expressed as
¯ψ(t ) = a(t ) ¯ψa + b(t ) ¯ψb with |a(t )|2 + |b(t )|2 = 1,
¯ ® ¯ ® ¯ ®
(2.4)
where the second conditions follows from ψ(t )|ψ(t )¯ψ(t )|ψ(t ) = 1.
¯ ®
The interaction operator Ŵ drives the nontrivial dynamics and induces changes. Using the
time dependent S CHRÖDINGER equation leads to
d ¯¯
ψ(t ) = Ĥ ¯ψ(t )
® ¯ ®
iℏ
dt
d¡
a(t ) ¯ψa + b(t ) ¯ψb = Ĥ (a(t ) ¯ψa + b(t ) ¯ψb )
¯ ® ¯ ®¢ ¯ ® ¯ ®
iℏ
dt
¯ ® d ¯ ® d
i ℏ ¯ψ a a(t ) + iℏ ¯ψb b(t ) = a(t ) Ĥ ¯ψa + b(t ) Ĥ ¯ψb .
¯ ® ¯ ®
(2.5)
dt dt
Multiplying from the left with ψa/b ¯ and using ψa |ψb ¯ψa |ψb = δab leads to
¯ ¯ ®
d
a = E a a + ψa |Ŵ |ψa a + ψa |Ŵ |ψb b
® ®
iℏ
dt
d (2.6)
iℏ b = E b b + ψb |Ŵ |ψb b + ψb |Ŵ |ψa a .
® ®
dt | {z }
interaction elements
The last integral just gives us the values at the boundary. For a decaying wave functions the
boundaries are zero, leading to vanishing interaction terms. In the last line we assume a real
valued ψa (r ).
17
2.1 Equations of motion of a two-level-system NLO
®∗
The second matrix element ψa |Ŵ |ψb = ψb |Ŵ |ψa is responsible for transitions between
®
the two states. It is usually nonzero , provided that an external field is present and both states
have different parity. For the evaluation we need a helpful commutator relation
· ¸ · ¸
1 2 1 2
[r̂ , Ĥ0 ] = r̂ p̂ + V (r ) − p̂ + V (r ) r̂
2m 2m
1 1 iℏ
= (r̂ p̂ 2 − p̂ 2 r̂ ) = ([r̂ , p̂]p̂ + p̂[r̂ , p̂]) = p̂. (2.8)
2m 2m m
Now we can calculate the matrix element as
e ® e D m E
ψa |Ŵ |ψb = A real ψa |p̂|ψb = A real ψa | (r̂ Ĥ0 − Ĥ0 r̂ )|ψb
®
m m iℏ
ie
= A real (E a − E b ) ψa |r̂ |ψb .
®
(2.9)
ℏ
The optical field couples to the states via the dipole moment of the transition
ˆ
ψa |r̂ |ψb = dV ψ∗a (r )r ψb (r ).
®
(2.10)
dI d¡ 2 ¢ (2.4) d ¡ d
|a| − |b|2 = 2|a|2 − 1 = 2 |a|2
¢
=
dt dt dt dt
d 2
= 2a ∗ a + c.c. = ψa |Ŵ |ψb a ∗ b + c.c.
®
dt iℏ
2e
= 2 A real (E a − E b ) ψa |r̂ |ψb a ∗ b + c.c.
®
(2.11)
ℏ
Now we want to do a frequency analysis of ab ∗ for a weak perturbation Ŵ and a resonant
excitation (E a − E b )/ℏ ≈ ω0 . If the energy of the state is much larger than the interaction
energy the solutions of a(t ) and b(t ) using (2.6) are
µ ¶ µ ¶
iE a iE b
a(t ) ∼ exp − , b(t ) ∼ exp − . (2.12)
ℏt ℏt
This leads to the following expression for ab ∗
E a − Eb
µ ¶
∗
ab ∼ exp −i t ≈ exp(−iω0 t ), (2.13)
ℏ
which means that ab ∗ is related to the optical polarization.
18
2.1 Equations of motion of a two-level-system NLO
b.) Polarization
The classical change of inversion is a result of absorbed optical power
dI absorbed power per system 2 1 ∂
=2 = E real P real . (2.14)
dt energy required to excite a system E a − E b N ∂t
This can be understood by noting that power is Voltage times current, where E real corre-
∂
sponds to the voltage and ∂t P real to current. The latter one can be understood that the polar-
ization describes a deviation of the electrons from their rest position in the lattice. A change
of polarization describes the electron oscillations around their rest position and thus a cur-
∂
rent. Then E real ∂t P real is the absorbed power per volume. It contains fast parts (oscillations
with ≈ 2ω0 ), however, they average out or can be neglected. In the following we only consider
the slow/almost stationary parts.
We decompose the fields into fast parts oscillating with the carrier frequency of the optical
field ω0 and slowly varying amplitudes
1£
P (t ) exp(−iω0 t ) + P ∗ (t ) exp(iω0 t )
¤
P real (t ) =
2
1£
E real (t ) = E (t ) exp(−iω0 t ) + E ∗ (t ) exp(iω0 t )
¤
(2.15)
2
1£
A real (t ) = A(t ) exp(−iω0 t ) + A ∗ (t ) exp(iω0 t ) .
¤
2
The amplitudes are assumed to be slow that their time derivatives are negligible compared
with the derivatives of the fast-varying exponents
¯∂ ¯ (2.2) 1
|ω0 A(t )| ≫ ¯ A(t )¯ ⇒ A(t ) ≈ E (t ) (2.16)
∂t iω0
¯∂ ∂ iω0 £
P (t ) exp(−iω0 t ) − P ∗ (t ) exp(iω0 t ) .
¯ ¤
|ω0 P (t )| ≫ ¯ P (t )¯ ⇒ P real (t ) ≈ − (2.17)
∂t ∂t 2
Then we can formulate the classical evolution of inversion
dI (2.14) 2 1 ∂
= E real P real
dt E a − Eb N ∂t
∗
2 1 E (t ) −iω0
= P (t ) + c.c. + fast oscillations
E a − Eb N 2 2
ω
= E ∗ (t )P (t ) + c.c. + fast oscillations. (2.18)
2i(E a − E b )N
Similarly we can formulate the quantum evolution of inversion
dI (2.11) 2e
A real (E a − E b ) ψa |r̂ |ψb a ∗ b + c.c.
®
=
dt ℏ2
(2.13) e i
(E a − E b ) ψb |r̂ |ψa E ∗ (t ) exp(iω0 t )ab ∗ + c.c. + fast oscillations.
®
= (2.19)
(2.15) ℏ ω0
2
If we now compare the slow parts containing E ∗ (t ) in the classical and quantum formula-
tions
e i ω0
(E a − E b ) ψb |r̂ |ψa E ∗ (t ) exp(iω0 t )ab ∗ ≈ E ∗ (t )P (t ).
®
(2.20)
ℏ ω0
2 2i(E a − E b )N
19
2.1 Equations of motion of a two-level-system NLO
E a − Eb 2
µ ¶
ψb |r̂ |ψa ab ∗ exp(iω0 t ).
®
P = −2eN (2.21)
ℏω0
| {z }
≈1
Here we assume a resonant excitation which means that the frequency detuning ∆ω is small
E a − Eb
∆ω = − ω ≪ ω0 . (2.22)
ℏ
The total or real valued polarization then is
1£
P (t ) exp(−iω0 t ) + c.c. = −eN ψb |r̂ |ψa ab ∗ + c.c.
¤ ®
P real = (2.23)
2
Now we also want to discuss the dynamics of the nonlinear polarization by using (2.6)
db ∗
µ ¶ µ ¶
d ¡ ∗¢ da ∗
iℏ ab = iℏ b − a iℏ
dt dt dt
= (E a − E b )ab − ψa |Ŵ |ψb |a|2 − |b|2
∗
®¡ ¢
ie
= (E a − E b )ab ∗ − A real (E a − E b ) ψa |r̂ |ψb I
®
ℏ
1 £ ¤ ie
≈ (E a − E b )ab ∗ − E (t )e−iω0 t − E ∗ (t )eiω0 t (E a − E b ) ψa |r̂ |ψb I
®
2iω0 ℏ
e(E a − E b )
= (E a − E b )ab ∗ − E (t )e−iω0 t − E ∗ (t )eiω0 t ψa |r̂ |ψb I .
£ ¤ ®
(2.24)
2ℏω0
The driving term of ab ∗ with the wrong frequency (∼ e+iω0 t ) is neglected. The evolution of
the polarization can be derived by using (2.24), (2.21) and (2.22)
dP Ne 2
= ∆ωP + I ψb |r̂ |ψa ( ψa |r̂ |ψb E ).
® ®
i (2.25)
dt ℏ
For an isotropic medium with linear polarization of the electric field, the polarization P is
parallel to E which leads to
Equations of motion
d Ne 2
P = ∆ωP + | ψb |r̂ |ψa |2 E I
®
i (2.27)
dt ℏ
dI 1
= E ∗ P + c.c. (2.28)
dt 2iℏN
Note that up to now the equations of motion contain no damping, because we treated the
optical field classically thus omitting spontaneous emission. Equations (2.27) and (2.28)
are strictly valid only in time domains which are small compared with classical relaxation
times.
20
2.2 Rabi oscillations NLO
ϑ
Furthermore we have the conservation of the norm
|a(t )|2 + |b(t )|2 = 1 and no interest in absolute phases Im b
which means we can set a to a real value. Now we may ϕ
plot the space of the two-level state as a function of
Re b
a, Re(b) and Im(b) on the surface of the Bloch-sphere
(sphere with radius 1).
We can now express the norm in terms of slowly varying
envelope of the polarization and inversion
(|a|2 + |b|2 )2 = |a|4 + |b|4 − 2|a|2 |b|2 + 4|a|2 |b|2 = (|a|2 − |b|2 )2 + 4|ab ∗ |2
1
= I2 + 2 2 ® |P |2 = 1 (2.30)
e N | ψb |r̂ |ψa |2
and find this conserved quantity connecting inversion and polarization. Here we used
Inversion and the scalar component of the normalized slowly varying envelope of the polar-
ization are on the surface of the Bloch sphere.
In what follows we solve the system of equations (2.27) and (2.28) for a linearly polarized
optical field oscillating with the frequency of the transition ∆ω = 0 with constant amplitude,
but a step-like envelope (rectangular pulse)
where Θ(t ) is the H EAVISIDE step function. We assume that for the initial conditions (t = −∞)
the system is purely in the ground state I = −1 and P = 0.
d d
I = 0 ⇒ I = −1, P = 0 ⇒ P = 0. (2.32)
dt dt
The system remains in its ground state.
21
2.2 Rabi oscillations NLO
b.) Now we look at 0 < t < t 1 with |E | = E 0 = const. (real valued). We can decompose the
polarization into its real and imaginary parts
P = [P R (t ) + iP i (t )]ê. (2.33)
d
P R = 0 ⇒ P R (t ) = 0
dt
d Ne 2
| ψb |r̂ |ψa |2 E 0 I
®
Pi = − (2.34)
dt ℏ
d E0
I= Pi . (2.35)
dt ℏN
Differentiating the third equation we find an equation for a harmonic motion
d2 e
I + Ω2 I = 0 Ω = E 0 | ψb |r̂ |ψa |
®
Rabi frequeny. (2.36)
dt 2 ℏ
The whole system performs oscillations with the Rabi frequency Ω. The solution of the sys-
tem is then
ℏN dI ℏN
I (t ) = − cos(Ωt ), Pi = = Ω sin(Ωt ). (2.37)
E 0 dt E0
d ℏN
P i = 0 ⇒ P i (t ) = Ω sin(Ωt 1 ) stationary envelope
dt E0
d
I = − cos(Ωt 1 ). (2.38)
dt
ê £
P i (t ) exp(−iω0 t ) − P i∗ (t ) exp(iω0 t ) .
¤
P real (t ) = i (2.39)
2
After the end of the excitation the polarization keeps on oscillating forever. But, any real
system is damped, an effect we have omitted by describing the optical field classically. Hence
above solution is only valid in a time range short compared with the relaxation time, which
can range from 100 fs in solids up to ms in dilute gases.
We conclude the results from the solution:
• Polarization and electric field need not to be proportional.
• The maximum strength of the polarization need not depend on the electric field strength
(see equation (2.37)).
• The polarization can be present, although the optical field is absent.
22
2.3 The influence of damping NLO
I /P
1 I
t
t1
Fig. 1: Visualization of the the inversion (blue) and polarization (orange) for an electric field present
at 0 ≤ t ≤ t 1 .
Damping is always present in real world systems. It is caused by the coupling to a continuum
of modes, which could be modes of the electromagnetic radiation, phonons or plasmons. Its
exact calculation requires a quantum mechanical description of the respective field. Here
we just introduce some phenomenological constants into respective evolution equations.
Ne 2
µ ¶
d i
i P = ∆ω − | ψb |r̂ |ψa |2 E I
®
P+ (2.40)
dt T2 ℏ
dI 1 I +1
= (E ∗ P − E P ∗ ) − , (2.41)
dt 2iℏN T0
where T2 describes the relaxation time of the polarization and T0 the life time of the upper
state. Note that the relaxation term of the inversion was chosen such that I relaxes towards
the equilibrium state −1.
We can describes different scenarios using these relaxations constants:
• Free atom/atom in a trap:
Here, relaxation is due to spontaneous emission with rather long relaxation times (ms).
The relaxation is strictly correlated T2 = 2T0 .
23
2.3 The influence of damping NLO
• Dense systems:
These are gases under normal pressure, liquids or solids which will be considered in
the following. They show a strong mutual interaction with many decay channels. The
polarization is dominantly affected by so-called “dephasing”. Therefore T2 is orders of
magnitude shorter than T0 , e. g. in a direct semiconductor: T2 ∼ 100 fs, T0 ∼ ns.
We can now distinguish different cases according to the actual pulse duration T p :
We can now neglect the time derivative of P (t ) which leads to an algebraic equation for the
polarization
Ne 2
µ ¶
i
0 = ∆ω − | ψb |r̂ |ψa |2 E I
®
P+
T2 ℏ
Ne 2
| ψb |r̂ |ψa |
®2
⇒ P (t ) = − ℏ I E (t ). (2.43)
∆ω − Ti2
The polarization is determined by the actual optical field. Now we can also find the effective
inversion dependent susceptibility
Ne 2 | ψb |r̂ |ψa |2
®
P (t ) = ε0 χ(t )E (t ) with χ(I ) = − I. (2.44)
ε0 ℏ ∆ω − Ti2
=: χ0
| {z }
The polarization has lost its independence. It has become enslaved. We can now also de-
scribe the evolution of the inversion
dI 1 I +1 ε0 I +1
= (E ∗ P − E P ∗ ) − = Im(χ(I ))|E |2 − (2.45)
dt 2iℏN T0 ℏN T0
2 | ψ |r̂ |ψ |2
®
Ne b a 1
with Im(χ(I )) = − I ∼ I.
ε0 ℏ ∆ω + 2 T2
2 1
T2
24
2.3 The influence of damping NLO
We can rewrite the differential equation as a rate equation of Im[χ(t )] by multiplying (2.45)
with the prefactor of Im[χ(t )]
d ε0 Im(χ) − Im(χ0 )
Im(χ) = − Im(χ0 ) Im(χ)|E |2 − , (2.46)
dt | ℏN {z }| T
{z0 }
(1) (2)
with χ0 being the susceptibility in the ground state. The first term (1) describes absorption
(Im χ0 > 0) induced bleaching of the transition by reducing Im χ. It switches off, when Im χ →
0. The second term (2) describes the relaxation towards the ground state.
The complete susceptibility induced by the two-level-system is then
∆ω2 + T12
χ= 2
i
T2 Im(χ). (2.47)
∆ω − T2
Finally lets look at the situation when the optical field is constant after the Rabi oscillations
have decayed (t → ∞):
• We are in a steady state with Im χ > 0.
• Almost zero inversion (I → ∞) can be reached where the system is almost transparent.
A complete transparency is impossible, as absorption is switched off close to trans-
parency. Carrier relaxation represented by T0 prevents the system from reaching zero
inversion.
• A state of positive inversion cannot be attained.
Real and imaginary part of the susceptibility are proportional to each other (c. f. equation (2.47))
χ0 (ω) Nℏc n
χ(ω) = IS = . (2.49)
1 + IIS0 2T0 Im(χ0 )
This is the simplest formulation describing the response of an optical transition. It is only
valid for time scales much longer than all relaxation times of the material. It describes a
single transition, to which the optical field is resonant. We only require two quantities for
the whole description:
25
2.3 The influence of damping NLO
• The ground state susceptibility χ0 (ω) of the two-level-system (usually only a small frac-
tion of the complete susceptibility)
• The saturation intensity I S .
Equation (2.48) is not equivalent to the expressions derived in chapter (1.5.2) (see e. g. equa-
tion (1.52)). It is not only highly dispersive, but also a Taylor expansion with respect to the
optical intensity yielding such quantities as e. g. n 2 will only converge in a limited range of
intensities. If the expansion is done at I 0 = 0
I0
χ ≈ χ0 − χ 0 (2.50)
IS
it will diverge for intensities higher than I s and will give wrong results much earlier.
26
3 Nonlinear optics and field propagation - a plane
wave approach
In the last chapter we have learned about the nonlinear response of the material. Together
with M AXWELL’s equations we can now in principle determine any kind of field propagation
in the presence of a nonlinear response. Unfortunately, the complete set of equations in-
cluding all nonlinear interaction is much too difficult to be solved in a self-consistent way.
Hence, approximations have to be made which are usually based on the assumption that
the nonlinear response of the material is small compared with the linear one and that some
peculiar features of the geometry allow for further simplifications. A few examples of such
treatments will be presented in that chapter.
• The material polarization at a particular frequency is split into a linear and a nonlinear
part representing the linear part with the dielectric tensor ε̂
Remark: An index profile as it is typical for waveguides can be treated in a similar man-
ner by taking waveguide modes instead of plane waves into account.
• Small perturbation due to nonlinear action: |P NL (r )| ≪ |ε0 ε̂(r )E (r )|.
We want to briefly state M AXWELL’s equations for monochromatic fields
⃗
∇ × E = iωB , ⃗ ∇ × B = −iωµ0 (ε0 ε̂E + P NL ), ⃗
∇ · ε0 ε̂E + P NL = 0, ⃗
¡ ¢
∇ · B = 0.
ω2
µ ¶
⃗ 1 1
∇⃗
³ ´
∇ · E − ∆E = 2 ε̂E + P NL with ε0 µ0 = 2
c ε0 c
ω 2 ¶
ω 1
2
µ
∆ −⃗
∇div + 2 ε̂ E = − 2 P NL . (3.4)
c c ε0
| {z } | {z }
linear operator small perturbation
27
3.2 Propagation in z-direction only NLO
First we want to look at solutions E L (r ) of the linear system. Then we can write E L (r ) as an
eigenvector of the eigenvalue problem
ω2
ω2
ε̂−1 (∆ − ⃗
∇ div)E L (r ) = − 2 E L (r ) with eigenvalue − . (3.5)
c c2
For a homogeneous infinite system we find the solution E L (r ) = E k eik·r . This is the starting
point to determine the nonlinear solution via perturbation theory.
Remark: One might be tempted to drop the term ⃗ ∇ · E L , but this is only possible in homo-
geneous isotropic media. In anisotropic media as they are typical for a quadratic nonlinear
response, the term ⃗
∇ · E L is not necessarily zero and longitudinal electric fields may occur
εx
0 0
∂E x ∂E y ∂E z
ε̂ = 0
εy 0 ⇒⃗
∇ · (ε̂E L ) = εx + εy + εz =0
∂x ∂y ∂z
0 0 εz
⃗ ∂E x ∂E y ∂E z
=⇒
̸ ∇ · EL = + + = 0. (3.6)
∂x ∂y ∂z
In what follows we discuss several cases for which further assumptions can be made.
In many relevant case (e. g. frequency conversion) very broad collimated beams, which
travel into a single direction are used and any dependence on transverse coordinates (x, y)
can be neglected or taken into account parametrically only. We assume for our solution to
take the form
The related magnetic induction is obtained by the linear solution E β exp iβz via the curl
¡ ¢
⃗
∇ × (E β exp iβz ) = (iβêz ) × (E β exp iβz ) = iωB β exp iβz
¡ ¢ ¡ ¢ ¡ ¢
β
Bβ = êz × E β . (3.9)
ω
28
3.2 Propagation in z-direction only NLO
Now we want to insert our ansatz (3.8) into (3.4) where we use that u(z) is a slowly varying
envelope with
¯ ¯
¯d ¯
¯ u(z)¯ ≪ |βu(z)|. (3.10)
¯ dz ¯
This leads us to
ω2
µ ¶
⃗
= ∆ − ∇ div + 2 ε̂ [E k u(r ) exp(ik · r )]
c
d2 ω2 β
· ¸
β β
= (E − E z êz ) 2 + 2 ε̂E (u(z) exp iβz )
¡ ¢
dz c
2
ω2 β ¢ d2
· ¸ µ ¶
β β d β β d
= u(z) (E − E z êz ) 2 + 2 ε̂E exp iβz + (E − E z êz ) exp iβz
¡ ¢ ¡
+ 2iβ u(z)
dz c dz 2 dz
ω2
µ ¶ µ ¶
β ¢d d
⃗ β β
= u(z) ∆ − ∇ div + 2 ε̂ E exp iβz + (E − E z êz ) exp iβz
¡ ¢ ¡
+ 2iβ u(z)
c dz dz
| {z }
=0 linear solution
ω2 1
µ ¶
β β d d
+ 2iβ u(z) = − 2 P NL .
¡ ¢
= (E − E z êz ) exp iβz (3.11)
dz dz c ε0
Assuming the amplitude to vary slowly we can simplify
µ ¶
d
+ 2iβ u(z) ≈ 2iβu(z) (3.12)
dz
and multiplying a suitable vector from the left
β d ω2 1
(E β − E z êz ) exp iβz 2iβ u(z) = − 2 P NL
¡ ¢
dz c ε0
β
ω2 exp −iβz (E β − E z êz )∗ NL
¡ ¢
d
i u(z) = − 2 β
P . (3.13)
dz c 2βε0 |E β − E êz |2
z
q
W
As the amplitude u(z) has currently no unit, we change the unit to m2
that the squared
amplitude represents the intensity flowing in z-direction
rD E
β
a(z) = s z u(z) (3.14)
β
D E
with s z being the z-component of the time averaged P OYNTING vector of the linear solu-
tion which we want to calculate in the following:
D E
β 1
sz = Re(E × B ∗ )z
2µ0
1
Re[E β iβz × B β∗
¡ ¢ ¡ ¢
= exp exp
−iβz ]z
2µ0
β 1 β h β2
· µ ¶¸
1 β β∗
i
= Re E × êz × E = Re |E | êz − E β (E β∗ êz )
2µ0 ω z 2µ0 ω z
1 β h
β
i 1 β¯ β ¯
β ¯2
¯
= |E β |2 − |E z |2 = ¯E − E z êz ¯ . (3.15)
2µ0 ω 2µ0 ω
29
3.3 Diffraction in linearly isotropic media NLO
We can summarize our findings of the nonlinearly driven evolution of a mode in z-direction
Evolution in z-direction
β
d ω2 1 (E β − E z êz )∗ NL ¡ ¢
i u(z) = − 2 P exp −iβz (3.16)
dz c 2βε0 |E β − E β êz |2
z
β
s
β
d ω3 (E − E z êz )∗ NL ¡ ¢
i a(z) = − P exp −iβz . (3.17)
dz 8βε0 c |E β − E β êz |
2
z
To describe effects of a nonlinearly induced change of the beam profile one has to take into
account the effect of transverse coordinates to account for diffraction, self-focusing and spa-
tial soliton formation. To simplify the analysis, we restrict to isotropic media, which are most
relevant for cubic nonlinearities.
The starting point is again equation (3.4) with a scalar valued dielectric function ε̂ = ε
ω ´ ω2 1
∆ −⃗
³
∇ div + 2 ε E = − 2 P NL . (3.18)
c c ε0
Again we make the assumption of a weak nonlinear action |P NL | ≪ |ε0 εE | and homogeneity
⃗
∇ ε = 0 and use the first of M AXWELL’s equations
0 =⃗
∇·D =⃗
∇ · ε0 εE (r ) + P NL
£ ¤
≈⃗
∇ · [ε0 εE (r )] = ε0 ε⃗
∇ · E (r ) = 0. (3.19)
ω1 ω2 1
µ ¶
∆ + 2 E = − 2 P NL . (3.20)
c ε c ε0
¯ ∂u ¯
¯ ¯
β
¡ ¢
E (r ) = E u(x, y, z) exp iβz with ¯¯ ¯¯ ≪ |β u(x, y, z)|, (3.21)
∂z
where u(x, y, z) is the slowly varying envelope. E β is the polarization of a linear wave propa-
gating in z-direction which is transverse (no z-component).
30
3.4 Pulse evolution in linearly isotropic media NLO
ω2
µ ¶
= ∆ + 2 ε [E β u(x, y, z) exp iβz ]
¡ ¢
c
à !
β 2 ∂ ∂2 ∂2 ∂2 ω2
= E exp iβz −β ε u(x, y, z)
¡ ¢
+ 2iβ ∂z + ∂z 2 + ∂x 2 + ∂y 2 + c
2
∂2 ∂2
· µ ¶ ¸
β
¡ ¢ d d
= E exp iβz + 2iβ + + u(x, y, z)
dz dz ∂x 2 ∂y 2
∂2 ∂2 ω2 1
· ¸
β d
+ 2 + 2 u(x, y, z) = − 2 P NL .
¡ ¢
= E exp iβz 2iβ (3.22)
dz ∂x ∂y c ε0
We again assumed a slow varying amplitude. Similarly to section 3.2 we can multiply a suit-
able vector from the left and rearrange
1 ∂2 ∂2 ω2 1 ¢ E β∗ NL
· µ ¶¸
d ¡
i + + u(x, y, z) = − exp −iβz P . (3.23)
dz 2β ∂x 2 ∂y 2 c 2 2βε0 |E β |2
rD E
β
Furthermore we can apply the transformation a = s z u to arrive at intensities based on
β
D E
the z-component of the time averaged P OYNTING vector of the linear solution s z (c. f. (3.15))
D
β
E 1 β β2
sz = |E | . (3.24)
2µ0 ω
β
Note that E z is zero for isotropic media. We now obtain the evolution equation for nonlin-
earity driven diffraction
To describe effects of a nonlinearly induced change of a pulse profile one has to take into
account a finite frequency spectrum and cannot restrict to monochromatic waves anymore
to account for dispersive spreading, self-compression and temporal soliton formation. To
simplify the analysis, we restrict to isotropic media and a propagation into a single direction
thus neglecting diffractive spreading.
We start with equation (3.20) for monochromatic waves propagating in an isotropic medium
with ε̂ = ε under the action of a weak nonlinear action at a fixed frequency
³ ω ´ ω2 1
∆ + 2 ε E (ω, z) = − 2 P NL (ω, z). (3.26)
c c ε0
31
3.4 Pulse evolution in linearly isotropic media NLO
As most of the quantities used in the further derivation are frequency dependent, we restrict
to a fixed frequency ω0 being in the middle of the investigated spectrum and representing
the carrier frequency of the pulse.
We assume that the field shape is similar to that of a linear wave propagating at ω0 for all
ω2
relevant frequency components E L (ω0 , r ) = E β0 exp iβ0 z with β20 = ε(ω0 ) c 20 . The nonlinear
¡ ¢
ω2 ∂ ω2
µ ¶ · ¸
β0 iβ0 z 2
∆ + 2 ε(ω) E (ω, z) ≈ E e 2iβ0 + ε(ω) − β0 u(ω, z)
c ∂z c 2
ω2 1
= − 2 P NL (ω, z). (3.28)
c ε0
Now we can multiply a suitable vector from the left and rearrange
h ∂ β0 ∗
1 ω2 ω2 1 −iβ0 z E
µ ¶i
2
i + ε(ω) − β 0 u(ω, z) = − e β
P NL (ω, z). (3.29)
∂z 2β0 c 2 c 2β0 ε0
2
|E 0 |2
| {z }
(1)
2
We prepare for a F OURIER-back-transform under the assumption that β(ω) = ωc 2 ε(ω) is close
to β0
β(ω) + β0 1 GVD
(1) = (β(ω) − β0 ) ≈ β(ω) − β0 = (ω − ω0 ) + (ω − ω0 )2 + ... (3.30)
2β0 vg 2
where we use the expressions of group velocity v g and group velocity dispersion GVD
¯ ¸−1
∂ d2
· ¯
β(ω) ¯¯ β(ω)
¯ ¯
vg = , GVD = ¯ . (3.31)
∂ω ω0 dω2 ¯
ω0
32
3.4 Pulse evolution in linearly isotropic media NLO
h ˆ ∞
i
β0 −i(β0 zω0 t )
E real (z, t ) = Re E dω u(ω, z)e−i(ω−ω0 )t e
| {z } | {z }
0 Fourier transform of SVE carrier wave
h i
= Re E β0 ũ(z, t )ei(β0 z−ω0 t ) , (3.33)
ˆ∞
ũ(z, t ) = dω u(ω, z)e−i(ω−ω0 )t . (3.34)
0
Now we can apply the F OURIER transform (3.34) on equation (3.29) using approximation (3.30)
ˆ∞
∂
· ¸
1 GVD
dω i + (ω − ω0 ) + (ω − ω0 ) + ... u(ω, z)e−i(ω−ω0 )t
2
∂z v g 2
0
∂ i ∂ GVD ∂2
· ¸
= i + − + ... ũ(t , z)
∂z v g ∂t 2 ∂t 2
ˆ∞
1 E β0 ∗
dω −ω2 P NL (ω, z)e−iωt e−i(β0 z−ω0 t )
¡ ¢
=
2β0 ε0 c 2 |E β0 |2
0
β0 ∗
E ∂2 NL
· ¸
1
= P̃ (z, t ) e−i(β0 z−ω0 t ) . (3.35)
2β0 ε0 c 2 |E β0 |2 ∂t 2
Remark: P̃ NL (z, t ) is not the real valued nonlinear polarization, but its complex counterpart,
the real part of which is the direct measurable polarization. It contains the positive frequency
part only, but not its complex conjugate counterpart on the negative frequency side.
rD E
β
Again we can apply the transformation a(z, t ) = s z 0 ũ(z, t ) to express the evolution in
β
D E
terms of intensities. As before s z 0 is the z-component of the time averaged P OYNTING
vector of the linear solution at the center frequency of the pulse ω0 . Again, we arrive at the
evolution equation for the slowly varying envelope of a pulse
· β0 ∗ 2
∂ i ∂ GVD ∂2 E ∂
· ¸ ¸
1
i + − + ... a(t , z) = p β0 |2 ∂t 2
P̃ (z, t ) e−i(β0 z−ω0 t ) .
NL
∂z v g ∂t 2 ∂t 2 8β0 c ε0 ω0
2 |E
(3.36)
33
4 Field evolution for quadratic nonlinearities
The main subject of this chapter will be Second Harmonic generation (SHG). The basic prin-
ciple is that a fundamental harmonic (FH) wave at ω is frequency doubled to a second har-
monic (SH) wave oscillating at 2ω.
The geometry consists of plane waves propagating in z- ω 2ω
direction (no transverse structure). Note that the orien-
tation of the coordinates is defined by the propagation χ(2)
ω ω
direction and not by the nonlinear crystal. Because typ-
ical tensors of the nonlinear response are given with re-
spect to the crystal axis, an additional coordinate transformation might be required.
3
P i(2) (ω) = ε0 K (−ω|2ω, −ω) χ(2) (−ω|2ω, −ω)E 2ω j E ω∗k
X
i jk
with K (−ω|2ω, −ω) = 1
j ,k=1
3 1
P i(2) (2ω) = ε0 K (−2ω|ω, ω) χ(2) (−2ω|ω, ω)E ω j E ωk with K (−2ω|ω, ω) = .
X
i jk
(4.1)
j ,k=1 2
The interacting waves are one at FH frequency and one at SH frequency. We adapt (3.17) and
obtain
ω∗
s
d ω3 E⊥ NL
(z, ω) exp −iβ(ω)z
¡ ¢
FH i a ω (z) = − 2 ω∗ P
dz 8β(ω)ε0 c |E ⊥ |
s (4.2)
2ω∗
d ω3 E⊥ NL
¡ ¢
SH i a 2ω (z) = − P (z, 2ω) exp −iβ(2ω)z .
dz β(2ω)ε0 c 2 |E ⊥
2ω∗
|
ω/2ω
Here E ⊥ = E ω/2ω − E zω/2ω êz is the transverse field polarization at FH/SH frequency.
The total electric field driving the nonlinear polarization can be found by using (3.8), (3.14)
and (3.15)
a ω (z)
E ω (r ) = E ω u ω (z)eiβ(ω)z = E ω q ® eiβ(ω)z
s zω
Eω
s
2ω
= ω a ω (z)eiβ(ω)z
|E ⊥ | ε0 c 2 β(ω)
s
E 2ω 4ω
E 2ω (r ) = 2ω a 2ω (z)eiβ(2ω)z . (4.3)
|E ⊥ | ε0 c β(2ω)
2
34
4.1 Equations of motion NLO
We can now insert the fields (4.3) and the polarization (4.1) into the evolution equations (4.2)
s
d ω3 e−iβ(ω)z X ω∗ (2)
i a ω (z) = − ω E P i (ω)
dz 8β(ω)c 2 ε0 |E ⊥ | i =x,y i
s
ω3 e−iβ(ω)z X ω∗ X 3
=− ω E χ(2) (−ω|2ω, −ω)E 2ω j E ω∗k
8β(ω)c 2 |E ⊥ | i =x,y i j ,k=1 i j k
s
ω3 e−iβ(ω)z X ω∗ X 3
=− 2 ω Ei χ(2)
i jk
(−ω|2ω, −ω)
8β(ω)c |E ⊥ | i =x,y j ,k=1
E 2ω E ω∗
s s
j 4ω iβ(2ω)z k 2ω
· 2ω a 2ω (z)e ω a ω (z)∗ e−iβ(ω)z
|E ⊥ | ε0 c β(2ω)
2 |E ⊥ | ε0 c β(ω)
2
ω∗ ω∗ E 2ω∗
s
ω5 3 E ⊥i E j
i[β(2ω)−2β(ω)]z
χi j k ω kω
(2) ∗
X
= −e a 2ω a ω . (4.4)
β(ω) β(2ω)ε0 c i , j ,k=1
2 6 2ω
|E ⊥ | |E ⊥ | |E ⊥ |
| {z }
=:χ(2)
FH
Now we assume that no further fields are involved and energy is conserved (Manley-Rowe
relations) which leads to the conservation of energy flux
µ ¶ µ ¶
d¡ 2 2
¢ ∗ d ∗ d
0= |a ω | + |a 2ω | = ia ω −i a ω + ia 2ω −i a 2ω + c.c.
dz dz dz
³ ´ ³ ´
∗ −i∆βz (2)
= ia ω e ∗
χFH a 2ω (z)a ω ∗
(z) + ia 2ω ei∆βz χ(2)
SH ω
a 2
(z) + c.c.
= ie−i∆βz χ(2) ∗2
a (z)a ω
FH 2ω
(z) − ie−i∆βz χ(2)∗ a a ∗ 2 (z) + c.c.
SH 2ω ω
h i
∗ 2 (2)
= −2 Im e−i∆βz a 2ω a ω (χFH − χ(2)∗
SH
) . (4.7)
35
4.2 Normalization NLO
Note that (4.7) must hold for every z. This can only be fulfilled by
χ(2)
FH
= χ(2)∗
SH
=: χ(2)
eff
. (4.9)
d
FH: i a −∆βb + χ(2)
eff
a∗b = 0 (4.11)
dz
d
SH: i b − ∆βb + χ(2)∗
eff
a 2 = 0. (4.12)
dz
4.2 Normalization
We now aim towards a further simplification of the evolution equations for SHG which might
provide us with simple estimations of relevant scales. For that we introduce dimensionless
quantities:
z = Z0 Z , a = I 0 A and b = I 0 B eiϕ ,
p p
(4.13)
where Z0 is the characteristic length, I 0 the relevant intensity level and ϕ a constant phase
term. We insert these definitions into (4.11) and (4.12) and find
d
A −∆βZ0 B + Z0 χ(2)
p iϕ ∗
FH: i eff
I0e A B = 0 (4.14)
dZ
d
B − ∆βZ0 B + Z0 χ(2)
p −iϕ 2
SH: i eff
I 0 e A = 0. (4.15)
dZ
Now we can actually choose the normalization constants Z0 and I 0 in such a way, that all
constants become unity
p
1 I0
Z0 = , χ(2)
eff
eiϕ = 1. (4.16)
|∆β| |∆β|
Z0 now describes length scales for which the FH and SH waves run out of phase. For a per-
fectly phase-matched crystal this length scale goes to infinity. For the intensity we obtain
¯ ∆β ¯2
¯ ¯
1
I0 = = ¯ (2) ¯ . (4.17)
¯ ¯
(2)
|Z0 χ | ¯ χ ¯
eff eff
36
4.3 Conservation laws - the Hamiltonian NLO
The characteristic induces a noticeable nonlinear action (it is reduced by a growing non-
linear coefficient, but growing with the squared mismatch). This explains why second har-
monic generation is such a rare phenomenon. Particular care has to be taken to reduce the
mismatch. Otherwise required power levels grow to infinity.
The phase ϕ is chosen as
³ ´
ϕ = −arg χ(2)
eff
. (4.18)
It removes the phase of the nonlinear coefficient. This means the phase of the coefficient of
the quadratic nonlinearity has no physical relevance.
Remark: Normalization is not a unique procedure. For example, in case of vanishing mis-
match the sample length is much smaller than the inverse mismatch therefore becoming
the relevant scale.
The final set of normalized equations can be summarized to
d
FH: i A −σB + A ∗ B = 0 (4.19)
dZ
d
SH: i B − σB + A 2 = 0. (4.20)
dZ
σ = sgn(∆β) is called the scaled mismatch and can have values −1, 0, +1. The original fields
for a given A(Z ) and B (Z ) are
Eω
s
z iβ(ω)z−iωt
µ ¶
p h 2ω
E real (z, t ) = I 0 Re ω A Z = e
β(ω)ε0 c 2 |E ⊥ | Z0
s
E 2ω z 2iβ(ω)z−2iωt +iϕ i
µ ¶
4ω
+ B Z= e . (4.21)
β(2ω)ε0 c 2 |E ⊥2ω
| Z0
H = σB ∗ B − A 2 B ∗ − A ∗2 B. (4.23)
37
4.4 General solution NLO
d d d d
H = σB ∗ B − A 2 B ∗ − 2AB ∗ A + c.c. (4.24)
dZ dz dz dz
= iσB ∗ (−σB + A 2 ) + iA 2 (−σB + A 2 )∗ − 2iAB ∗ A ∗ B + c.c.
= −iσ|B |2 + iσ|A|4 − 2i|A|2 |B |2 + ∗2 ∗2
iσB A −
iσB A + c.c.
= 0. (4.25)
Since the other terms are purely imaginary they will cancel with they complex conjugated
terms. The quantum relation of the conserved quantities is the conserved norm of the wave
function which is equivalent to total intensity and the expectation value of the Energy (Hamil-
tonian operator) which is equivalent to the Hamiltonian.
Now we want to find a general solution for the set of normalized evolution equations for SHG.
We start with a transformation of complex variables into real ones by separating amplitudes
and phases
giving us four independent variables. We can insert this ansatz (4.26) into equations (4.19)
and (4.20) and divide by eiϕ A (Z ) and eiϕB (Z ) which yields
d ∂
FH i |A(Z )| − |A(Z )| ϕ A (Z ) + |A(Z )| · |B (Z )|eiϕB (Z )−2iϕ A (Z ) = 0
dZ ∂Z (4.27)
d ∂
SH i |B (Z )| − |B (Z )| ϕB (Z ) − σ|B (Z )| + |A(Z )|2 e2iϕ A (Z )−iϕB (Z ) = 0.
dZ ∂Z
Now we only consider the imaginary part of equation (4.27)
d
|A(Z )| + |A(Z )| · |B (Z )| sin ∆ϕ(Z ) = 0
¡ ¢
FH
dZ (4.28)
d
|B (Z )| − |A(Z )|2 sin ∆ϕ(Z ) = 0,
¡ ¢
SH
dZ
where ∆ϕ(Z ) = ϕB (Z ) − 2ϕ A (Z ) is the phase difference between SH and FH fields. We want
to investigate the evolution equation of the phase difference as well by using the real part of
equation (4.27)
d ∂ ∂
∆ϕ(Z ) = ϕB (Z ) − 2 ϕ A (Z )
dZ ∂Z ∂Z
|A(Z )|2
µ ¶
− 2|B (Z )| cos ∆ϕ(Z ) .
¡ ¢
= −σ + (4.29)
|B (Z )|
In all coherent systems without a fixed time reference the absolute phase does not matter
and only a phase difference plays a role. Thus we only have three independent real valued
variables ∆ϕ(Z ), |A(Z )|, |B (Z )|.
38
4.5 Phase-matched SHG NLO
We can reduce the number of free variables further by using the energy conservation (4.22)
d
|B (Z )| = Q − |B (Z )|2 sin ∆ϕ(Z )
¡ ¢ ¡ ¢
dZ
(4.30)
Q − |B (Z )|2
µ ¶
d
∆ϕ(Z ) = −σ + − 2|B (Z )| cos ∆ϕ(Z ) .
¡ ¢
dZ |B (Z )|
We can once more reduce the number of free variables by using the conserved Hamilto-
nian (4.23)
We can use this result to extract and eliminate sin ∆ϕ(Z ) = 1 − cos2 (∆ϕ(Z )) from equa-
¡ ¢ p
tion (4.30)
s
¶2
H − σ|B |2
µ
d 2 2
|B | = (Q − |B | ) sin ∆ϕ = ±(Q − |B | ) 1 −
¡ ¢
. (4.32)
dZ 2|B |(Q − |B |2 )
Now we transform to the SH power I = |B |2 by multiplying equation (4.32) with 2|B | because
dI
dZ
= 2|B | d|B
dZ
|
v
¶2
H − σI
u µ
d u
I = ±(Q − I ) 4I − 4I p
t
dZ 2 I (Q − I )
q
= ± 4I (Q − I )2 − (H − σI )2 . (4.33)
The formal solution can be obtained via integration. Using the initial conditions for given
SH power at Z = 0 : I (0) = I0
ˆ (L)
I ˆL
dI
p =± dZ = ±L. (4.34)
4I (Q − I )2 − (H − σI )2
I0 0
A practicable solution of equation (4.34) still requires solving the integral and inverting the
resulting function. This finally leads to so-called elliptic functions.
A nonlinear system with N real valued variables requires N − 1 conserved quantities or sym-
metries for an analytical solution. The existence of conserved quantities is essential for an
analytical solution.
Now we want to find solutions of (4.34) for the case of complete phase matching (σ = 0) and
no SH power at Z = 0
|B (Z = 0)|2 = I0 = 0. (4.35)
39
4.6 Down conversion (phase matched case) NLO
Looking at the conserved quantities we find that the total energy is Q = |A(Z = 0)|2 and the
Hamiltonian is simply zero H = 0 by comparing with (4.23). This simplifies (4.34) dramati-
cally and leads to
ˆ (L)
I
dI
p = ±L,
2(Q − I ) I
0
p !¯
Q + I ¯¯I (L)
Ãp
1
⇒ p ln p p = ±L. (4.36)
2 Q Q− I 0
¯
We conclude that total conversion is only reached for infinity propagation length. The con-
version efficiency η = I
p
Q depends only on the product QL. The required power scales in-
verse to the squared propagation length (see normalization in chapter 4.2).
1
SH power I
0.5
0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
p (2) −1
propagation length in ( P in |χeff |)
40
4.7 Stability of the second harmonic field NLO
Only the initial value of I (I = Q instead of I = 0) has changed compared with SHG. The
derivation is analogous to section 4.5 but we have changed the initial conditions
p !¯
Q + I ¯¯I (L)
Ãp
1
p ln p p = ±L. (4.40)
2 Q Q− I Q
¯
This leads to some problems. The function on the left side of equation (4.40) becomes in-
finity on the lower boundary. The intensity of the SH wave can only deviate from Q, if the
propagation length L tends to infinity. This is consistent with an inversion of the SH process
discussed above. The latter one also needs an infinite propagation length to obtain complete
conversion.
A direct solution of equation (4.36) yields
d p
I = ±2(Q − I ) I ⇒ I (L) = Q. (4.41)
dZ
The system can stay in a state with all the power in the SH field and formally requires an
infinite propagation length to obtain noticeable down-conversion.
1
SH power I
0.5
0
−5 −4 −3 −2 −1 0 1 2 3 4 5
p (2) −1
propagation length in ( P in |χeff |)
Fig. 3: Down-conversion for the phase-matched case. The down conversion process starts from
quantum fluctuations.
In real systems any small perturbation (e. g. quantum noise) starts the down-conversion pro-
cess after a finite length. In the phase-matched case complete down-conversion is obtained
(inverse behaviour to SHG) and SHG starts (see previous chapter).
We have already noticed that down-conversion is started by field noise. To investigate that
phenomenon in more detail we perform a stability analysis. We take the results from sec-
tion 4.2
d
FH: i A −σB + A ∗ B = 0
dZ (4.42)
d 2
SH: i B − σB + A = 0.
dZ
41
4.7 Stability of the second harmonic field NLO
p −iσZ
We want to look for stationary solutions B 0 = Qe and A 0 = 0 with small fluctuations
δA(Z ) and δB (Z )
σ
³p ´
A(Z ) = δA(Z )e−i 2 Z and B (Z ) = Q + δB (Z ) e−iσZ . (4.43)
d σ ³p ´
FH: i δA + δA + δA ∗ Q + δB
=0
dZ 2 (4.44)
d
SH: i δB δA2
+ = 0.
dZ
We assume δA and δB to be small, so that quadratic terms can be neglected. Then the SH
perturbation δB remains
¡ in first¢ order approximation a constant (as long as δA(Z ) remains
small). Now we apply −i dz + σ2 to the first part of (4.44)
d
d σ σ
µ ¶·µ ¶ ¸
d ∗
p
0 = −i + i + δA + δA Q
dz 2 dZ 2
¶ ¸∗
σ 2 d σ
µ 2 ³ ´ ¶ ·µ
d p
= + δA + Q i + δA
dz 2 2 dz 2
σ 2
µ 2 ³ ´ ¶
(4.44) d
= + −Q δA. (4.45)
dz 2 2
the general solution of this equation can be expressed in terms of exponential functions
Ãr ! Ã r !
³ σ ´2 ³ σ ´2
δA(Z ) = δA + exp Q − Z + δA − exp − Q − Z . (4.46)
2 2
The two coefficients δA ± are fixed and determined by the initial fluctuations δA(Z = 0) and
their first derivative
σ
¯
d
= i δA(Z = 0) + i QδA ∗ (Z = 0).
p
δA(Z ) ¯¯
¯
(4.47)
dZ Z =0 2
d
r ³ σ ´2
δA(Z = 0) = δA + + δA − and δA(Z = 0) = Q− [δA + − δA − ] (4.48)
dZ 2
to find an expression for δA ±
à !
σ/2
p
1 Q
δA ± = ±ip δA(Z = 0) ± i p δA ∗ (Z = 0). (4.49)
2 4Q − σ2 4Q − σ2
42
4.7 Stability of the second harmonic field NLO
³ σ ´2
a.) SH power is below a threshold Q ≤
2
This occurs for a mismatched sample (σ = ±1) only: Q ≤ 41 . The arguments in the exponential
terms become imaginary
Ãr ! Ã r !
1 1
δA(Z ) = δA + exp i −Q Z + δA − exp −i −Q Z . (4.50)
4 4
This leads to an oscillatory behaviour. The initial perturbation is conserved (which is typ-
ical for Hamiltonic systems). We do not observe growth of the SH wave or further down-
conversion. The system is (marginally) stable.
³ σ ´2
b.) SH power is above a threshold Q >
2
Then the fluctuations will grow exponentially
Ãr ! Ã r !
³ σ ´2 ³ σ ´2
δA(Z ) = δA + exp Q − Z + δA − exp − Q − Z . (4.51)
2 2
| {z }
≈ 0 for large Z
"Ã ! # Ãr !
σ/2 ³ σ ´2
p
1 Q
+ip δA(Z = 0) + i p δA ∗ (Z = 0) exp Q − Z . (4.52)
2 4Q − σ2 4Q − σ2 2
We observe that this solution is unstable. Also the initial perturbation δA and its complex
conjugate δA ∗ are amplified (so-called parametric amplifier).
We also want to discuss the case of large SH power with 4Q ≫ σ2 = 1. Then we can approxi-
mate δA(Z ) as
1¡ ¢ p
δA(Z ) ≈ δA(Z = 0) + iδA ∗ (Z = 0) e Q Z . (4.53)
2
The phase of the amplified perturbation for large Z is
Im[δA(Z = 0) + iδA ∗ (Z = 0)]
µ ¶
ϕ = arctan
Re[δA(Z = 0) + iδA ∗ (Z = 0)]
π
µ ¶
Im[δA(Z = 0)] + Re[δA(Z = 0)]
= arctan = arctan(1) = . (4.54)
Re[δA(Z = 0)] + Im[δA(Z = 0)] 4
This means that after a certain propagation length the phase of the perturbation is fixed and
independent of δA(Z = 0). This means the phase fluctuations (even of quantum origin) are
suppressed which leads to squeezed (non-classical) light.
Lastly we want to mention the critical intensity for down conversion in real world units Q <
1
4
¯ ¯2
³ σ ´2 1 ¯¯ ∆β ¯¯
I thresh. = I 0 = ¯ (2) ¯ . (4.55)
2 4¯χ ¯
eff
43
4.8 Visualization of the field dynamics using the Hamiltonian NLO
We aim for a demonstration of the field dynamics in an intuitive way. We start again with
the conserved quantities, namely the Hamiltonian H given in (4.31) and the total intensity Q
given by (4.22). Then we can find another conserved quantity
where |B ′ | is the scaled amplitude of the SH field with 0 ≤ |B |′ ≤ 1, ∆ϕ the phase difference
between SH and FH fields and σ′ the rescaled mismatch.
The system now evolves along the lines/contours of fixed rescaled Hamiltonian H ′ . The con-
served direction of motion on each contour line is determined by
Q
µ ¶
d
∆ϕ(Z ) = −σ + − 3|B (Z )| cos ∆ϕ(Z )
¡ ¢
dZ |B (Z )|
µ ¶
1 d ′ 1 ′
∆ϕ(Z ) = −σ + − 3|B | cos ∆ϕ(Z ) .
¡ ¢
p ′
(4.57)
Q dZ |B |
We can draw the contour lines for several values of the Hamiltonian H ′ and different mis-
matches in figure 4.
Fig. 4: Contour plots of the scaled Hamiltonian H ′ for zero mismatch (left) and σ′ = 0.5 (right). The
radial coordinate |B ′ | describes the scaled SH power. For our initial conditions we start at con-
tour lines going through the origin.
For SHG we find a contour line of H ′ = 0 which starts at the origin (at zero SH power). For a
zero mismatch 4 (left) we see that |B ′ | approaches 1 (complete conversion) after infinite time.
However, if we introduce a mismatch 4 (right) then the inversion is incomplete, because the
44
4.9 Low-depletion-limit NLO
1
0.5
H′
H′
0
0
−0.5 1 1
0.5 ¡ ¢ 0.5 ¡ ¢
−1 −0.5 0 0
sin ∆ϕ −1 −0.5 0 sin ∆ϕ
¡ 0 ¢ 0.5 −0.5
1 −1 ¡ ¢ 0.5 −0.5
1 −1
cos ∆ϕ cos ∆ϕ
Fig. 5: 3D visualization of the scaled Hamiltonian for zero mismatch (left) and σ′ = 0.5 (right).
system evolves on a contour line which never reaches the border of |B ′ | = 1. We observe a
periodic motion corresponding to up- and down-conversion.
If we start outside the origin (e. g. simultaneously inject FH and SH fields), the phase between
SH and FH fields matters. We see a periodic motion around the minimum or maximum of
the scaled Hamiltonian. If we start at one of the extrema of H ′ we will have a stationary state
consisting of SH and FH fields.
The scaled mismatch σ′ = pσ controls the dynamics of the system. Increasing the power
Q
corresponds to a reduction of mismatch.
4.9 Low-depletion-limit
We now want to consider a system with a large mismatch with the unscaled equations (4.12).
We will find a low conversion efficiency, which means that the FH intensity is nearly constant
(undepleted pump approximation). Then
p
a(z) ≈ I FH = const. (4.58)
Then we may remove a(z) from the equation of the SH field (4.12) yielding
d
SH: i b − ∆βb + χ(2)∗ I = 0.
eff FH
(4.59)
dz
Using the initial condition of no SH intensity b(z = 0) = 0 we find the solution
χ(2)∗
eff ¡ ¢
b(z) = − I FH (exp −i∆βz − 1). (4.60)
∆β
The intensity of the SH wave is then
2
|χ(2)∗
eff
|2 ¯ ∆βz ¯2 ¯ ∆βz
2 ¯ i 2 ¯ ¯ −i 2
∆βz ¯2
¯
I SH = |b(z)| = I FH − ei 2 ¯
∆β2
¯e ¯ ¯e
´ 2
∆β
³
sin 2 z
= 4|χ(2)∗ | 2 2
I FH
. (4.61)
eff ∆β
45
4.9 Low-depletion-limit NLO
I SH
distance z in 1/∆β
1 2 3 4 5 6 7 8 9 10 11 12
L
Fig. 6: SH power in the undepleted pump approximation with a(z) = const. L describes the oscillation
period of SH-generation.
We observe (figure 6) that the intensity decays quickly with growing mismatch and oscillates
periodically with a period of
¯ ¯
¯ ¯
λ λ
¯ ¯ ¯ ¯ ¯ ¯
2π ¯ 2π ¯ ¯ 1 ¯ ¯ FH SH ¯
L= = ¯¯ ¯=¯
n n
¯=¯ ¯
|∆β| β(2ω) − 2β(ω) ¯
¯ SH − 2 FH ¯
¯ ¯ ¯ λFH n SH − 2λSH n FH ¯
¯λ λFH ¯
SH
¯ λSH
¯ ¯
¯
= ¯¯ ¯, (4.62)
n −n ¯
FH SH
where n FH , n SH are the refractive indices at FH/SH frequency and λFH , λSH the vacuum wave-
length of the FH/SH wave.
46
5 Diffraction, dispersion and nonlinearities
Now we want to discuss effects of further (transverse) coordinates which leads us to partial
differential equations instead of ODEs. Usually we will find no analytical solutions. However,
we need the transverse coordinates in order to account for diffraction.
The monochromatic field implies that no Third-Harmonic Generation takes place and we
only have to deal with self-phase modulation. Since the electric field is linearly polarized, the
nonlinearity induces a refractive index change ∆n = n 2 I (I -intensity). Then the nonlinear
induced polarization takes the form
∂ 1 ∂2 ∂2 ω
· µ ¶ ¸
2
i + + + n 2 |a(x, y, z)| a(x, y, z) = 0. (5.4)
∂z 2β ∂x 2 ∂y 2 c
47
5.1 Beam evolution in isotropic Kerr materials NLO
We can now choose the normalization constants in such a way that the prefactors become
one
Z0 2πnW02
=1 ⇒ Z0 = βW02 =
βW02 λ0
ω λ0 1 λ0 λ0 λ0 2
µ ¶
1
n 2 Z0 I 0 = 1 ⇒ I0 = = = . (5.6)
c 2π |n 2 |Z0 2π 2πn|n 2 |W02 n|n 2 | 2πW0
We see that W0 is a free scaling parameter, so we choose it e. g. as the initial beam width. This
means that the principal effects we discuss will not depend on the beam width.
Since we chose W0 to be the beam width its advantageous to introduce the total beam power
ˆ∞ ˆ∞ ˆ∞ ˆ∞
2
P= dx dy |a(x, y, z)| = W02 I 0 dX dY |A(X , Y , Z )|2 , (5.7)
−∞ −∞ −∞ −∞
µ ¶2
1 λ0
with P 0 = W02 I 0 = as a new scaling quantity.
n|n 2 | 2π
Note that P 0 only depends on material properties and wavelengths, not the beam width.
Qualitative changes are now to be expected to depend on the beam power in relation to P 0 ,
but not on the beam width, as it can be scaled away.
The final equation we have to solve is a nonlinear 2+1 dimensional S CHRÖDINGER equation
(NLS)
∂ 1 ∂2 ∂2
· µ ¶ ¸
2
i + + + κ|A| A = 0 with κ = sgn(n 2 ) = ±1. (5.8)
∂z 2 ∂X 2 ∂Y 2
However, this equation has no general solution (not integrable). The nonlinearity creates an
effective potential κ|A|2 ∼ −V . For κ = +1(n 2 > 0) high intensities correspond to potential
dips.
Conserved quantities
For our further discussions we want to investigate the conserved quantities. We start with
48
5.1 Beam evolution in isotropic Kerr materials NLO
The second term corresponds to a potential energy since |A|4 can be interpreted as a multi-
plication of particle density ∼ |A|2 times Potential ∼ |A|2 . We can proof that the Hamiltonian
is conserved in a similar way as for the total power.
Second moment
A third quantity we want to introduce is called the second moment defined by
ˆ∞ ˆ∞
S(Z ) = dX dY (X 2 + Y 2 )|A(X , Y , Z )|2 . (5.11)
−∞ −∞
Despite the other quantities
q the second moment is not conserved. It is a measure of the
beam width 〈W (Z )〉 = S(Z )
Q .
We can use the definition of S(Z ) to derive its evolution equation. However, we only want to
mention the result:
d2
S = 4H conserved. (5.12)
dZ 2
In the following we want to distinguish two cases:
49
5.1 Beam evolution in isotropic Kerr materials NLO
ˆ∞ ˆ∞
¯ ∂A ¯2 ¯ ∂A ¯2
µ¯ ¯ ¯ ¯ ¶
4
a.) H = dX dY ¯
¯ ¯ + ¯ ¯ − κ|A| > 0.
∂X ¯ ¯ ∂Y ¯
−∞ −∞
We fulfill this case always for κ = −1(n 2 < 0). Here n 2 corresponds to a defocussing nonlin-
earity. The case κ = +1 also works for prevailing kinetic energy (e. g. curved phase fronts).
Then we find
d2
S = 4H > 0. (5.13)
dZ 2
The solutions of this equations are parabolas which diverge for Z → ∞
s
S(Z )
〈W (Z )〉 = → ∞. (5.14)
Q
S S
Z
Z
ZC
Fig. 7: Solutions of the second moment S(Z ) for H > 0 (left) and H < 0 (right). Since S(Z ) > 0 per
definition, the evolution path is not possible after S(Z ) = 0.
ˆ∞ ˆ∞
¯ ∂A ¯2 ¯ ∂A ¯2
µ¯ ¯ ¯ ¯ ¶
4
b.) H = dX dY ¯
¯ ¯ + ¯ ¯ − κ|A| < 0.
∂X ¯ ¯ ∂Y ¯
−∞ −∞
This case only works for κ = +1(n 2 > 0). Using the evolution equation of the second moment
we find that the beam will collapse within a finite propagation distance
The collapse is unphysical because of zero width (paraxial approximation fails) and infinite
intensities (power expansion of P NL cannot be stopped after the third order). However, the
evolution is often very well described up to the damage threshold of the material.
Furthermore we want to calculate the critical power for self-focusing for κ = 1. We assume
that the initial field has a fixed beam profile f (X , Y ) but a varying amplitude: A(X , Y , Z =
50
5.2 Pulse evolution in a two-level medium NLO
ˆ∞ ˆ∞
¯ ∂A ¯2 ¯ ∂A ¯2
µ¯ ¯ ¯ ¯ ¶
H= dX dY ¯
¯ ¯ +¯
¯ ¯ − |A| = |A 0 |2 h kin − |A 0 |4 h int
4
(5.16)
∂X ¯ ∂Y ¯
−∞ −∞
ˆ∞ ˆ∞ ˆ∞ ˆ∞
¯ ∂ f ¯2 ¯ ∂ f
µ¯ ¯ ¯ ¯2 ¶
1 1
dY | f |4 .
¯
with h kin = dX dY ¯
¯ ¯ +¯ ¯ and h int = dX
2 ∂X ¯ ¯ ∂Y ¯ 2
−∞ −∞ −∞ −∞
Both h kin and h int are only shape dependent. We can easily find the critical power as
h kin
H = |A 0 |2 h kin − |A 0 |4 h int < 0 for |A 0 |2 > |A crit | = . (5.17)
h int
Every beam collapses above a certain (profile dependent) power level. The beam profile
with the minimum power for collapse is called the Townes soliton. It is a metastable (quasi-
stationary) state between spreading and self-focusing. It is cylindrically symmetric, has no
analytical solution and thus must be numerically determined. It looks similar to a Gaussian
profile (except of a simple exponential decay in the tails). Numerically we find Q Townes as
ˆ∞ ˆ∞
Q Townes = dX dY |A Townes (X , Y , Z )|2 = 5.85043. (5.18)
−∞ −∞
Lastly we want to consider the pulse evolution in a two-level-system. First we have to make
some assumptions:
• The propagation is entirely in z-direction (no diffraction).
• The medium is isotropic: β0 = n(ω0 ) ωc0 .
• The two-level-system is not damped (pulse duration ≪ damping time).
• The electric field is linearly polarized.
• Group velocity dispersion (GVD) of the background medium can be neglected.
51
5.2 Pulse evolution in a two-level medium NLO
The evolution equation for the slowly varying envelope of a pulse is (see (3.36))
· β0 ∗ 2
∂ i ∂ E ∂
· ¸ ¸
1
i + a(t , z) = p β0 |2 ∂t 2
P̃ (z, t ) e−i(β0 z−ω0 t )
NL
∂z v g ∂t 8β0 c ε0 ω0
2 |E
· β0 ∗ 2 ³
1 E ∂ i(β0 z−ω0 t )
´¸
=p β0 |2 ∂t 2
TLS
P (z, t )e e−i(β0 z−ω0 t )
8β0 c ε0 ω0
2 |E
E β0 ∗
à !
1 2 ∂ ∂2 TLS
=p −ω0 − 2iω0 + 2 P (z, t )
8β0 c 2 ε0 ω0 |E β0 |2 ∂t ∂t
s
ω20 E β0 ∗ TLS
≈− P (z, t ). (5.23)
8ncε0 |E β0 |2
Now we use the evolution equations of the polarization for zero detuning ℏω0 = E a − E b
(equations (2.27) and (2.28))
d TLS Ne 2
| ψb |r̂ |ψa |2 E I
®
i P =
dt ℏ (5.24)
dI 1
= E ∗ P TLS + c.c.
dt 2iℏN
q
Again we introduce a normalization z = Z0 Z , t = T0 T, a(Z , T ) = Iˆ0 A(Z , T ) and
E β0
P TLS (z, t ) = −i Pˆ0 P (Z , T ). (5.25)
|E β0 |
Then equation (5.23) becomes
s
∂ i Z0 ∂ Z0 ω20 E β0 ∗ E β0
· ¸ · ¸
i + A(Z , T ) = − q −i β P̂ 0 P (Z , T ) . (5.26)
∂z v g T0 ∂T ˆ 8ncε0 |E β0 | |E 0 |
I0
Inserting all the expressions for P TLS (5.25) and E (5.22) into the evolution equations we
find
s
E β0 ® 2 E β0
" #
Ne 2
· ¸
d 2
q
i −i β Pˆ0 P (Z , T ) = T0 Iˆ0 | ψb |r̂ |ψa | A(Z , T ) I (5.27)
dT |E 0 | ℏ |E β0 | ε0 cn(ω0 )
q
#∗ ·
T0 Iˆ0 E β0
s
E β0
" ¸
dI 2
= A(Z , T ) ˆ
−i β P 0 P (Z , T ) + c.c.
dT 2iℏN |E β0 | ε0 cn(ω0 ) |E 0 |
(5.28)
52
5.2 Pulse evolution in a two-level medium NLO
We choose our scale parameters in such a way, that the constants vanish. This leads us to the
following conditions:
s
(5.26)Z0 Z0 ω20
⇒ = 1, q P̂ 0 = 1,
v g T0 8ncε0
Iˆ0
q
Iˆ0 Ne 2
s
(5.27) T0 2
| ψb |r̂ |ψa |2
®
⇒ = 1, (5.29)
P̂ 0 ℏ ε0 cn(ω0 )
q
Iˆ0
s
(5.28) T0 2
⇒ P̂ 0 = 1.
ℏN ε0 cn(ω0 )
Some algebraic rearrangements result in the following solutions for the scale parameters
with µ := | ψb |r̂ |ψa |:
®
s
1 2n c ℏ ε0 1
P̂ 0 = N eµ, T0 = , Iˆ0 = ℏω0 v g N . (5.30)
eµ v g ω0 N 4
∂
P (Z , T ) = A I (5.31a)
∂T
∂ 1¡
I (Z , T ) = − A ∗ P + AP ∗
¢
(5.31b)
∂T µ2
∂ ∂
¶
P (Z , T ) = + A(Z , T ). (5.31c)
∂Z ∂T
|P |2 + I 2 = 1 (5.32)
∂ ¡ 2 ∂P ∂P ∗
∂I
· ¸
1
|P | + I 2 = P ∗ = P ∗ AI + P A ∗ I + 2I − (A ∗ P + AP ∗ ) = 0.
¢
+P + 2I
∂T ∂T ∂T ∂T 2
In the following we restrict to real valued solutions for A, P, I which simplifies (5.31b) to
∂
I (Z , T ) = −A(Z , T )P (Z , T ). (5.33)
∂T
Since (5.32) looks similar to the trigonometric P YTHAGOREAN, we introduce a substitution
∂Θ(Z , T )
P (Z , T ) = sin Θ(Z , T ), I (Z , T ) = − cos Θ(Z , T ), A(Z , T ) = − . (5.34)
∂T
53
5.2 Pulse evolution in a two-level medium NLO
This inherently fulfils the conservation law |P |2 +I 2 = 1, but also equations (5.31a) and (5.31b).
An evaluation of equation (5.31c) leads us to
∂ ∂ ∂
µ ¶
+ Θ(Z , T ) = − sin Θ(Z , T ). (5.35)
∂Z ∂T ∂T
For further simplification we transform into a moving spatial frame Z̃ = 2Z −T which changes
the derivatives accordingly
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂2
µ ¶ µ 2 ¶
⇒2 , ⇒ − , + ⇒ + . (5.36)
∂Z ∂ Z̃ ∂T ∂T ∂ Z̃ ∂Z ∂T ∂T ∂T 2 ∂ Z̃ 2
∂2 ∂2
µ ¶
+ Θ( Z̃ , T ) = sin Θ( Z̃ , T ) (5.37)
∂ Z̃ 2 ∂T 2
Z̃ − vT
· µ ¶¸
Θ( Z̃ , T ) = 4 arctan exp ± p +δ , (5.38)
1 − v2
where 0 ≤ v < 1 is the velocity of the soliton. We observe that the pendula make a full turn
in angular space changing their phase from 0 to 2π (kink, the case of + sign) or from 2π to 0
(anti-kink, the case of − sign). For our description this means that the inversion goes from
−1 to +1 and returns to −1. The given solution is displayed in figure 8.
54
5.2 Pulse evolution in a two-level medium NLO
Fig. 8: Basic single soliton solution, simulated in Mathematica as systems of spring-coupled torsional
pendula, see V.G. Ivancevic and T.T. Ivancevig, JGSP 31(2013) 1–56
55