✴️

Nonlinear Optics: SHG, Phase Matching and Solitons

Solve the coupled-amplitude equations for second-harmonic generation with pump depletion, tune birefringent phase matching in BBO and KDP or quasi-phase matching in PPLN, then send a pulse through a Kerr fibre with a converged split-step Fourier solver to see self-phase modulation and optical solitons.

P4 specialist extension Prerequisites: EM waves in media, polarization and birefringence, dispersion and ultrashort pulses, phasor addition.

Type exact values in the number boxes • Hover or tap a plot for a cursor readout • Share the setup with “Copy link”
η — SH conversion
ΔkL — Phase mismatch (rad)
ΓL — Nonlinear drive
N — Soliton order
φNL — Peak SPM phase (rad)

Second-harmonic generation

Pump and SH intensity along the crystal normalised to I₁(0)

Solid: RK4 solution of the coupled equations, I₁/I₁(0) (pump) and I₂/I₁(0) (second harmonic). Dashed: closed-form reference for the same parameters: tanh²(Γz) when Δk = 0, otherwise the undepleted (Γz)² sinc²(Δkz/2). For a poled crystal Γ → (2/π)Γ and Δk → Δk − 2π/Λ. The upper panel rescales the SH axis by the stated power of ten when η is small.

Efficiency versus phase mismatch η = I₂(L)/I₁(0)

Solid: numerical η over a sweep of Δk·L at the current L and I₁(0). Dashed: undepleted (ΓL)² sinc²(ΔkL/2). Marker: the current crystal setting.

Phase-matching diagram indices vs θ

—

n(ω) / n(2ω)
—
Δk = k₂ − 2k₁
—
Coherence length Lc = π/|Δk|
—
Phase-matching angle θpm / θ
—
Angular acceptance (FWHM of sinc²)
—
Walk-off ρ of the SH (not modelled)
—
deff / coupling κ
—
Γ = κ√I₁(0); ΓL
—
η numerical (RK4)
—
η closed form
—
η with 2× steps (convergence)
—
Energy: max |I₁ + I₂ − I₁(0)|/I₁(0)
—
Manley–Rowe: N₁(0) → N₁(L) + 2N₂(L)
—
QPM vs ideal Δk = 0 (same d)
—
Cursor
—

Kerr effect: self-phase modulation and solitons

Pulse evolution |A(z, T)|² peak-normalised to P₀

Retarded time T = t − z/vg across, distance z along the fibre upwards. Colour: power |A|²/P₀ (inferno, linear). Computed with the split-step solver at the chosen step count.

Temporal profile |A|²/P₀

Dashed: input. Solid: output at z = Lf. Dotted: the analytic fundamental soliton √P₀ sech(T/T₀), shown when β₂ < 0.

Power spectrum input peak = 1

|Ã(ν)|² against the offset ν − ν₀ from the carrier. SPM creates new frequencies. Dispersion alone only rephases them, so the spectrum stays unchanged.

Step-count convergence log–log

Relative L2 error of the output field versus the number of split steps. The reference is the analytic soliton for an N = 1 sech input with β₂ < 0, and otherwise a run with 4× the finest step count. The dashed guide has slope −2. Vertical line: your step count.

Dispersion length LD = T₀²/|β₂|
—
Nonlinear length LNL = 1/(γP₀)
—
Soliton order N = √(LD/LNL)
—
Soliton period z₀ = πLD/2
—
Peak SPM phase γP₀Lf
—
Numerical phase at the peak (β₂ = 0)
—
Output peak / P₀
—
RMS width out/in (time)
—
RMS width out/in (spectrum)
—
Energy out/in
—
Error at this step count
—
Observed order of convergence
—
Grid
—
💡 How to Use

Learn with this tool

Learning objectives

  • Derive and solve the SHG coupled-amplitude equations, and explain why the undepleted sinc² law and the depleted tanh² law are two limits of one energy-conserving model.
  • Compute a phase-matching angle from Sellmeier data, and a quasi-phase-matching period and its (2/π)² efficiency penalty.
  • Use LD, LNL and N to predict whether a pulse in a Kerr fibre broadens spectrally, broadens temporally or propagates as a soliton, and show that a numerical solver has converged.

Prerequisites (P4 specialist topic)

Model

