physicskit.optics #

Returns:

ndarray

physicskit.optics.fabry_perot_free_spectral_range(length, n=1.0, wavelength=None)[source]#

Free spectral range \(c/2nL\) in frequency, or \(\lambda^2/2nL\) in wavelength.

Parameters:

length (float) – Mirror spacing \(L\).
n (float) – Index between the mirrors.
wavelength (float | None) – If given, return the range in wavelength units near this wavelength; otherwise in units of \(c\) per length unit.

Return type:

Returns:

float

physicskit.optics.field_grid(shape, dx)[source]#

Centered physical-coordinate grid for an array of the given shape.

Parameters:

shape (tuple of int) – (Ny, Nx) array shape.
dx (float) – Grid spacing (same units as wavelength elsewhere in this module).

Returns:

X, Y (ndarray of shape (Ny, Nx)) – Cartesian coordinates, with the origin at index (Ny//2, Nx//2).

Examples

>>> X, Y = field_grid((4, 4), dx=1.0)
>>> X[0].tolist()
[-2.0, -1.0, 0.0, 1.0]

physicskit.optics.finesse(R)[source]#

Reflectance finesse \(\mathcal F = \pi\sqrt R/(1 - R)\): free spectral range over linewidth.

Parameters:: R (ArrayLike) – Mirror power reflectance.
Return type:: NDArray
Returns:: ndarray

Examples

>>> round(float(finesse(0.9)), 2)
29.8

physicskit.optics.flat_interface(n1, n2)[source]#

ABCD matrix for refraction at a flat interface, index n1 to n2.

\[\begin{split}M = \begin{pmatrix} 1 & 0 \\ 0 & n_1/n_2 \end{pmatrix}\end{split}\]

Parameters:

n1 (float) – Refractive index of the incident medium.
n2 (float) – Refractive index of the transmitted medium.

Returns:

ndarray of shape (2, 2)

Examples

>>> flat_interface(1.0, 1.5)
array([[1.        , 0.        ],
       [0.        , 0.66666667]])

physicskit.optics.fock_state(n, cutoff)[source]#

The number (Fock) state \(\lvert n\rangle\).

Parameters:

n (int) – Photon number.
cutoff (int) – Dimension of the truncated Fock space.

Returns:

ndarray of shape (cutoff,) – Complex state vector, all zeros except a 1 at index n.

Raises:

ValueError – If n is negative or n >= cutoff.

Examples

>>> fock_state(1, 4)
array([0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j])

physicskit.optics.fraunhofer_diffraction(aperture, wavelength, z, dx)[source]#

Far-field (Fraunhofer) diffraction pattern of aperture at distance z.

Uses the standard single-Fourier-transform Fraunhofer formula (Goodman convention),

\[U(x',y') = \frac{e^{ikz}}{i\lambda z} e^{i\frac{k}{2z}(x'^2+y'^2)} \iint U_0(x,y)\, e^{-i\frac{2\pi}{\lambda z}(x x' + y y')}\,dx\,dy,\]

with \(k = 2\pi/\lambda\). The double integral is exactly a 2D Fourier transform of aperture evaluated at spatial frequency \((f_x, f_y) = (x'/(\lambda z), y'/(\lambda z))\), computed here with a single FFT (the discrete sum is converted to a continuum integral by the dx**2 pixel-area factor). The output array lives on its own natural grid of spacing \(\lambda z / (N\,dx)\), not the input grid.

Parameters:

aperture (ndarray of shape (Ny, Nx)) – Input field (real or complex) immediately after the aperture.
wavelength (float) – Wavelength \(\lambda\), in the same length units as dx.
z (float) – Propagation distance to the observation plane.
dx (float) – Input grid spacing.

Returns:

ndarray of shape (Ny, Nx), complex – The far-field complex amplitude.

See also

fresnel_diffraction: The near-field approximation this is a further limit of.

Examples

>>> ap = circular_aperture((64, 64), dx=0.01, radius=0.05)
>>> U = fraunhofer_diffraction(ap, wavelength=0.5e-3, z=2.0, dx=0.01)
>>> U.shape
(64, 64)
>>> bool(intensity(U)[32, 32] == intensity(U).max())
True

physicskit.optics.free_space(d)[source]#

ABCD matrix for propagation through a distance d of free space (or any homogeneous medium).

\[\begin{split}M = \begin{pmatrix} 1 & d \\ 0 & 1 \end{pmatrix}\end{split}\]

A ray’s height changes in proportion to its angle and the distance traveled; its angle is unchanged.

Parameters:: d (float) – Propagation distance.
Returns:: ndarray of shape (2, 2)

Examples

>>> free_space(2.0)
array([[1., 2.],
       [0., 1.]])

physicskit.optics.fresnel_diffraction(aperture, wavelength, z, dx)[source]#

Near-field (Fresnel) diffraction pattern of aperture at distance z.

The single-Fourier-transform (“one-step”) Fresnel propagation method: multiply the input field by the Fresnel quadratic phase, Fourier transform, then multiply by the corresponding output quadratic phase and prefactor,

\[U(x',y') = \frac{e^{ikz}}{i\lambda z} e^{i\frac{k}{2z}(x'^2+y'^2)}\, \mathcal{F}\!\left[U_0(x,y)\, e^{i\frac{k}{2z}(x^2+y^2)}\right]_{f_x=x'/(\lambda z),\ f_y=y'/(\lambda z)},\]

the paraxial approximation to exact scalar diffraction. As with fraunhofer_diffraction(), the output lives on its own natural grid of spacing \(\lambda z/(N\,dx)\).

Parameters:

aperture (ndarray of shape (Ny, Nx)) – Input field (real or complex).
wavelength (float) – Wavelength \(\lambda\).
z (float) – Propagation distance.
dx (float) – Input grid spacing.

Returns:

ndarray of shape (Ny, Nx), complex – The propagated complex amplitude.

See also

angular_spectrum_propagate: The exact (non-paraxial) alternative, on the same input grid.

Examples

>>> ap = circular_aperture((64, 64), dx=0.01, radius=0.05)
>>> U = fresnel_diffraction(ap, wavelength=0.5e-3, z=5.0, dx=0.01)
>>> U.shape
(64, 64)

physicskit.optics.grin_medium(n0, n2_coeff, d)[source]#

ABCD matrix for a graded-index (GRIN) medium of length d.

For the standard quadratic radial index profile

\[n(r) = n_0\left(1 - \frac{n_2 r^2}{2}\right),\]

paraxial rays oscillate sinusoidally about the axis, giving

\[A = D = \cos(\sqrt{n_2}\,d), \qquad B = \frac{\sin(\sqrt{n_2}\,d)}{n_0\sqrt{n_2}}, \qquad C = -n_0\sqrt{n_2}\,\sin(\sqrt{n_2}\,d).\]

Input and output ray angles are measured outside the rod, in a medium of index 1: this is the inside-the-rod solution (\(A=D=\cos\), \(B=\sin/\sqrt{n_2}\), \(C=-\sqrt{n_2}\sin\)) sandwiched between the flat entrance and exit faces, flat_interface(n0, 1) @ M_inside @ flat_interface(1, n0), which is where the \(n_0\) factors come from. As n2_coeff \(\to 0\), \(A=D\to 1\), \(C\to 0\), and \(B \to d/n_0\) – the familiar reduced thickness of a homogeneous slab of index \(n_0\) in air, equal to free_space() only for \(n_0 = 1\).

Parameters:

n0 (float) – On-axis refractive index.
n2_coeff (float) – Quadratic index-gradient coefficient \(n_2\) (units of 1/length^2). Must be non-negative; 0 gives a homogeneous medium.
d (float) – Length of the GRIN medium.

Returns:

ndarray of shape (2, 2)

Examples

A tiny gradient barely distinguishable from a homogeneous medium of index 1 reduces to plain free-space propagation:

>>> np.allclose(grin_medium(n0=1.0, n2_coeff=1e-8, d=2.0), free_space(2.0), atol=1e-4)
True

physicskit.optics.hermite_gaussian_mode(x, y, z, beam, m, n)[source]#

Hermite-Gaussian \(\mathrm{TEM}_{mn}\) mode amplitude.

\[u_{mn}(x, y, z) = \frac{w_0}{w(z)}\, H_m\!\left(\frac{\sqrt2\,x}{w(z)}\right) H_n\!\left(\frac{\sqrt2\,y}{w(z)}\right) \exp\!\left[-\frac{x^2+y^2}{w(z)^2}\right] \exp\!\left[-\frac{ik(x^2+y^2)}{2R(z)}\right] \exp\!\left[i(m+n+1)\zeta(z)\right] \exp(-ikz)\]

with \(k = 2\pi/\lambda\). \(m=n=0\) reduces to the fundamental Gaussian mode carried by beam.

Parameters:

x (array_like) – Transverse coordinates (broadcastable).
y (array_like) – Transverse coordinates (broadcastable).
z (float) – Axial position.
beam (GaussianBeam) – The underlying fundamental-mode beam (sets \(w_0\), \(z_0\), \(\lambda\)).
m (int) – Transverse mode indices along \(x\) and \(y\).
n (int) – Transverse mode indices along \(x\) and \(y\).

Returns:

ndarray – Complex field amplitude, broadcast shape of x and y.

Examples

>>> import numpy as np
>>> beam = GaussianBeam(wavelength=1.0, w0=1.0, z0=0.0)
>>> u00 = hermite_gaussian_mode(0.0, 0.0, 0.0, beam, 0, 0)
>>> round(float(abs(u00)), 6)
1.0

physicskit.optics.intensity(U)[source]#

Optical intensity \(|U|^2\) of a complex field.

Parameters:: U (ndarray) – Complex (or real) field amplitude.
Returns:: ndarray – Real-valued intensity, same shape as U.

Examples

>>> intensity(np.array([3.0 + 4.0j]))
array([25.])

physicskit.optics.interactive_wigner_surface(W, x_grid, p_grid, title=None)[source]#

An interactive 3D Plotly surface plot of a Wigner quasi-probability distribution.

Parameters:

W (ndarray of shape (Nx, Np)) – Wigner function, e.g. from physicskit.optics.quantum_optics.compute_wigner_function().
x_grid (ndarray of shape (Nx,)) – Grid of x quadrature values.
p_grid (ndarray of shape (Np,)) – Grid of p quadrature values.
title (str, optional) – Plot title.

Returns:

plotly.graph_objects.Figure

Examples

>>> import numpy as np
>>> from physicskit.optics.quantum_optics import coherent_state, compute_wigner_function
>>> x = np.linspace(-4, 4, 40)
>>> p = np.linspace(-4, 4, 40)
>>> W = compute_wigner_function(coherent_state(1.0, 20), x, p)
>>> fig = interactive_wigner_surface(W, x, p)
>>> isinstance(fig, go.Figure)
True

physicskit.optics.laguerre_gaussian_mode(r, phi, z, beam, l, p)[source]#

Laguerre-Gaussian \(\mathrm{LG}_p^l\) mode amplitude.

\[u_{lp}(r, \phi, z) = \frac{w_0}{w(z)} \left(\frac{r\sqrt2}{w(z)}\right)^{|l|} L_p^{|l|}\!\left(\frac{2r^2}{w(z)^2}\right) \exp\!\left[-\frac{r^2}{w(z)^2}\right] \exp\!\left[-\frac{ikr^2}{2R(z)}\right] \exp(il\phi) \exp\!\left[i(|l|+2p+1)\zeta(z)\right] \exp(-ikz)\]

with \(k = 2\pi/\lambda\) and \(L_p^{|l|}\) the associated Laguerre polynomial.

Parameters:

r (array_like) – Polar transverse coordinates (broadcastable); r is the radial distance from the axis, phi the azimuthal angle.
phi (array_like) – Polar transverse coordinates (broadcastable); r is the radial distance from the axis, phi the azimuthal angle.
z (float) – Axial position.
beam (GaussianBeam) – The underlying fundamental-mode beam.
l (int) – Azimuthal (orbital) mode index.
p (int) – Radial mode index.

Returns:

ndarray – Complex field amplitude, broadcast shape of r and phi.

Examples

>>> import numpy as np
>>> beam = GaussianBeam(wavelength=1.0, w0=1.0, z0=0.0)
>>> u00 = laguerre_gaussian_mode(0.0, 0.0, 0.0, beam, 0, 0)
>>> round(float(abs(u00)), 6)
1.0

physicskit.optics.m2_beam_waist(z, wavelength, w0, M2=1.0)[source]#

Beam radius of a non-ideal (\(M^2 > 1\)) real laser beam, waist fixed at \(z=0\).

\[w(z) = w_0 \sqrt{1 + \left(\frac{z}{z_{R,\text{eff}}}\right)^2}, \qquad z_{R,\text{eff}} = \frac{\pi w_0^2}{M^2 \lambda}\]

A beam-quality factor \(M^2 \ge 1\) (equal to 1 for an ideal diffraction-limited Gaussian beam) reduces the effective Rayleigh range, so a real beam of the same waist diverges faster than the ideal case.

Parameters:

z (float or array_like) – Axial position(s), measured from the waist at \(z=0\).
wavelength (float) – Wavelength.
w0 (float) – Waist radius.
M2 (float, default=1.0) – Beam-quality factor, \(M^2 \ge 1\).

Returns:

float or ndarray – Beam radius w(z).

Examples

>>> round(float(m2_beam_waist(0.0, wavelength=0.5e-3, w0=0.1, M2=2.0)), 6)
0.1

physicskit.optics.multilayer_response(n_layers, d_layers, wavelength, n_incident=1.0, n_substrate=1.5, angle=0.0, polarization='s')[source]#

Reflection and transmission of a stack of thin layers by the characteristic-matrix method.

Parameters:

n_layers (ArrayLike) – Refractive index of each layer, from the incident side; complex values \(n + i\kappa\) with \(\kappa > 0\) describe absorption (time dependence \(e^{-i\omega t}\)).
d_layers (ArrayLike) – Physical thickness of each layer, in the units of wavelength.
wavelength (ArrayLike) – Vacuum wavelength(s).
n_incident (complex) – Index of the incident medium.
n_substrate (complex) – Index of the exit medium.
angle (float) – Angle of incidence in radians.
polarization (str)

Return type:

Returns:

dict – Complex amplitude coefficients "r", "t" and the power reflectance "R" and transmittance "T", each shaped like wavelength.

Examples

Bare glass reflects \(((1 - 1.5)/(1 + 1.5))^2 = 4\%\); a quarter-wave film of index \(\sqrt{1.5}\) removes it at the design wavelength:

>>> round(float(multilayer_response([], [], 550.0)["R"]), 4)
0.04
>>> n_ar = np.sqrt(1.5)
>>> float(multilayer_response([n_ar], [550.0 / (4 * n_ar)], 550.0)["R"]) < 1e-12
True

physicskit.optics.plot_beam_envelope(beam, z_range, n_points=200, ax=None)[source]#

Plot a Gaussian beam’s \(\pm w(z)\) waist envelope over a propagation range.

Parameters:

beam (physicskit.optics.gaussian.GaussianBeam) – The beam to plot.
z_range (tuple(float, float)) – (z_min, z_max) propagation-axis range to plot over.
n_points (int, default=200) – Number of sampled points.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> from physicskit.optics.gaussian import GaussianBeam
>>> beam = GaussianBeam(wavelength=0.5e-3, w0=0.1, z0=0.0)
>>> fig, ax = plot_beam_envelope(beam, z_range=(-10, 10))
>>> isinstance(fig, plt.Figure)
True
>>> bool(ax.lines[0].get_ydata().min() >= 0)
True

physicskit.optics.plot_diffraction_pattern(U, dx, ax=None, log_scale=False)[source]#

Plot the intensity pattern of a complex diffracted field as a 2D heatmap.

Parameters:

U (ndarray of shape (Ny, Nx)) – Complex field amplitude (e.g. from physicskit.optics.wave.angular_spectrum_propagate()).
dx (float) – Grid spacing, used to set physical axis extents.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.
log_scale (bool, default=False) – If True, plot log10(intensity + eps) to reveal faint diffraction fringes alongside the bright central peak.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> from physicskit.optics.wave import circular_aperture
>>> U = circular_aperture((64, 64), dx=1e-3, radius=5e-3).astype(complex)
>>> fig, ax = plot_diffraction_pattern(U, dx=1e-3)
>>> isinstance(fig, plt.Figure)
True

physicskit.optics.plot_ray_trace(system, y0, theta0, ax=None)[source]#

Plot a ray’s height through each element of an OpticalSystem.

Parameters:

system (physicskit.optics.ray.OpticalSystem) – The optical system to trace the ray through.
y0 (float) – Initial ray height.
theta0 (float) – Initial ray angle.
ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

Examples

>>> from physicskit.optics.ray import OpticalElement, OpticalSystem, free_space, thin_lens
>>> system = OpticalSystem([
...     OpticalElement(free_space(1.0), name="d1", length=1.0),
...     OpticalElement(thin_lens(1.0), name="lens"),
...     OpticalElement(free_space(1.0), name="d2", length=1.0),
... ])
>>> fig, ax = plot_ray_trace(system, y0=0.5, theta0=0.0)
>>> isinstance(fig, plt.Figure)
True
>>> len(ax.lines[0].get_xdata())
4

physicskit.optics.propagate_q(q, M)[source]#

Propagate a complex beam parameter through a paraxial optical system.

Applies the same bilinear transformation used for ray-transfer (ABCD) matrices,

\[q_{\text{out}} = \frac{A q + B}{C q + D},\]

with \(A, B, C, D\) the entries of M in the convention M = [[A, B], [C, D]].

Parameters:

q (complex) – Beam parameter before the system.
M (ndarray of shape (2, 2)) – Ray-transfer (ABCD) matrix of the optical system.

Returns:

complex – Beam parameter q_out after the system.

Examples

>>> import numpy as np
>>> free_space = np.array([[1.0, 2.0], [0.0, 1.0]])
>>> propagate_q(1j, free_space)
(2+1j)

physicskit.optics.q_to_beam_params(q, wavelength)[source]#

Extract the beam radius and radius of curvature from a complex beam parameter.

Inverts

\[\frac{1}{q} = \frac{1}{R} - i\,\frac{\lambda}{\pi w^2}.\]

At the beam waist \(\mathrm{Re}(1/q) = 0\), so R would formally diverge; this is handled explicitly and np.inf is returned there instead of raising or silently producing nan.

Parameters:

q (complex) – Complex beam parameter.
wavelength (float) – Wavelength (in the same length units as q).

Returns:

w (float) – Beam radius (\(1/e^2\) intensity radius).
R (float) – Radius of curvature of the wavefront; np.inf at the waist.

Examples

>>> import numpy as np
>>> zR = 3.0
>>> w, R = q_to_beam_params(1j * zR, wavelength=0.5e-3)
>>> round(float(w), 6), R
(0.021851, inf)

physicskit.optics.quarter_wave_band_gap(n_high, n_low)[source]#

Relative width \(\Delta\omega/\omega_0\) of the first gap of a quarter-wave stack.

\[\frac{\Delta\omega}{\omega_0} = \frac4\pi\arcsin\frac{n_H - n_L}{n_H + n_L}\]

(Yariv and Yeh, Optical Waves in Crystals, 1984, sec. 6.3).

Parameters:

n_high (float)
n_low (float)

Return type:

Returns:

float

Examples

>>> round(quarter_wave_band_gap(2.3, 1.38), 4)
0.3217

physicskit.optics.quarter_wave_stack(n_high, n_low, n_pairs, design_wavelength, cavity=False)[source]#

Layer indices and thicknesses of a quarter-wave (Bragg) mirror, or of a Fabry-Pérot cavity.

Parameters:

n_high (float) – Indices of the two materials; the stack starts with n_high.
n_low (float) – Indices of the two materials; the stack starts with n_high.
n_pairs (int) – Number of high/low pairs.
design_wavelength (float) – Wavelength at which every layer is a quarter wave thick.
cavity (bool) – If True, return two such mirrors separated by a half-wave (\(\lambda/2n_{\text{low}}\)) spacer of the low index: a Fabry-Pérot cavity with a resonance at design_wavelength.

Return type:

tuple[NDArray, NDArray]

Returns:

n_layers, d_layers (ndarray)

Examples

>>> n, d = quarter_wave_stack(2.3, 1.38, 2, 600.0)
>>> n.real.tolist(), (n.real * d).tolist()
([2.3, 1.38, 2.3, 1.38], [150.0, 150.0, 150.0, 150.0])

physicskit.optics.shg_coupled_amplitudes(z, delta_k=0.0, kappa=1.0, A0=1.0, qpm_period=None)[source]#

Integrate the SHG coupled-amplitude equations, with optional quasi-phase matching.

Parameters:

z (ArrayLike) – Increasing output positions, starting at 0.
delta_k (float) – Phase mismatch \(\Delta k\).
kappa (float) – Nonlinear coupling.
A0 (complex) – Fundamental amplitude at \(z = 0\); the harmonic starts at 0.
qpm_period (float | None) – If given, \(\kappa\) changes sign every half period (periodically poled crystal); first-order QPM uses \(2\pi/\Delta k\).

Return type:

Returns:

dict – Complex "A1" and "A2" at z.

Examples

>>> out = shg_coupled_amplitudes(np.linspace(0, 2, 5))
>>> bool(np.allclose(np.abs(out["A2"]) ** 2, np.tanh(np.linspace(0, 2, 5)) ** 2, atol=1e-6))
True

physicskit.optics.shg_phase_matched_efficiency(z, kappa=1.0, A0=1.0)[source]#

Conversion efficiency \(|A_2|^2/|A_0|^2 = \tanh^2(\kappa A_0 z)\) with pump depletion.

Parameters:

z (ArrayLike) – Propagation distance(s).
kappa (float)
A0 (float) – Input fundamental amplitude.

Return type:

Returns:

ndarray

Examples

>>> round(float(shg_phase_matched_efficiency(10.0)), 6)
1.0

physicskit.optics.shg_undepleted_power(L, delta_k, kappa=1.0, A1=1.0)[source]#

Second-harmonic power for an undepleted pump.

\[|A_2(L)|^2 = \kappa^2 |A_1|^4 L^2 \operatorname{sinc}^2(\Delta k L / 2) = \frac{4\kappa^2|A_1|^4}{\Delta k^2}\sin^2\frac{\Delta k L}{2}.\]

Parameters:

L (ArrayLike) – Crystal length(s).
delta_k (float) – Phase mismatch.
kappa (float) – Coupling constant.
A1 (float) – Pump amplitude.

Return type:

Returns:

ndarray

Examples

>>> round(float(shg_undepleted_power(np.pi, delta_k=1.0)), 12)  # one coherence length: 4/dk^2
4.0

physicskit.optics.single_slit_aperture(shape, dx, width)[source]#

A single slit of the given width along x, open along the full height in y.

Parameters:

shape (tuple of int) – (Ny, Nx) array shape.
dx (float) – Grid spacing.
width (float) – Full width of the slit opening (in the x direction).

Returns:

ndarray of shape (Ny, Nx) – Real-valued array of 1.0 (open) and 0.0 (blocked).

Examples

>>> ap = single_slit_aperture((3, 7), dx=1.0, width=2.5)
>>> ap[0].tolist()
[0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0]

physicskit.optics.spherical_mirror(R)[source]#

ABCD matrix for reflection from a spherical mirror of radius of curvature R.

\[\begin{split}M = \begin{pmatrix} 1 & 0 \\ -2/R & 1 \end{pmatrix}\end{split}\]

Parameters:: R (float) – Radius of curvature, positive for a concave mirror as seen by the incoming ray (i.e. a focusing mirror).
Returns:: ndarray of shape (2, 2)

Examples

>>> spherical_mirror(2.0)
array([[ 1.,  0.],
       [-1.,  1.]])

physicskit.optics.squeezed_state(xi, alpha=0.0, cutoff=30)[source]#

The displaced squeezed vacuum state, using the squeeze-then-displace convention.

Constructs

\[\lvert\xi, \alpha\rangle = \hat D(\alpha)\, \hat S(\xi)\, \lvert 0\rangle, \qquad \hat S(\xi) = \exp\!\left[\frac{\xi^* \hat a^2 - \xi \hat a^{\dagger 2}}{2}\right], \qquad \hat D(\alpha) = \exp\!\left(\alpha \hat a^\dagger - \alpha^* \hat a\right)\]

i.e. the vacuum is squeezed first and then displaced. (The opposite order, displace-then-squeeze, gives a different, also physically valid, state – squeeze-then-displace is the convention adopted here and is the more common one in the literature.) The ladder operators are built on the same truncated Fock basis, so this is only accurate while the squeezed/displaced amplitude is small relative to cutoff.

Parameters:

xi (complex) – Squeezing parameter \(\xi = r e^{i\theta}\).
alpha (complex, default=0.0) – Displacement amplitude.
cutoff (int, default=30) – Dimension of the truncated Fock space.

Returns:

ndarray of shape (cutoff,) – Complex, L2-normalized state vector.

Examples

>>> psi = squeezed_state(0.0, alpha=0.0, cutoff=6)
>>> psi
array([1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])

physicskit.optics.thick_lens(R1, R2, t, n, n_ext=1.0)[source]#

ABCD matrix for a thick lens: two curved interfaces separated by thickness t.

Composes, in the order light passes through them, refraction into the lens at the first surface, propagation across the lens body, and refraction back out at the second surface:

\[\begin{split}M = M_{R_2}\,M_t\,M_{R_1}, \qquad M_{R_1} = \text{curved\_interface}(R_1, n_{\text{ext}}, n), \\ M_t = \text{free\_space}(t) \ \text{(in medium } n\text{)}, \qquad M_{R_2} = \text{curved\_interface}(R_2, n, n_{\text{ext}}).\end{split}\]

Parameters:

R1 (float) – Radius of curvature of the first (entrance) surface.
R2 (float) – Radius of curvature of the second (exit) surface.
t (float) – Center thickness of the lens.
n (float) – Refractive index of the lens material.
n_ext (float, default=1.0) – Refractive index of the surrounding medium (air, by default).

Returns:

ndarray of shape (2, 2)

Examples

A symmetric biconvex lens in air; the resulting matrix has unit determinant, as any lossless ABCD system must:

>>> M = thick_lens(R1=0.1, R2=-0.1, t=0.01, n=1.5)
>>> round(float(np.linalg.det(M)), 8)
1.0

physicskit.optics.thin_lens(f)[source]#

ABCD matrix for a thin lens of focal length f.

\[\begin{split}M = \begin{pmatrix} 1 & 0 \\ -1/f & 1 \end{pmatrix}\end{split}\]

Parameters:: f (float) – Focal length (positive for a converging lens, negative for diverging).
Returns:: ndarray of shape (2, 2)

Examples

>>> thin_lens(0.5)
array([[ 1.,  0.],
       [-2.,  1.]])

physicskit.optics.type_i_phase_matching_angle(n_o_w, n_o_2w, n_e_2w)[source]#

Type-I (ooe) phase-matching angle of a negative uniaxial crystal.

Solves \(n_e^{2\omega}(\theta) = n_o^{\omega}\):

\[\sin^2\theta_m = \frac{(n_o^{\omega})^{-2} - (n_o^{2\omega})^{-2}} {(n_e^{2\omega})^{-2} - (n_o^{2\omega})^{-2}}.\]

Parameters:

n_o_w (float) – Ordinary index at the fundamental.
n_o_2w (float) – Ordinary and extraordinary indices at the harmonic.
n_e_2w (float) – Ordinary and extraordinary indices at the harmonic.

Return type:

Returns:

float – \(\theta_m\) in radians.

Raises:

ValueError – If the birefringence cannot compensate the dispersion.

Examples

KDP doubling 1064 nm (Sellmeier indices from Zernike, JOSA 54, 1215 (1964)):

>>> round(float(np.degrees(type_i_phase_matching_angle(1.4942, 1.5129, 1.4709))), 1)
41.3

physicskit.optics.wigner_negativity(W, x_grid, p_grid)[source]#

Wigner negativity: the integrated negative volume of a Wigner function.

\[N_W = \iint \max(0, -W(x,p))\; dx\, dp\]

A standard non-classicality diagnostic: zero for classical/Gaussian states (coherent, thermal, squeezed vacuum), strictly positive for genuinely non-classical states such as Fock states with \(n \ge 1\).

Parameters:

W (ndarray of shape (Nx, Np)) – Wigner function, e.g. from compute_wigner_function().
x_grid (ndarray of shape (Nx,)) – Grid of \(x\) quadrature values.
p_grid (ndarray of shape (Np,)) – Grid of \(p\) quadrature values.

Returns:

float – The (non-negative) integrated negative volume.

Examples

>>> import numpy as np
>>> vac = fock_state(0, 8)
>>> xg = np.linspace(-5, 5, 61)
>>> pg = np.linspace(-5, 5, 61)
>>> W = compute_wigner_function(vac, xg, pg)
>>> wigner_negativity(W, xg, pg)
0.0

Geometric ray optics: paraxial ray transfer (ABCD) matrices and optical systems.

In the paraxial approximation, a light ray at a given plane along the optical axis is fully described by two numbers – its height above the axis, \(y\), and the angle it makes with the axis, \(\theta\) (in radians, small-angle/paraxial regime) – collected into a state vector \((y, \theta)^T\). Every simple optical element (propagation through free space, refraction at an interface, a thin or thick lens, a curved mirror, a graded-index medium, …) acts on this vector as a linear map, a \(2\times 2\) “ABCD” matrix:

\[\begin{split}\begin{pmatrix} y_{\text{out}} \\ \theta_{\text{out}} \end{pmatrix} = \begin{pmatrix} A & B \\ C & D \end{pmatrix} \begin{pmatrix} y_{\text{in}} \\ \theta_{\text{in}} \end{pmatrix}.\end{split}\]

This formalism – developed piecemeal through the 19th century (Gauss’s theory of optical systems) and formalized for laser resonator design by Kogelnik and Li in the 1960s – reduces the analysis of an arbitrarily long chain of lenses, mirrors, and gaps to ordinary matrix multiplication: the matrix of a compound system is just the product of the matrices of its elements, applied in the order light encounters them (rightmost first). The same matrices reappear in physicskit.optics.gaussian to propagate the complex beam parameter of a Gaussian laser beam, and their trace controls whether a laser cavity is a stable resonator.

class physicskit.optics.ray.OpticalElement(matrix, name='', length=0.0)[source]#

Bases: object

A single named optical element wrapping one ABCD matrix.

Parameters:

matrix (array_like, shape (2, 2)) – The element’s ray transfer matrix, e.g. from thin_lens() or free_space().
name (str, default="") – Human-readable label (e.g. "f=50mm lens").
length (float, default=0.0) – Physical length occupied by this element along the optical axis (zero for a “thin” element such as a lens or mirror).

Examples

>>> elem = OpticalElement(thin_lens(0.05), name="focusing lens", length=0.0)
>>> elem.name
'focusing lens'

class physicskit.optics.ray.OpticalSystem(elements)[source]#

Bases: object

An ordered sequence of OpticalElement forming a compound optical system.

Parameters:: elements (list of OpticalElement) – The elements in the order light passes through them: elements[0] is hit first.

Examples

A single thin lens followed by propagation over its focal length focuses any parallel ray bundle back to the axis:

>>> f = 0.1
>>> sys = OpticalSystem([
...     OpticalElement(thin_lens(f), name="lens"),
...     OpticalElement(free_space(f), name="propagate to focus"),
... ])
>>> trajectory = sys.trace_ray(y0=0.01, theta0=0.0)
>>> abs(float(trajectory[-1, 0])) < 1e-12
True

is_stable()[source]#

Whether the system satisfies the resonator stability condition \(|A+D| \le 2\).

Returns:: bool

property stability_parameter#

Resonator stability parameter \((A+D)/2\) of the system matrix.

Returns:: float

system_matrix()[source]#

Total ABCD matrix of the system, \(M = M_n \cdots M_2 M_1\).

elements[0] is applied first (it is the rightmost factor), so it acts on the incoming ray state before any later element.

Returns:: ndarray of shape (2, 2)

trace_ray(y0, theta0)[source]#

Trace a single ray through every element, recording its state at each step.

Parameters:

y0 (float) – Initial height.
theta0 (float) – Initial angle, in radians.

Returns:

ndarray of shape (n_elements + 1, 2) – Row 0 is the input state [y0, theta0]; row i (for i >= 1) is the state after passing through elements[0], ..., elements[i-1].

physicskit.optics.ray.cavity_round_trip_matrix(elements)[source]#

Round-trip ABCD matrix of a laser cavity, given its elements in traversal order.

Equivalent to OpticalSystem(elements).system_matrix(); provided as a standalone function for cavities analyzed without constructing a full OpticalSystem.

Parameters:: elements (list of OpticalElement) – The elements encountered over one full round trip, in order.
Returns:: ndarray of shape (2, 2)

See also

cavity_stability: Test the resulting matrix for resonator stability.

Examples

>>> M = cavity_round_trip_matrix([
...     OpticalElement(spherical_mirror(2.0), name="M1"),
...     OpticalElement(free_space(1.0), name="gap"),
...     OpticalElement(spherical_mirror(2.0), name="M2"),
...     OpticalElement(free_space(1.0), name="gap"),
... ])
>>> cavity_stability(M)
True

physicskit.optics.ray.cavity_stability(M)[source]#

Resonator stability test \(|A+D| \le 2\) (equivalently \(|\operatorname{tr} M| \le 2\)).

A laser cavity with round-trip matrix M supports stable, non-diverging paraxial ray bundles if and only if this holds.

Parameters:: M (array_like, shape (2, 2)) – A round-trip ABCD matrix, e.g. from cavity_round_trip_matrix() or OpticalSystem.system_matrix().
Returns:: bool

Examples

>>> cavity_stability(np.array([[1.0, 0.0], [0.0, 1.0]]))
True
>>> cavity_stability(np.array([[3.0, 0.0], [0.0, 3.0]]))
False

physicskit.optics.ray.curved_interface(R, n1, n2)[source]#

ABCD matrix for refraction at a spherical interface of radius R, index n1 to n2.

\[\begin{split}M = \begin{pmatrix} 1 & 0 \\ (n_1-n_2)/(R n_2) & n_1/n_2 \end{pmatrix}\end{split}\]

Parameters:

R (float) – Radius of curvature of the interface, positive if the center of curvature lies on the outgoing (transmitted) side of the surface.
n1 (float) – Refractive index of the incident medium.
n2 (float) – Refractive index of the transmitted medium.

Returns:

ndarray of shape (2, 2)

See also

flat_interface: The \(R \to \infty\) limit.

Examples

>>> np.allclose(curved_interface(1.0e12, 1.0, 1.5), flat_interface(1.0, 1.5), atol=1e-6)
True

physicskit.optics.ray.flat_interface(n1, n2)[source]#

ABCD matrix for refraction at a flat interface, index n1 to n2.

\[\begin{split}M = \begin{pmatrix} 1 & 0 \\ 0 & n_1/n_2 \end{pmatrix}\end{split}\]

Parameters:

n1 (float) – Refractive index of the incident medium.
n2 (float) – Refractive index of the transmitted medium.

Returns:

ndarray of shape (2, 2)

Examples

>>> flat_interface(1.0, 1.5)
array([[1.        , 0.        ],
       [0.        , 0.66666667]])

physicskit.optics.ray.free_space(d)[source]#

ABCD matrix for propagation through a distance d of free space (or any homogeneous medium).

\[\begin{split}M = \begin{pmatrix} 1 & d \\ 0 & 1 \end{pmatrix}\end{split}\]

A ray’s height changes in proportion to its angle and the distance traveled; its angle is unchanged.

Parameters:: d (float) – Propagation distance.
Returns:: ndarray of shape (2, 2)

Examples

>>> free_space(2.0)
array([[1., 2.],
       [0., 1.]])

physicskit.optics.ray.grin_medium(n0, n2_coeff, d)[source]#

ABCD matrix for a graded-index (GRIN) medium of length d.

For the standard quadratic radial index profile

\[n(r) = n_0\left(1 - \frac{n_2 r^2}{2}\right),\]

paraxial rays oscillate sinusoidally about the axis, giving

\[A = D = \cos(\sqrt{n_2}\,d), \qquad B = \frac{\sin(\sqrt{n_2}\,d)}{n_0\sqrt{n_2}}, \qquad C = -n_0\sqrt{n_2}\,\sin(\sqrt{n_2}\,d).\]

Input and output ray angles are measured outside the rod, in a medium of index 1: this is the inside-the-rod solution (\(A=D=\cos\), \(B=\sin/\sqrt{n_2}\), \(C=-\sqrt{n_2}\sin\)) sandwiched between the flat entrance and exit faces, flat_interface(n0, 1) @ M_inside @ flat_interface(1, n0), which is where the \(n_0\) factors come from. As n2_coeff \(\to 0\), \(A=D\to 1\), \(C\to 0\), and \(B \to d/n_0\) – the familiar reduced thickness of a homogeneous slab of index \(n_0\) in air, equal to free_space() only for \(n_0 = 1\).

Parameters:

n0 (float) – On-axis refractive index.
n2_coeff (float) – Quadratic index-gradient coefficient \(n_2\) (units of 1/length^2). Must be non-negative; 0 gives a homogeneous medium.
d (float) – Length of the GRIN medium.

Returns:

ndarray of shape (2, 2)

Examples

A tiny gradient barely distinguishable from a homogeneous medium of index 1 reduces to plain free-space propagation:

>>> np.allclose(grin_medium(n0=1.0, n2_coeff=1e-8, d=2.0), free_space(2.0), atol=1e-4)
True

physicskit.optics.ray.spherical_mirror(R)[source]#

ABCD matrix for reflection from a spherical mirror of radius of curvature R.

\[\begin{split}M = \begin{pmatrix} 1 & 0 \\ -2/R & 1 \end{pmatrix}\end{split}\]

Parameters:: R (float) – Radius of curvature, positive for a concave mirror as seen by the incoming ray (i.e. a focusing mirror).
Returns:: ndarray of shape (2, 2)

Examples

>>> spherical_mirror(2.0)
array([[ 1.,  0.],
       [-1.,  1.]])

physicskit.optics.ray.thick_lens(R1, R2, t, n, n_ext=1.0)[source]#

ABCD matrix for a thick lens: two curved interfaces separated by thickness t.

Composes, in the order light passes through them, refraction into the lens at the first surface, propagation across the lens body, and refraction back out at the second surface:

\[\begin{split}M = M_{R_2}\,M_t\,M_{R_1}, \qquad M_{R_1} = \text{curved\_interface}(R_1, n_{\text{ext}}, n), \\ M_t = \text{free\_space}(t) \ \text{(in medium } n\text{)}, \qquad M_{R_2} = \text{curved\_interface}(R_2, n, n_{\text{ext}}).\end{split}\]

Parameters:

R1 (float) – Radius of curvature of the first (entrance) surface.
R2 (float) – Radius of curvature of the second (exit) surface.
t (float) – Center thickness of the lens.
n (float) – Refractive index of the lens material.
n_ext (float, default=1.0) – Refractive index of the surrounding medium (air, by default).

Returns:

ndarray of shape (2, 2)

Examples

A symmetric biconvex lens in air; the resulting matrix has unit determinant, as any lossless ABCD system must:

>>> M = thick_lens(R1=0.1, R2=-0.1, t=0.01, n=1.5)
>>> round(float(np.linalg.det(M)), 8)
1.0

physicskit.optics.ray.thin_lens(f)[source]#

ABCD matrix for a thin lens of focal length f.

\[\begin{split}M = \begin{pmatrix} 1 & 0 \\ -1/f & 1 \end{pmatrix}\end{split}\]

Parameters:: f (float) – Focal length (positive for a converging lens, negative for diverging).
Returns:: ndarray of shape (2, 2)

Examples

>>> thin_lens(0.5)
array([[ 1.,  0.],
       [-2.,  1.]])

Scalar wave optics: FFT-based Fresnel, Fraunhofer, and angular-spectrum diffraction.

Treats light as a scalar complex field \(U(x,y)\) obeying the Helmholtz equation, and propagates it between parallel planes using the three workhorse methods of Fourier optics (Goodman, Introduction to Fourier Optics):

angular_spectrum_propagate() – the exact scalar-diffraction solution, decomposing the field into plane waves (its 2D Fourier transform), advancing each by its own propagation phase \(e^{ik_z z}\), and re-synthesizing. Valid at any distance, including deep into the near field, as long as the grid resolves the field’s spatial frequencies.
fresnel_diffraction() – the paraxial (parabolic-wave) near-field approximation, computable as a single Fourier transform of the aperture times a quadratic phase.
fraunhofer_diffraction() – the further far-field approximation valid once \(z\) is large enough that the quadratic phase across the aperture itself is negligible; the diffraction pattern becomes simply the (scaled, phase-prefactored) Fourier transform of the aperture, the basis for classic single- and double-slit interference patterns.

All propagation routines expect a field sampled on a uniform square-pixel grid of spacing dx (same length units as wavelength), and treat array index (Ny//2, Nx//2) as the on-axis origin.

physicskit.optics.wave.angular_spectrum_propagate(U0, wavelength, z, dx)[source]#

Exact scalar diffraction propagation via the angular spectrum method.

Decomposes the input field into plane-wave components (its 2D Fourier transform), advances each by its own longitudinal propagation phase, and re-synthesizes:

\[\begin{split}U(x,y;z) = \mathcal{F}^{-1}\!\left[ \mathcal{F}[U_0](f_x,f_y)\; e^{ik_z z} \right], \qquad k_z = \begin{cases} \sqrt{k^2 - k_x^2 - k_y^2} & k_x^2+k_y^2 \le k^2 \quad \text{(propagating)}\\ i\sqrt{k_x^2+k_y^2 - k^2} & k_x^2+k_y^2 > k^2 \quad \text{(evanescent)} \end{cases}\end{split}\]

with \(k=2\pi/\lambda\), \(k_x = 2\pi f_x\), \(k_y = 2\pi f_y\). Evanescent orders get a purely imaginary \(k_z\), so \(e^{ik_z z}\) decays exponentially rather than producing nan/inf. Unlike fresnel_diffraction() and fraunhofer_diffraction(), this is not a paraxial approximation and the output remains on the same grid (shape and spacing dx) as the input.

Parameters:

U0 (ndarray of shape (Ny, Nx)) – Input complex (or real) field.
wavelength (float) – Wavelength \(\lambda\).
z (float) – Propagation distance.
dx (float) – Grid spacing (both input and output).

Returns:

ndarray of shape (Ny, Nx), complex – The propagated field, on the same grid as the input.

Examples

>>> ap = circular_aperture((64, 64), dx=0.01, radius=0.05)
>>> U = angular_spectrum_propagate(ap, wavelength=0.5e-3, z=0.05, dx=0.01)
>>> U.shape
(64, 64)
>>> bool(np.isclose(intensity(U).sum(), intensity(ap).sum(), rtol=0.02))
True

physicskit.optics.wave.circular_aperture(shape, dx, radius)[source]#

A circular (disk) aperture, transmittance 1 inside, 0 outside.

Parameters:

shape (tuple of int) – (Ny, Nx) array shape.
dx (float) – Grid spacing.
radius (float) – Aperture radius.

Returns:

ndarray of shape (Ny, Nx) – Real-valued array of 1.0 (open) and 0.0 (blocked).

Examples

>>> ap = circular_aperture((5, 5), dx=1.0, radius=1.5)
>>> float(ap[2, 2])
1.0
>>> float(ap[0, 0])
0.0

physicskit.optics.wave.double_slit_aperture(shape, dx, width, separation)[source]#

Two parallel slits of the given width, centers separated by separation.

Parameters:

shape (tuple of int) – (Ny, Nx) array shape.
dx (float) – Grid spacing.
width (float) – Full width of each slit opening.
separation (float) – Center-to-center distance between the two slits.

Returns:

ndarray of shape (Ny, Nx) – Real-valued array of 1.0 (open) and 0.0 (blocked).

Examples

>>> ap = double_slit_aperture((3, 11), dx=1.0, width=1.5, separation=6.0)
>>> ap[0].tolist()
[0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0]

physicskit.optics.wave.field_grid(shape, dx)[source]#

Centered physical-coordinate grid for an array of the given shape.

Parameters:

shape (tuple of int) – (Ny, Nx) array shape.
dx (float) – Grid spacing (same units as wavelength elsewhere in this module).

Returns:

X, Y (ndarray of shape (Ny, Nx)) – Cartesian coordinates, with the origin at index (Ny//2, Nx//2).

Examples

>>> X, Y = field_grid((4, 4), dx=1.0)
>>> X[0].tolist()
[-2.0, -1.0, 0.0, 1.0]

physicskit.optics.wave.fraunhofer_diffraction(aperture, wavelength, z, dx)[source]#

Far-field (Fraunhofer) diffraction pattern of aperture at distance z.

Uses the standard single-Fourier-transform Fraunhofer formula (Goodman convention),

\[U(x',y') = \frac{e^{ikz}}{i\lambda z} e^{i\frac{k}{2z}(x'^2+y'^2)} \iint U_0(x,y)\, e^{-i\frac{2\pi}{\lambda z}(x x' + y y')}\,dx\,dy,\]

with \(k = 2\pi/\lambda\). The double integral is exactly a 2D Fourier transform of aperture evaluated at spatial frequency \((f_x, f_y) = (x'/(\lambda z), y'/(\lambda z))\), computed here with a single FFT (the discrete sum is converted to a continuum integral by the dx**2 pixel-area factor). The output array lives on its own natural grid of spacing \(\lambda z / (N\,dx)\), not the input grid.

Parameters:

aperture (ndarray of shape (Ny, Nx)) – Input field (real or complex) immediately after the aperture.
wavelength (float) – Wavelength \(\lambda\), in the same length units as dx.
z (float) – Propagation distance to the observation plane.
dx (float) – Input grid spacing.

Returns:

ndarray of shape (Ny, Nx), complex – The far-field complex amplitude.

See also

fresnel_diffraction: The near-field approximation this is a further limit of.

Examples

>>> ap = circular_aperture((64, 64), dx=0.01, radius=0.05)
>>> U = fraunhofer_diffraction(ap, wavelength=0.5e-3, z=2.0, dx=0.01)
>>> U.shape
(64, 64)
>>> bool(intensity(U)[32, 32] == intensity(U).max())
True

physicskit.optics.wave.fresnel_diffraction(aperture, wavelength, z, dx)[source]#

Near-field (Fresnel) diffraction pattern of aperture at distance z.

The single-Fourier-transform (“one-step”) Fresnel propagation method: multiply the input field by the Fresnel quadratic phase, Fourier transform, then multiply by the corresponding output quadratic phase and prefactor,

\[U(x',y') = \frac{e^{ikz}}{i\lambda z} e^{i\frac{k}{2z}(x'^2+y'^2)}\, \mathcal{F}\!\left[U_0(x,y)\, e^{i\frac{k}{2z}(x^2+y^2)}\right]_{f_x=x'/(\lambda z),\ f_y=y'/(\lambda z)},\]

the paraxial approximation to exact scalar diffraction. As with fraunhofer_diffraction(), the output lives on its own natural grid of spacing \(\lambda z/(N\,dx)\).

Parameters:

aperture (ndarray of shape (Ny, Nx)) – Input field (real or complex).
wavelength (float) – Wavelength \(\lambda\).
z (float) – Propagation distance.
dx (float) – Input grid spacing.

Returns:

ndarray of shape (Ny, Nx), complex – The propagated complex amplitude.

See also

angular_spectrum_propagate: The exact (non-paraxial) alternative, on the same input grid.

Examples

>>> ap = circular_aperture((64, 64), dx=0.01, radius=0.05)
>>> U = fresnel_diffraction(ap, wavelength=0.5e-3, z=5.0, dx=0.01)
>>> U.shape
(64, 64)

physicskit.optics.wave.intensity(U)[source]#

Optical intensity \(|U|^2\) of a complex field.

Parameters:: U (ndarray) – Complex (or real) field amplitude.
Returns:: ndarray – Real-valued intensity, same shape as U.

Examples

>>> intensity(np.array([3.0 + 4.0j]))
array([25.])

physicskit.optics.wave.single_slit_aperture(shape, dx, width)[source]#

A single slit of the given width along x, open along the full height in y.

Parameters:

shape (tuple of int) – (Ny, Nx) array shape.
dx (float) – Grid spacing.
width (float) – Full width of the slit opening (in the x direction).

Returns:

ndarray of shape (Ny, Nx) – Real-valued array of 1.0 (open) and 0.0 (blocked).

Examples

>>> ap = single_slit_aperture((3, 7), dx=1.0, width=2.5)
>>> ap[0].tolist()
[0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0]

Gaussian beam propagation: the complex beam parameter, ABCD transformation, and higher-order modes.

A paraxial, monochromatic beam with a Gaussian transverse intensity profile is completely described, at any position \(z\) along its propagation axis, by a single complex number – the beam parameter \(q(z)\), introduced by Kogelnik & Li (Appl. Opt. 5, 1550, 1966) – through

\[\frac{1}{q(z)} = \frac{1}{R(z)} - i\,\frac{\lambda}{\pi w(z)^2},\]

where \(w(z)\) is the \(1/e^2\) intensity radius and \(R(z)\) the radius of curvature of the wavefronts. Propagation through any paraxial optical system described by a ray-transfer (ABCD) matrix updates \(q\) by the same bilinear (Mobius) transformation that acts on ray-tracing matrices, which is what makes the complex-\(q\) formalism so useful: free-space propagation, lenses, mirrors, and interfaces all become 2x2 matrix multiplications, exactly as in ordinary ray optics. This module also covers the free-space eigenmodes of the paraxial wave equation built on top of the fundamental Gaussian – the Hermite-Gaussian and Laguerre-Gaussian mode families – and the empirical \(M^2\) beam-quality factor used to describe real, non-diffraction-limited laser beams.

The ABCD matrix convention used throughout (shared with physicskit.optics.ray) is M = [[A, B], [C, D]] acting on a ray state [y, theta] as state_out = M @ state_in; e.g. free space of length d is [[1, d], [0, 1]] and a thin lens of focal length f is [[1, 0], [-1/f, 1]].

class physicskit.optics.gaussian.GaussianBeam(wavelength, w0, z0=0.0)[source]#

Bases: object

A fundamental (\(\mathrm{TEM}_{00}\)) Gaussian laser beam.

Parameters:

wavelength (float) – Wavelength, in the same length units as w0 and z0.
w0 (float) – Waist radius (\(1/e^2\) intensity radius at the narrowest point).
z0 (float, default=0.0) – Axial position of the waist.

property divergence_angle#: Far-field half-angle divergence, \(\theta = \lambda/(\pi w_0)\).

gouy_phase(z)[source]#

Gouy phase at axial position z, \(\zeta(z) = \arctan[(z-z_0)/z_R]\).

Parameters:: z (float or array_like) – Axial position(s).
Returns:: float or ndarray

q_parameter(z)[source]#

Complex beam parameter at axial position z, \(q(z) = (z-z_0) + i z_R\).

Parameters:: z (float) – Axial position.
Returns:: complex

radius_of_curvature(z)[source]#

Wavefront radius of curvature at axial position z.

\[R(z) = (z - z_0)\left[1 + \left(\frac{z_R}{z - z_0}\right)^2\right]\]

with \(R(z_0) = \infty\) (the wavefront is flat at the waist), handled explicitly rather than relying on the indeterminate-form arithmetic of the formula above.

Parameters:: z (float or array_like) – Axial position(s).
Returns:: float or ndarray

property rayleigh_range#: Rayleigh range, \(z_R = \pi w_0^2/\lambda\).

waist(z)[source]#

Beam radius at axial position z.

\[w(z) = w_0 \sqrt{1 + \left(\frac{z - z_0}{z_R}\right)^2}\]

Parameters:: z (float or array_like) – Axial position(s).
Returns:: float or ndarray

physicskit.optics.gaussian.hermite_gaussian_mode(x, y, z, beam, m, n)[source]#

Hermite-Gaussian \(\mathrm{TEM}_{mn}\) mode amplitude.

\[u_{mn}(x, y, z) = \frac{w_0}{w(z)}\, H_m\!\left(\frac{\sqrt2\,x}{w(z)}\right) H_n\!\left(\frac{\sqrt2\,y}{w(z)}\right) \exp\!\left[-\frac{x^2+y^2}{w(z)^2}\right] \exp\!\left[-\frac{ik(x^2+y^2)}{2R(z)}\right] \exp\!\left[i(m+n+1)\zeta(z)\right] \exp(-ikz)\]

with \(k = 2\pi/\lambda\). \(m=n=0\) reduces to the fundamental Gaussian mode carried by beam.

Parameters:

x (array_like) – Transverse coordinates (broadcastable).
y (array_like) – Transverse coordinates (broadcastable).
z (float) – Axial position.
beam (GaussianBeam) – The underlying fundamental-mode beam (sets \(w_0\), \(z_0\), \(\lambda\)).
m (int) – Transverse mode indices along \(x\) and \(y\).
n (int) – Transverse mode indices along \(x\) and \(y\).

Returns:

ndarray – Complex field amplitude, broadcast shape of x and y.

Examples

>>> import numpy as np
>>> beam = GaussianBeam(wavelength=1.0, w0=1.0, z0=0.0)
>>> u00 = hermite_gaussian_mode(0.0, 0.0, 0.0, beam, 0, 0)
>>> round(float(abs(u00)), 6)
1.0

physicskit.optics.gaussian.laguerre_gaussian_mode(r, phi, z, beam, l, p)[source]#

Laguerre-Gaussian \(\mathrm{LG}_p^l\) mode amplitude.

\[u_{lp}(r, \phi, z) = \frac{w_0}{w(z)} \left(\frac{r\sqrt2}{w(z)}\right)^{|l|} L_p^{|l|}\!\left(\frac{2r^2}{w(z)^2}\right) \exp\!\left[-\frac{r^2}{w(z)^2}\right] \exp\!\left[-\frac{ikr^2}{2R(z)}\right] \exp(il\phi) \exp\!\left[i(|l|+2p+1)\zeta(z)\right] \exp(-ikz)\]

with \(k = 2\pi/\lambda\) and \(L_p^{|l|}\) the associated Laguerre polynomial.

Parameters:

r (array_like) – Polar transverse coordinates (broadcastable); r is the radial distance from the axis, phi the azimuthal angle.
phi (array_like) – Polar transverse coordinates (broadcastable); r is the radial distance from the axis, phi the azimuthal angle.
z (float) – Axial position.
beam (GaussianBeam) – The underlying fundamental-mode beam.
l (int) – Azimuthal (orbital) mode index.
p (int) – Radial mode index.

Returns:

ndarray – Complex field amplitude, broadcast shape of r and phi.

Examples

>>> import numpy as np
>>> beam = GaussianBeam(wavelength=1.0, w0=1.0, z0=0.0)
>>> u00 = laguerre_gaussian_mode(0.0, 0.0, 0.0, beam, 0, 0)
>>> round(float(abs(u00)), 6)
1.0

physicskit.optics.gaussian.m2_beam_waist(z, wavelength, w0, M2=1.0)[source]#

Beam radius of a non-ideal (\(M^2 > 1\)) real laser beam, waist fixed at \(z=0\).

\[w(z) = w_0 \sqrt{1 + \left(\frac{z}{z_{R,\text{eff}}}\right)^2}, \qquad z_{R,\text{eff}} = \frac{\pi w_0^2}{M^2 \lambda}\]

A beam-quality factor \(M^2 \ge 1\) (equal to 1 for an ideal diffraction-limited Gaussian beam) reduces the effective Rayleigh range, so a real beam of the same waist diverges faster than the ideal case.

Parameters:

z (float or array_like) – Axial position(s), measured from the waist at \(z=0\).
wavelength (float) – Wavelength.
w0 (float) – Waist radius.
M2 (float, default=1.0) – Beam-quality factor, \(M^2 \ge 1\).

Returns:

float or ndarray – Beam radius w(z).

Examples

>>> round(float(m2_beam_waist(0.0, wavelength=0.5e-3, w0=0.1, M2=2.0)), 6)
0.1

physicskit.optics.gaussian.propagate_q(q, M)[source]#

Propagate a complex beam parameter through a paraxial optical system.

Applies the same bilinear transformation used for ray-transfer (ABCD) matrices,

\[q_{\text{out}} = \frac{A q + B}{C q + D},\]

with \(A, B, C, D\) the entries of M in the convention M = [[A, B], [C, D]].

Parameters:

q (complex) – Beam parameter before the system.
M (ndarray of shape (2, 2)) – Ray-transfer (ABCD) matrix of the optical system.

Returns:

complex – Beam parameter q_out after the system.

Examples

>>> import numpy as np
>>> free_space = np.array([[1.0, 2.0], [0.0, 1.0]])
>>> propagate_q(1j, free_space)
(2+1j)

physicskit.optics.gaussian.q_to_beam_params(q, wavelength)[source]#

Extract the beam radius and radius of curvature from a complex beam parameter.

Inverts

\[\frac{1}{q} = \frac{1}{R} - i\,\frac{\lambda}{\pi w^2}.\]

At the beam waist \(\mathrm{Re}(1/q) = 0\), so R would formally diverge; this is handled explicitly and np.inf is returned there instead of raising or silently producing nan.

Parameters:

q (complex) – Complex beam parameter.
wavelength (float) – Wavelength (in the same length units as q).

Returns:

w (float) – Beam radius (\(1/e^2\) intensity radius).
R (float) – Radius of curvature of the wavefront; np.inf at the waist.

Examples

>>> import numpy as np
>>> zR = 3.0
>>> w, R = q_to_beam_params(1j * zR, wavelength=0.5e-3)
>>> round(float(w), 6), R
(0.021851, inf)

Quantum optics in the truncated Fock basis: states, Wigner functions, and Jaynes-Cummings dynamics.

Light is quantized in terms of the harmonic-oscillator (photon-number, or Fock) states \(\lvert n\rangle\) of a single cavity mode. This module works throughout in a Fock basis truncated to a finite dimension cutoff – large enough that the state’s amplitude on the highest retained level is negligible – and provides:

constructors for the number, coherent, and squeezed states (built from the ladder operators of physicskit.quantum.core.operators);
the Wigner quasi-probability distribution, computed from its closed-form matrix elements in the Fock basis (Cahill & Glauber, Phys. Rev. 177, 1857, 1969), using the dimensionless quadratures

\[\hat x = \frac{\hat a + \hat a^\dagger}{\sqrt2}, \qquad \hat p = \frac{\hat a - \hat a^\dagger}{i\sqrt2};\]
the Jaynes-Cummings model (Jaynes & Cummings, Proc. IEEE 51, 89, 1963), the fully quantum treatment of a two-level atom coupled to a single cavity mode in the rotating-wave approximation, famous for predicting vacuum Rabi oscillations and collapse-and-revival – effects with no semiclassical analogue.

class physicskit.optics.quantum_optics.JaynesCummingsModel(omega_c, omega_a, g, cutoff=10)[source]#

Bases: object

The Jaynes-Cummings model: a two-level atom coupled to a single quantized cavity mode.

In the rotating-wave approximation, the Hamiltonian is

\[\hat H = \omega_c\, \hat a^\dagger \hat a \otimes \hat I + \frac{\omega_a}{2}\, \hat I \otimes \hat\sigma_z + g\left(\hat a \otimes \hat\sigma_+ + \hat a^\dagger \otimes \hat\sigma_-\right)\]

with the cavity (dimension cutoff) as the left tensor factor and the atom (dimension 2, basis \(\{\lvert e\rangle, \lvert g\rangle\}\), so \(\hat\sigma_z = \mathrm{diag}(+1,-1)\), \(\hat\sigma_+ = \lvert e\rangle\langle g\rvert\), \(\hat\sigma_- = \lvert g\rangle\langle e\rvert\)) as the right factor, i.e. H = kron(H_cav, H_atom) terms throughout. A combined basis state \(\lvert n, e/g\rangle\) therefore sits at flat index 2*n (excited) or 2*n+1 (ground) of the 2*cutoff-dimensional Hilbert space.

Parameters:

omega_c (float) – Cavity mode frequency.
omega_a (float) – Atomic transition frequency.
g (float) – Atom-cavity coupling strength.
cutoff (int, default=10) – Truncation of the cavity Fock space.

evolve(psi0, t_array)[source]#

Time-evolve a state under the Jaynes-Cummings Hamiltonian.

Diagonalizes \(\hat H\) once and reconstructs \(\lvert\psi(t)\rangle = e^{-i\hat H t}\lvert\psi_0\rangle\) at every requested time from the eigendecomposition, which is much faster than exponentiating \(\hat H\) separately for each t.

Parameters:

psi0 (ndarray of shape (2*cutoff,)) – Initial state.
t_array (ndarray of shape (T,)) – Times at which to evaluate the evolved state.

Returns:

ndarray of shape (T, 2*cutoff) – Complex states \(\lvert\psi(t)\rangle\) at each requested time.

excited_state_population(t_array, n_photons=0)[source]#

Atomic excited-state population \(P_e(t)\) starting from \(\lvert e, n\rangle\).

Starts from the initial state \(\lvert e, n_\text{photons}\rangle\) (atom excited, n_photons photons in the cavity) and returns

\[P_e(t) = \sum_n \left\lvert \langle e, n \rvert \psi(t)\rangle \right\rvert^2.\]

On resonance (\(\omega_c = \omega_a\)) with n_photons=0, this reduces to the textbook vacuum Rabi formula \(P_e(t) = \cos^2(gt)\).

Parameters:

t_array (ndarray of shape (T,)) – Times at which to evaluate the population.
n_photons (int, default=0) – Initial photon number (with the atom excited).

Returns:

ndarray of shape (T,) – Real-valued excited-state population at each time.

Examples

>>> import numpy as np
>>> jc = JaynesCummingsModel(omega_c=1.0, omega_a=1.0, g=0.5, cutoff=10)
>>> t = np.array([0.0, np.pi / 2 / 0.5])
>>> Pe = jc.excited_state_population(t, n_photons=0)
>>> np.round(Pe, 6)
array([1., 0.])

hamiltonian()[source]#

The Jaynes-Cummings Hamiltonian as a dense matrix.

Returns:: ndarray of shape (2*cutoff, 2*cutoff) – The Hamiltonian \(\hat H\), in the cavity (x) atom basis ordering documented on the class.

Examples

>>> import numpy as np
>>> jc = JaynesCummingsModel(omega_c=1.0, omega_a=1.0, g=0.1, cutoff=3)
>>> H = jc.hamiltonian()
>>> H.shape
(6, 6)
>>> bool(np.allclose(H, H.conj().T))
True

physicskit.optics.quantum_optics.coherent_state(alpha, cutoff)[source]#

The coherent state \(\lvert\alpha\rangle\), truncated to a finite Fock basis.

\[\lvert\alpha\rangle = e^{-\lvert\alpha\rvert^2/2} \sum_{n=0}^{\infty} \frac{\alpha^n}{\sqrt{n!}}\,\lvert n\rangle\]

The sum is truncated to cutoff terms and the result renormalized, so it is a good approximation only while \(\lvert\alpha\rvert^2 \ll\) cutoff (i.e. the mean photon number is well within the truncated space). Amplitudes are built iteratively in log-space via scipy.special.gammaln() to avoid overflow for large cutoff.

Parameters:

alpha (complex) – Coherent-state amplitude.
cutoff (int) – Dimension of the truncated Fock space.

Returns:

ndarray of shape (cutoff,) – Complex, L2-normalized state vector.

Examples

>>> psi = coherent_state(0.0, 8)
>>> psi
array([1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])

physicskit.optics.quantum_optics.compute_wigner_function(state_vector, x_grid, p_grid)[source]#

Wigner quasi-probability distribution of a Fock-basis state.

Computed from the density matrix \(\rho = \lvert\psi\rangle\langle\psi\rvert\) and the closed-form Fock-basis Wigner matrix elements,

\[W(x,p) = \sum_{n,m} \rho_{nm}\, K_{nm}(x,p),\]

with \(K_{nm}\) built from associated Laguerre polynomials of \(2(x^2+p^2)\) (see scipy.special.eval_genlaguerre()). Uses dimensionless quadratures \(x = (a+a^\dagger)/\sqrt2\), \(p = (a-a^\dagger)/(i\sqrt2)\), for which \(\iint W(x,p)\,dx\,dp = 1\) for any normalized state.

Parameters:

state_vector (ndarray of shape (cutoff,)) – Fock-basis ket.
x_grid (ndarray of shape (Nx,)) – Grid of \(x\) quadrature values.
p_grid (ndarray of shape (Np,)) – Grid of \(p\) quadrature values.

Returns:

ndarray of shape (Nx, Np) – Real-valued Wigner function \(W(x,p)\).

Examples

>>> import numpy as np
>>> vac = fock_state(0, 8)
>>> xg = np.linspace(-4, 4, 41)
>>> pg = np.linspace(-4, 4, 41)
>>> W = compute_wigner_function(vac, xg, pg)
>>> W.shape
(41, 41)
>>> round(float(W[20, 20] * np.pi), 6)
1.0

physicskit.optics.quantum_optics.fock_state(n, cutoff)[source]#

The number (Fock) state \(\lvert n\rangle\).

Parameters:

n (int) – Photon number.
cutoff (int) – Dimension of the truncated Fock space.

Returns:

ndarray of shape (cutoff,) – Complex state vector, all zeros except a 1 at index n.

Raises:

ValueError – If n is negative or n >= cutoff.

Examples

>>> fock_state(1, 4)
array([0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j])

physicskit.optics.quantum_optics.squeezed_state(xi, alpha=0.0, cutoff=30)[source]#

The displaced squeezed vacuum state, using the squeeze-then-displace convention.

Constructs

\[\lvert\xi, \alpha\rangle = \hat D(\alpha)\, \hat S(\xi)\, \lvert 0\rangle, \qquad \hat S(\xi) = \exp\!\left[\frac{\xi^* \hat a^2 - \xi \hat a^{\dagger 2}}{2}\right], \qquad \hat D(\alpha) = \exp\!\left(\alpha \hat a^\dagger - \alpha^* \hat a\right)\]

i.e. the vacuum is squeezed first and then displaced. (The opposite order, displace-then-squeeze, gives a different, also physically valid, state – squeeze-then-displace is the convention adopted here and is the more common one in the literature.) The ladder operators are built on the same truncated Fock basis, so this is only accurate while the squeezed/displaced amplitude is small relative to cutoff.

Parameters:

xi (complex) – Squeezing parameter \(\xi = r e^{i\theta}\).
alpha (complex, default=0.0) – Displacement amplitude.
cutoff (int, default=30) – Dimension of the truncated Fock space.

Returns:

ndarray of shape (cutoff,) – Complex, L2-normalized state vector.

Examples

>>> psi = squeezed_state(0.0, alpha=0.0, cutoff=6)
>>> psi
array([1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])

physicskit.optics.quantum_optics.wigner_negativity(W, x_grid, p_grid)[source]#

Wigner negativity: the integrated negative volume of a Wigner function.

\[N_W = \iint \max(0, -W(x,p))\; dx\, dp\]

A standard non-classicality diagnostic: zero for classical/Gaussian states (coherent, thermal, squeezed vacuum), strictly positive for genuinely non-classical states such as Fock states with \(n \ge 1\).

Parameters:

W (ndarray of shape (Nx, Np)) – Wigner function, e.g. from compute_wigner_function().
x_grid (ndarray of shape (Nx,)) – Grid of \(x\) quadrature values.
p_grid (ndarray of shape (Np,)) – Grid of \(p\) quadrature values.

Returns:

float – The (non-negative) integrated negative volume.

Examples

>>> import numpy as np
>>> vac = fock_state(0, 8)
>>> xg = np.linspace(-5, 5, 61)
>>> pg = np.linspace(-5, 5, 61)
>>> W = compute_wigner_function(vac, xg, pg)
>>> wigner_negativity(W, xg, pg)
0.0

Interference in layered media: thin films, Fabry-Pérot cavities and 1D photonic crystals.

A plane wave crossing a stack of homogeneous layers is described exactly by the characteristic matrix of each layer (Abelès, Ann. Phys. (Paris) 5, 596 (1950); Born and Wolf, Principles of Optics, sec. 1.6),

\[\begin{split}M_j = \begin{pmatrix} \cos\delta_j & i\sin\delta_j/\eta_j \\ i\eta_j\sin\delta_j & \cos\delta_j \end{pmatrix}, \qquad \delta_j = \frac{2\pi}{\lambda} n_j d_j \cos\theta_j,\end{split}\]

with tilted admittance \(\eta_j = n_j\cos\theta_j\) (s polarization) or \(n_j/\cos\theta_j\) (p). The product of the layer matrices gives the reflection and transmission of the whole stack. Three classic results follow:

a single quarter-wave film of index \(\sqrt{n_0 n_s}\) cancels reflection at one wavelength (antireflection coating);
two mirrors of reflectance \(R\) facing each other transmit the Airy function \(T = 1/(1 + F\sin^2(\delta/2))\), \(F = 4R/(1-R)^2\) (Fabry and Pérot, 1899);
a periodic stack has Bloch modes with \(\cos(K\Lambda) = \tfrac12\operatorname{tr} M_{\text{period}}\), and where \(|\tfrac12\operatorname{tr} M| > 1\) light cannot propagate: a photonic band gap (Rayleigh 1887; Yablonovitch 1987).

physicskit.optics.thin_films.airy_transmission(delta, R)[source]#

Fabry-Pérot (Airy) transmission of two lossless mirrors of reflectance \(R\).

\[T = \frac{1}{1 + F\sin^2(\delta/2)}, \qquad \delta = \frac{4\pi n L\cos\theta}{\lambda},\]

where \(\delta\) is the round-trip phase.

Parameters:

delta (ArrayLike) – Round-trip phase.
R (float) – Mirror reflectance.

Return type:

NDArray

Returns:

ndarray

Examples

>>> float(airy_transmission(0.0, 0.99))  # every resonance transmits fully
1.0

physicskit.optics.thin_films.bloch_wavenumber(n_layers, d_layers, wavelength)[source]#

Bloch wavenumber \(K\) of an infinite periodic stack at normal incidence.

For one period of layers, \(\cos(K\Lambda) = \tfrac12 \operatorname{tr} M_{\text{period}}\); for two layers this is

\[\cos(K\Lambda) = \cos\delta_1\cos\delta_2 - \tfrac12\left(\frac{n_1}{n_2} + \frac{n_2}{n_1}\right)\sin\delta_1\sin\delta_2.\]

Parameters:

n_layers (ArrayLike) – Indices and thicknesses of the layers of one period.
d_layers (ArrayLike) – Indices and thicknesses of the layers of one period.
wavelength (ArrayLike) – Vacuum wavelength(s).

Return type:

NDArray[cdouble]

Returns:

ndarray of complex – \(K\Lambda\), with real part in \([0, \pi]\) and a nonzero imaginary part (the decay per period) inside band gaps.

Examples

>>> KL = bloch_wavenumber([1.0, 1.0], [0.25, 0.25], 1.0)  # uniform medium: K = 2 pi n / lambda
>>> round(float(KL.real), 6) == round(np.pi, 6)
True

physicskit.optics.thin_films.coefficient_of_finesse(R)[source]#

\(F = 4R/(1 - R)^2\).

Parameters:: R (ArrayLike) – Mirror power reflectance, \(0 \le R < 1\).
Return type:: NDArray
Returns:: ndarray

physicskit.optics.thin_films.fabry_perot_free_spectral_range(length, n=1.0, wavelength=None)[source]#

Free spectral range \(c/2nL\) in frequency, or \(\lambda^2/2nL\) in wavelength.

Parameters:

length (float) – Mirror spacing \(L\).
n (float) – Index between the mirrors.
wavelength (float | None) – If given, return the range in wavelength units near this wavelength; otherwise in units of \(c\) per length unit.

Return type:

Returns:

float

physicskit.optics.thin_films.finesse(R)[source]#

Reflectance finesse \(\mathcal F = \pi\sqrt R/(1 - R)\): free spectral range over linewidth.

Parameters:: R (ArrayLike) – Mirror power reflectance.
Return type:: NDArray
Returns:: ndarray

Examples

>>> round(float(finesse(0.9)), 2)
29.8

physicskit.optics.thin_films.multilayer_response(n_layers, d_layers, wavelength, n_incident=1.0, n_substrate=1.5, angle=0.0, polarization='s')[source]#

Reflection and transmission of a stack of thin layers by the characteristic-matrix method.

Parameters:

n_layers (ArrayLike) – Refractive index of each layer, from the incident side; complex values \(n + i\kappa\) with \(\kappa > 0\) describe absorption (time dependence \(e^{-i\omega t}\)).
d_layers (ArrayLike) – Physical thickness of each layer, in the units of wavelength.
wavelength (ArrayLike) – Vacuum wavelength(s).
n_incident (complex) – Index of the incident medium.
n_substrate (complex) – Index of the exit medium.
angle (float) – Angle of incidence in radians.
polarization (str)

Return type:

Returns:

dict – Complex amplitude coefficients "r", "t" and the power reflectance "R" and transmittance "T", each shaped like wavelength.

Examples

Bare glass reflects \(((1 - 1.5)/(1 + 1.5))^2 = 4\%\); a quarter-wave film of index \(\sqrt{1.5}\) removes it at the design wavelength:

>>> round(float(multilayer_response([], [], 550.0)["R"]), 4)
0.04
>>> n_ar = np.sqrt(1.5)
>>> float(multilayer_response([n_ar], [550.0 / (4 * n_ar)], 550.0)["R"]) < 1e-12
True

physicskit.optics.thin_films.quarter_wave_band_gap(n_high, n_low)[source]#

Relative width \(\Delta\omega/\omega_0\) of the first gap of a quarter-wave stack.

\[\frac{\Delta\omega}{\omega_0} = \frac4\pi\arcsin\frac{n_H - n_L}{n_H + n_L}\]

(Yariv and Yeh, Optical Waves in Crystals, 1984, sec. 6.3).

Parameters:

n_high (float)
n_low (float)

Return type:

Returns:

float

Examples

>>> round(quarter_wave_band_gap(2.3, 1.38), 4)
0.3217

physicskit.optics.thin_films.quarter_wave_stack(n_high, n_low, n_pairs, design_wavelength, cavity=False)[source]#

Layer indices and thicknesses of a quarter-wave (Bragg) mirror, or of a Fabry-Pérot cavity.

Parameters:

n_high (float) – Indices of the two materials; the stack starts with n_high.
n_low (float) – Indices of the two materials; the stack starts with n_high.
n_pairs (int) – Number of high/low pairs.
design_wavelength (float) – Wavelength at which every layer is a quarter wave thick.
cavity (bool) – If True, return two such mirrors separated by a half-wave (\(\lambda/2n_{\text{low}}\)) spacer of the low index: a Fabry-Pérot cavity with a resonance at design_wavelength.

Return type:

tuple[NDArray, NDArray]

Returns:

n_layers, d_layers (ndarray)

Examples

>>> n, d = quarter_wave_stack(2.3, 1.38, 2, 600.0)
>>> n.real.tolist(), (n.real * d).tolist()
([2.3, 1.38, 2.3, 1.38], [150.0, 150.0, 150.0, 150.0])

Second-harmonic generation and phase matching.

In a medium with a second-order susceptibility \(\chi^{(2)}\), a field at frequency \(\omega\) drives a polarization at \(2\omega\) (Franken et al., Phys. Rev. Lett. 7, 118 (1961)). For collinear plane waves with slowly varying amplitudes \(A_1\) (fundamental) and \(A_2\) (second harmonic), normalized so that \(|A_1|^2 + |A_2|^2\) is the conserved power, the coupled-amplitude equations are (Armstrong, Bloembergen, Ducuing and Pershan, Phys. Rev. 127, 1918 (1962); Boyd, Nonlinear Optics, ch. 2)

\[\frac{dA_1}{dz} = -i\kappa A_2 A_1^* e^{-i\Delta k z}, \qquad \frac{dA_2}{dz} = -i\kappa A_1^2 e^{i\Delta k z}, \qquad \Delta k = k_{2\omega} - 2k_\omega.\]

Three regimes are solved in closed form:

undepleted pump, \(|A_2(L)|^2 = \kappa^2|A_1|^4 L^2 \operatorname{sinc}^2(\Delta k L/2)\), which oscillates with the coherence length \(L_c = \pi/\Delta k\);
perfect phase matching with pump depletion, \(|A_2|^2 = |A_0|^2\tanh^2(\kappa A_0 z)\);
quasi-phase matching, where the sign of \(\kappa\) is flipped every \(L_c\) and the harmonic grows linearly at \(2/\pi\) of the phase-matched rate.

Birefringent phase matching uses the extraordinary index \(n_e(\theta)\) to make \(n_e^{2\omega}(\theta) = n_o^{\omega}\) (Giordmaine 1962; Maker et al. 1962).

physicskit.optics.nonlinear.extraordinary_index(theta, n_o, n_e)[source]#

Index of the extraordinary wave at angle \(\theta\) to the optic axis of a uniaxial crystal.

\[\frac{1}{n_e(\theta)^2} = \frac{\cos^2\theta}{n_o^2} + \frac{\sin^2\theta}{n_e^2}\]

Parameters:

theta (ArrayLike) – Angle in radians.
n_o (float) – Ordinary and principal extraordinary indices.
n_e (float) – Ordinary and principal extraordinary indices.

Return type:

Returns:

ndarray

physicskit.optics.nonlinear.shg_coupled_amplitudes(z, delta_k=0.0, kappa=1.0, A0=1.0, qpm_period=None)[source]#

Integrate the SHG coupled-amplitude equations, with optional quasi-phase matching.

Parameters:

z (ArrayLike) – Increasing output positions, starting at 0.
delta_k (float) – Phase mismatch \(\Delta k\).
kappa (float) – Nonlinear coupling.
A0 (complex) – Fundamental amplitude at \(z = 0\); the harmonic starts at 0.
qpm_period (float | None) – If given, \(\kappa\) changes sign every half period (periodically poled crystal); first-order QPM uses \(2\pi/\Delta k\).

Return type:

Returns:

dict – Complex "A1" and "A2" at z.

Examples

>>> out = shg_coupled_amplitudes(np.linspace(0, 2, 5))
>>> bool(np.allclose(np.abs(out["A2"]) ** 2, np.tanh(np.linspace(0, 2, 5)) ** 2, atol=1e-6))
True

physicskit.optics.nonlinear.shg_phase_matched_efficiency(z, kappa=1.0, A0=1.0)[source]#

Conversion efficiency \(|A_2|^2/|A_0|^2 = \tanh^2(\kappa A_0 z)\) with pump depletion.

Parameters:

z (ArrayLike) – Propagation distance(s).
kappa (float)
A0 (float) – Input fundamental amplitude.

Return type:

Returns:

ndarray

Examples

>>> round(float(shg_phase_matched_efficiency(10.0)), 6)
1.0

physicskit.optics.nonlinear.shg_undepleted_power(L, delta_k, kappa=1.0, A1=1.0)[source]#

Second-harmonic power for an undepleted pump.

\[|A_2(L)|^2 = \kappa^2 |A_1|^4 L^2 \operatorname{sinc}^2(\Delta k L / 2) = \frac{4\kappa^2|A_1|^4}{\Delta k^2}\sin^2\frac{\Delta k L}{2}.\]

Parameters:

L (ArrayLike) – Crystal length(s).
delta_k (float) – Phase mismatch.
kappa (float) – Coupling constant.
A1 (float) – Pump amplitude.

Return type:

Returns:

ndarray

Examples

>>> round(float(shg_undepleted_power(np.pi, delta_k=1.0)), 12)  # one coherence length: 4/dk^2
4.0

physicskit.optics.nonlinear.type_i_phase_matching_angle(n_o_w, n_o_2w, n_e_2w)[source]#

Type-I (ooe) phase-matching angle of a negative uniaxial crystal.

Solves \(n_e^{2\omega}(\theta) = n_o^{\omega}\):

\[\sin^2\theta_m = \frac{(n_o^{\omega})^{-2} - (n_o^{2\omega})^{-2}} {(n_e^{2\omega})^{-2} - (n_o^{2\omega})^{-2}}.\]

Parameters:

n_o_w (float) – Ordinary index at the fundamental.
n_o_2w (float) – Ordinary and extraordinary indices at the harmonic.
n_e_2w (float) – Ordinary and extraordinary indices at the harmonic.

Return type: