physicskit.optics

physicskit.optics hub · History · Examples · API reference

physicskit.optics: classical and quantum optics – rays, waves, Gaussian beams, and quantized light.

Typical usage:

import physicskit as pk
import numpy as np

beam = pk.optics.GaussianBeam(wavelength=1.064e-3, w0=0.05)
system = pk.optics.OpticalSystem([
    pk.optics.OpticalElement(pk.optics.free_space(1.0), length=1.0),
    pk.optics.OpticalElement(pk.optics.thin_lens(0.5)),
])
q_out = pk.optics.propagate_q(beam.q_parameter(0.0), system.system_matrix())

• physicskit.optics.ray – geometric ray optics and ABCD matrices: free space, thin/thick lenses, curved and flat interfaces, spherical mirrors, GRIN media, cascaded optical systems, and laser-cavity stability.

• physicskit.optics.wave – scalar wave optics: circular/slit apertures and Fraunhofer, Fresnel, and angular-spectrum diffraction.

• physicskit.optics.gaussian – Gaussian beam propagation via the complex beam parameter \(q\), Hermite-Gaussian and Laguerre-Gaussian higher-order modes, and \(M^2\) beam-quality propagation.

• physicskit.optics.quantum_optics – quantum optical states (Fock, coherent, squeezed) in a truncated photon-number basis, Wigner quasi-probability distributions, and the Jaynes-Cummings model of a two-level atom coupled to a quantized cavity mode.

• physicskit.optics.visualizers – ray-trace, beam-envelope, diffraction-pattern, and interactive Wigner-surface plots.

class physicskit.optics.GaussianBeam(wavelength, w0, z0=0.0)

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)

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)

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)

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)

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

class physicskit.optics.JaynesCummingsModel(omega_c, omega_a, g, cutoff=10)

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)

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)

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()

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

class physicskit.optics.OpticalElement(matrix, name='', length=0.0)

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.OpticalSystem(elements)

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()

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()

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)

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.angular_spectrum_propagate(U0, wavelength, z, dx)

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.animate_diffraction_propagation(aperture, wavelength, z_values, dx, log_scale=False, interval=100, ax=None)

Animate the diffraction pattern developing as propagation distance z increases.

Steps physicskit.optics.wave.angular_spectrum_propagate() over each distance in z_values (treating z as the animation’s “time” axis – the standard way to visualize the Fresnel-to-Fraunhofer development of a diffraction pattern), and animates the resulting intensity as an imshow heatmap: a double-slit aperture’s near-field wavefronts visibly evolve into the far-field interference fringes as z grows.

Parameters:

• aperture (ndarray of shape (Ny, Nx)) – Input field immediately after the aperture (e.g. from physicskit.optics.wave.double_slit_aperture()).

• wavelength (float) – Wavelength.

• z_values (ndarray) – Sequence of propagation distances to sweep over.

• dx (float) – Grid spacing (same for input and output, since angular_spectrum_propagate() stays on the same grid).

• log_scale (bool, default=False) – If True, plot log10(intensity + eps) to reveal faint fringes.

• interval (int, default=100) – Delay between frames, in milliseconds.

• ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.optics.wave import double_slit_aperture
>>> ap = double_slit_aperture((64, 64), dx=1e-3, width=2e-3, separation=1e-2).astype(complex)
>>> z_values = np.linspace(0.01, 2.0, 5)
>>> anim = animate_diffraction_propagation(ap, wavelength=0.5e-3, z_values=z_values, dx=1e-3)
>>> isinstance(anim, FuncAnimation)
True

physicskit.optics.cavity_round_trip_matrix(elements)

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.cavity_stability(M)

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.circular_aperture(shape, dx, radius)

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.coherent_state(alpha, cutoff)

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.compute_wigner_function(state_vector, x_grid, p_grid)

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.curved_interface(R, n1, n2)

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.double_slit_aperture(shape, dx, width, separation)

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.field_grid(shape, dx)

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.flat_interface(n1, n2)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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.plot_beam_envelope(beam, z_range, n_points=200, ax=None)

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)

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)

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)

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)

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.single_slit_aperture(shape, dx, width)

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)

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)

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)

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)

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.wigner_negativity(W, x_grid, p_grid)

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)

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)

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()

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()

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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()

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)

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)

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)

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)

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)

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

Plotting helpers: ray traces, Gaussian beam envelopes, diffraction patterns, Wigner surfaces.

Matplotlib is used for 2D plots (ray traces through an optical system, beam waist envelopes, diffraction-pattern heatmaps), and Plotly for the interactive 3D Wigner phase-space surface. Every function returns its figure object rather than calling show(), so it composes cleanly into larger figures or headless pipelines.

physicskit.optics.visualizers.animate_diffraction_propagation(aperture, wavelength, z_values, dx, log_scale=False, interval=100, ax=None)

Animate the diffraction pattern developing as propagation distance z increases.

Steps physicskit.optics.wave.angular_spectrum_propagate() over each distance in z_values (treating z as the animation’s “time” axis – the standard way to visualize the Fresnel-to-Fraunhofer development of a diffraction pattern), and animates the resulting intensity as an imshow heatmap: a double-slit aperture’s near-field wavefronts visibly evolve into the far-field interference fringes as z grows.

Parameters:

• aperture (ndarray of shape (Ny, Nx)) – Input field immediately after the aperture (e.g. from physicskit.optics.wave.double_slit_aperture()).

• wavelength (float) – Wavelength.

• z_values (ndarray) – Sequence of propagation distances to sweep over.

• dx (float) – Grid spacing (same for input and output, since angular_spectrum_propagate() stays on the same grid).

• log_scale (bool, default=False) – If True, plot log10(intensity + eps) to reveal faint fringes.

• interval (int, default=100) – Delay between frames, in milliseconds.

• ax (matplotlib.axes.Axes, optional) – Axes to draw into; a new figure is created if omitted.

Returns:

matplotlib.animation.FuncAnimation

Examples

>>> import numpy as np
>>> from physicskit.optics.wave import double_slit_aperture
>>> ap = double_slit_aperture((64, 64), dx=1e-3, width=2e-3, separation=1e-2).astype(complex)
>>> z_values = np.linspace(0.01, 2.0, 5)
>>> anim = animate_diffraction_propagation(ap, wavelength=0.5e-3, z_values=z_values, dx=1e-3)
>>> isinstance(anim, FuncAnimation)
True

physicskit.optics.visualizers.interactive_wigner_surface(W, x_grid, p_grid, title=None)

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.visualizers.plot_beam_envelope(beam, z_range, n_points=200, ax=None)

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.visualizers.plot_diffraction_pattern(U, dx, ax=None, log_scale=False)

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.visualizers.plot_ray_trace(system, y0, theta0, ax=None)

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