SHG. A monochromatic plane-wave pump at ω and its second harmonic at 2ω travel collinearly along z through a lossless crystal. The fields are Ej = Re{Aj(z) ei(kjz − ωjt)}, with Aj the peak amplitude. The tool integrates the intensity-normalised amplitudes aj = √(njε₀c/2) Aj, so |aj|² = Ij exactly:

da₁/dz = iκ s(z) a₂ a₁* eiΔkz,   da₂/dz = iκ s(z) a₁² e−iΔkz,   κ = ω deff √(2/(n₁²n₂ε₀c³)),   Δk = k₂ − 2k₁.

With the same κ in both equations, d(I₁ + I₂)/dz = 0 identically. Pump depletion therefore conserves energy by construction, and photon fluxes satisfy Manley–Rowe: N₁ + 2N₂ = const (two pump photons make one SH photon). s(z) = ±1 is the sign of d(z): always +1 in a bulk crystal, a square wave of period Λ in a poled one. The solver is fixed-step RK4 (core.rk4Step). A poled crystal is integrated domain by domain so that no step crosses a sign flip.

Kerr. A pulse envelope A(z, T) in watts½ obeys the nonlinear Schrödinger equation ∂A/∂z = −(iβ₂/2)∂²A/∂T² + iγ|A|²A, with T = t − z/vg. The symmetric split-step method applies half a dispersion step exactly in the Fourier domain, then the Kerr phase exp(iγ|A|²h) exactly, then the other half dispersion step. Its global error is O(h²).

I₁, I₂
pump and SH intensities (W/m²); η = I₂(L)/I₁(0)
deff
effective nonlinear coefficient (pm/V) for the chosen polarisations and angles
κ, Γ
coupling (m⁻¹(W/m²)^−½) and gain coefficient Γ = κ√I₁(0) (m⁻¹)
Δk, Lc
wave-vector mismatch 2ω(n₂ − n₁)/c; coherence length π/|Δk|
θ, θpm, ρ
angle between k and the optic axis; phase-matching angle; walk-off angle of the e-wave
Λ
poling period; first-order QPM needs Λ = 2π/Δk = 2Lc
β₂, γ
group-velocity dispersion (ps²/km) and Kerr coefficient γ = n₂ω₀/(cAeff) (W⁻¹km⁻¹)
T₀, P₀
pulse half-width (sech: TFWHM = 1.763 T₀; Gaussian: 1.665 T₀) and peak power
LD, LNL, N
T₀²/|β₂|, 1/(γP₀), soliton order √(LD/LNL)

Assumptions and validity: infinite plane waves (no diffraction, no walk-off), CW or quasi-CW pump for SHG (no group-velocity mismatch), no absorption, Kleinman symmetry, and the slowly varying envelope approximation (|d²A/dz²| ≪ k|dA/dz|). The NLSE keeps β₂ and instantaneous Kerr response only: no β₃, Raman or self-steepening. The fibre is lossless and single-mode, with Aeff absorbed into γ.

Derivation: coupled equations, sinc², tanh², Manley–Rowe and QPM

Driven wave equation. With PNL(t) = 2ε₀deffE(t)², the 2ω part of E² is ½A₁²e2i(k₁z−ωt) + c.c., so the peak polarisation phasor is P2ω = ε₀deffA₁²e2ik₁z. Inserting E₂ into ∇²E − (n²/c²)∂²E/∂t² = μ₀∂²PNL/∂t² and dropping d²A₂/dz² (slowly varying envelope) gives dA₂/dz = (iω₂/(2n₂cε₀))P2ωe−ik₂z = i(ωdeff/(n₂c))A₁²e−iΔkz. The ω part of E² couples A₂A₁*, which gives dA₁/dz = i(ωdeff/(n₁c))A₂A₁*eiΔkz.

Normalisation. Substituting Aj = aj√(2/(njε₀c)) gives the symmetric form above with κ = ωdeff√(2/(n₁²n₂ε₀c³)). Then d(|a₁|² + |a₂|²)/dz = 2Re[iκX] + 2Re[iκX*] = 0 with X = a₁*²a₂eiΔkz, which is energy conservation. Since Nj = Ij/(ħωj) and ω₂ = 2ω₁, this is the same statement as N₁ + 2N₂ = const (Manley–Rowe).

Undepleted limit. If a₁ ≈ a₁(0), then a₂(L) = iκa₁(0)² ∫₀L e−iΔkzdz, so |a₂|² = κ²I₁(0)²L² sinc²(ΔkL/2): η = (ΓL)² sinc²(ΔkL/2). This formula predicts η > 1 when ΓL > 1, which signals that its assumption has failed.

Depleted, Δk = 0. Take a₁ real and a₂ = iu. Then da₁/dz = −κua₁ and du/dz = κa₁², with a₁² + u² = I₁(0). This gives du/dz = κ(I₁(0) − u²), so u = √I₁(0) tanh(Γz) and η = tanh²(ΓL) (Armstrong, Bloembergen, Ducuing and Pershan, 1962).

Birefringent phase matching. In a negative uniaxial crystal (ne < no), an extraordinary SH wave at angle θ sees 1/ne(θ)² = cos²θ/no² + sin²θ/ne². Setting ne(2ω, θ) = no(ω) gives sin²θpm = [no(ω)⁻² − no(2ω)⁻²]/[ne(2ω)⁻² − no(2ω)⁻²]. Near θpm, Δk ≈ (∂Δk/∂θ)δθ, so the sinc² in ΔkL sets an angular acceptance δθFWHM = 5.566/(L|∂Δk/∂θ|).

Quasi-phase matching. Write s(z) = Σm odd (4/(πm)) sin(2πmz/Λ). The m = 1 term contains (2/π)ei2πz/Λ/i. When 2π/Λ = Δk it cancels e−iΔkz, so on average the crystal behaves like a phase-matched one with deff → (2/π)deff, and the efficiency is multiplied by (2/π)² ≈ 0.405. Within each domain the SH still grows along an arc, which is the staircase visible in the growth plot.

Solitons. In normalised units U = A/√P₀, ξ = z/LD and τ = T/T₀, the NLSE becomes iUξ = (sgn β₂/2)Uττ − N²|U|²U. For β₂ < 0 and N = 1, U = sech(τ)eiξ/2 is an exact solution. The chirp from dispersion is cancelled by the SPM chirp at every z.

Exercise 1 — Coherence length without phase matching

  1. In unpoled LiNbO₃ at 1064 nm, ne(ω) = 2.1555 and ne(2ω) = 2.2336. Predict Lc = λ/(4Δn) and the largest SH intensity reachable anywhere in the crystal, relative to a phase-matched crystal of length Lc.
  2. Run No poling (L = 0.05 mm).
  3. Read Lc and find the z of the first SH maximum and minimum with the cursor. Compare the peak value with η(Δk = 0) at z = Lc.
  4. Why does the SH return to zero, and where does its energy go?
Show answer

Δn = 0.0781, so Lc = 1.064/(4 × 0.0781) µm ≈ 3.41 µm. The SH peaks at z = Lc, 2π/Δk × ½, and vanishes at 2Lc. At the peak, sinc²(π/2) = (2/π)² ≈ 0.405 of the phase-matched value at the same length. Beyond Lc, SH generated earlier arrives out of phase with the local driving polarisation, and the energy flows back to the pump. I₁ + I₂ stays constant to about 10⁻¹³.

Exercise 2 — Limiting case: when does the undepleted formula fail?

  1. For BBO at θpm (Δk = 0), L = 10 mm, κ = 1.557 × 10⁻⁴ m⁻¹(W/m²)−½. Predict η at 1 MW/cm² and at 100 MW/cm² from (ΓL)² and from tanh²(ΓL).
  2. Run BBO, weak pump, then type 100 into the intensity box.
  3. Read η numerical, η closed form and the energy error. Note how far apart the dashed and solid curves are in the efficiency sweep.
  4. Show that tanh²(x) → x² as x → 0. What is the smallest ΓL at which the two formulas differ by 10 %?
Show answer

At 1 MW/cm² = 10¹⁰ W/m², Γ = 15.57 m⁻¹ and ΓL = 0.156. Then (ΓL)² = 0.0242 and tanh² = 0.0238: they agree to 1.6 %. At 100 MW/cm², ΓL = 1.557. The undepleted formula gives 2.42, an impossible 242 %, while tanh² = 0.837, matching the solver. tanh x = x − x³/3 + …, so tanh²x ≈ x²(1 − 2x²/3). A 10 % difference occurs at x ≈ 0.39 (η ≈ 14 %).

Exercise 3 — Quasi-phase matching costs (2/π)²

  1. PPLN, d33 = 25 pm/V, L = 1 mm, 1 MW/cm². Predict the first-order period and ηQPM/η(Δk = 0, same d).
  2. Run PPLN QPM. Then change Λ by +1 % and by +5 %.
  3. Read the QPM ratio readout and look at the staircase in the growth plot. Measure η at Λ + 1 %.
  4. Why is a 1 % period error serious in a 1 mm crystal but a 1 % error in L is not?
Show answer

Λ = 2π/Δk = 6.818 µm and the ratio is (2/π)² = 0.405. For this 1 mm crystal (147 periods) the tool gives 0.408; the ratio tends to 0.4053 for long crystals (0.4053 in the node tests at 400 periods). At Λ + 1 % the residual mismatch is Δk − 2π/Λ ≈ 0.0099Δk ≈ 9.1 mm⁻¹, so ΔkL ≈ 9.1 rad: past the first zero of the sinc² (2π), and η drops to a few percent of its matched value. A period error accumulates over all ≈ 150 periods, while a length error only changes where the growth stops.

Exercise 4 — Balancing dispersion with SPM: the fundamental soliton

  1. With β₂ = −20 ps²/km, γ = 1.3 W⁻¹km⁻¹ and T₀ = 1 ps, what P₀ gives N = 1? What happens after 5LD at that power and at P₀ = 0?
  2. Run Soliton N = 1, then Dispersion only, and finally raise P₀ by 50 %.
  3. Read the output peak/P₀, the rms width ratios and the error at this step count against the analytic soliton.
  4. Why does a slightly wrong P₀ still produce a soliton-like pulse after reshaping?
Show answer

P₀ = |β₂|/(γT₀²) = 20/1.3 = 15.4 W and LD = 50 m. The N = 1 pulse leaves 250 m unchanged (peak 1.000, both widths 1.00), with a field error of about 3 × 10⁻⁴ against the analytic soliton at 200 steps. Without γ (Dispersion only, 2LD) the rms width grows by √5 ≈ 2.24 (by ≈ 5.1 after 5LD) and the spectrum is unchanged. At 1.5 P₀ (N = 1.22) the pulse first narrows and oscillates, then sheds a little dispersive radiation and settles to a soliton whose amplitude matches its new width. Solitons are attractors of the integrable NLSE.

Exercise 5 — SPM only: counting spectral peaks and checking convergence

  1. A Gaussian with φmax = γP₀L = 4.5π and β₂ = 0. How many spectral peaks (M ≈ φmax/π + ½)? What rms broadening, √(1 + 4φmax²/(3√3))?
  2. Run SPM only, then set the split steps to 1.
  3. Count the peaks, then read the spectral rms ratio, the numerical phase at the peak and the error estimate.
  4. Why is one step already exact here, and why is it not when β₂ ≠ 0?
Show answer

M ≈ 5 peaks and an rms ratio of √(1 + 0.770 × 199.9) ≈ 12.4. With β₂ = 0 only the nonlinear operator acts and it is applied exactly: |A| never changes, so exp(iγ|A|²L) is the exact solution, and 1 step gives the same answer as 2000. When both operators act they do not commute. Symmetric splitting leaves an O(h³) error per step, which is O(h²) globally: the convergence plot has slope −2.

Worked example — frequency doubling a Nd:YAG laser in BBO

A 1064 nm beam of peak intensity 100 MW/cm² is doubled in a 10 mm BBO crystal cut for type I (o + o → e). Find θpm, the conversion efficiency, the angular tolerance, and whether walk-off matters for a 1 mm diameter beam.

  1. Eimerl Sellmeier: no(1064) = 1.6545, no(532) = 1.6742, ne(532) = 1.5547. sin²θ = (1.6545⁻² − 1.6742⁻²)/(1.5547⁻² − 1.6742⁻²) = 0.1502, so θpm = 22.8°.
  2. deff = d₃₁ sinθ + d₂₂ cosθ = 0.04 × 0.387 + 2.2 × 0.922 = 2.04 pm/V. κ = ωdeff√(2/(n₁²n₂ε₀c³)) = 1.557 × 10⁻⁴ m⁻¹(W/m²)−½.
  3. Γ = κ√(10¹² W/m²) = 155.7 m⁻¹ and ΓL = 1.557, so η = tanh²(1.557) = 0.837. The undepleted formula would give 2.42: it is invalid here.
  4. ∂Δk/∂θ = −1.09 × 10⁶ rad⁻¹m⁻¹, so δθFWHM = 5.566/(0.01 × 1.09 × 10⁶) = 0.51 mrad (internal angle). Tilting the crystal by 0.26 mrad halves the weak-pump efficiency.
  5. Walk-off ρ = 55.6 mrad (3.19°). The e-polarised SH drifts sideways by ρL = 0.56 mm over 10 mm, comparable to the 0.5 mm beam radius. The plane-wave η is an upper bound: a Boyd–Kleinman treatment would give a noticeably lower value. Check steps 1–4 with the BBO, depletion preset set to 100 MW/cm².

When the model fails

  • Finite beams: real beams diffract and the extraordinary SH walks off (ρ above). The plane-wave η overestimates focused-beam SHG. Use the Boyd–Kleinman focusing function h(ξ, B) for Gaussian beams. Spatial and temporal pump profiles also average η over I(r, t).
  • Pulsed pumps: group-velocity mismatch between ω and 2ω separates the pulses after LGVM = τ/|1/vg1 − 1/vg2|, which is millimetres for 100 fs pulses in BBO. The CW equations then overestimate η and miss spectral filtering.
  • High intensity: back-conversion when Δk ≠ 0, two-photon absorption, photorefractive damage (LiNbO₃ in the visible), thermal dephasing, and crystal damage thresholds (≈ 1–10 GW/cm² for ns pulses) are not modelled. Above those levels the slider shows physics the crystal would not survive.
  • Sellmeier data are valid only in their fit ranges and at room temperature. Phase-matching angles and QPM periods shift by ~0.1° or ~0.01 µm per few kelvin.
  • NLSE: femtosecond pulses or kilowatt powers need β₃, Raman scattering (soliton self-frequency shift, fission) and self-steepening. Fibre loss makes N decrease along the fibre. A finite time window with periodic boundaries wraps any radiation that reaches its edge; the tool warns you when more than 10⁻⁶ of the energy reaches the outer 5 % of the window.

References

  • J. A. Armstrong, N. Bloembergen, J. Ducuing and P. S. Pershan, “Interactions between light waves in a nonlinear dielectric,” Phys. Rev. 127, 1918–1939 (1962).
  • R. W. Boyd, Nonlinear Optics, 4th ed. (Academic Press, 2020), ch. 2 (coupled-wave equations, phase matching, QPM, Manley–Rowe).
  • D. Eimerl, L. Davis, S. Velsko, E. K. Graham and A. Zalkin, “Optical, mechanical, and thermal properties of barium borate,” J. Appl. Phys. 62, 1968–1983 (1987).
  • F. Zernike, “Refractive indices of ammonium dihydrogen phosphate and potassium dihydrogen phosphate between 2000 Å and 1.5 µ,” J. Opt. Soc. Am. 54, 1215–1220 (1964).
  • D. E. Zelmon, D. L. Small and D. Jundt, “Infrared corrected Sellmeier coefficients for congruently grown lithium niobate and 5 mol.% magnesium oxide–doped lithium niobate,” J. Opt. Soc. Am. B 14, 3319–3322 (1997).
  • D. N. Nikogosyan, Nonlinear Optical Crystals: A Complete Survey (Springer, 2005): d-coefficients, walk-off, acceptance.
  • M. M. Fejer, G. A. Magel, D. H. Jundt and R. L. Byer, “Quasi-phase-matched second harmonic generation: tuning and tolerances,” IEEE J. Quantum Electron. 28, 2631–2654 (1992).
  • G. D. Boyd and D. A. Kleinman, “Parametric interaction of focused Gaussian light beams,” J. Appl. Phys. 39, 3597–3639 (1968).
  • G. P. Agrawal, Nonlinear Fiber Optics, 6th ed. (Academic Press, 2019), ch. 2 (split-step Fourier), ch. 4 (SPM, eq. 4.1.13) and ch. 5 (solitons).
  • G. Strang, “On the construction and comparison of difference schemes,” SIAM J. Numer. Anal. 5, 506–517 (1968) (second-order symmetric splitting).