Matrix Quantum Optics

• Matrix quantum optics is a framework that expresses optical phenomena—such as propagation and mode transformation—as structured matrices including transfer, unitary, and symplectic forms.
• It enables practical applications like quantum state tomography and emulation of quantum operations through physically realizable linear-optical devices and tensor-network simulations.
• The approach spans various formulations, from Gaussian and density matrix representations to tensor-network methods, offering insights into both classical and nonlinear photonic systems.

Searching arXiv for recent and foundational uses of “matrix quantum optics” and closely related matrix-based optical formalisms. Using the available arXiv corpus in the provided materials, the article synthesizes foundational and recent usages spanning matrix optics, Gaussian/symplectic quantum optics, optical matrix multiplication, tomography of linear optical devices, and tensor-network formulations. Matrix quantum optics denotes a family of matrix-centered formulations in which optical propagation, photonic mode transformations, Gaussian-state dynamics, and quantum-state evolution are written directly in terms of transfer matrices, unitary mode matrices, symplectic matrices, density matrices, or tensor-network objects. In the cited literature, passive linear optical circuits are represented by unitary mode transformations, paraxial and Maxwell propagation are recast as operator-valued matrix evolutions, continuous-variable Gaussian operations are organized by structured matrix decompositions, structured light realizes matrix–vector implementations of quantum gates, and broadband or multimode nonlinear optics is simulated with matrix product states (Laing et al., 2012, Houde et al., 2024, Koni et al., 2024, Yanagimoto et al., 2021, Cilluffo et al., 3 Feb 2025).

1. Linear-optical networks as matrix transformations

At the most standard level, matrix quantum optics is the description of a passive linear optical device by a matrix acting on optical modes. An ideal, lossless, passive linear optical circuit with $m$ modes is described by an $m\times m$ unitary matrix $U$, and the Heisenberg-picture mode transformation is $$\hat a_j^{(\mathrm{out})}=\sum_{k=1}^m U_{jk}\hat a_k^{(\mathrm{in})}$$. For $n$ indistinguishable photons in $m$ modes, the relevant Hilbert space dimension is $M=\binom{m+n-1}{n}$, and the same physical interferometer induces a unitary $U\in U(M)$ acting on the $n$-photon sector (Laing et al., 2012, Garcia-Escartin et al., 2019).

This matrix viewpoint is not merely descriptive. The tomography scheme of Laing and O’Brien reconstructs the unitary matrix of an arbitrary linear optical device from one-photon count rates and two-photon visibilities, while remaining independent of input and output port losses and stable on the arbitrarily increasable length scale of the photon packet. In the same matrix language, real devices are modeled by an effective non-unitary matrix $E=L^{(\mathrm{out})} U L^{(\mathrm{in})}$, with diagonal loss matrices accounting for input and output transmission or detection inefficiency (Laing et al., 2012).

A complementary problem is implementability: given a target unitary on the $m\times m$0-photon Hilbert space, when does it arise from linear optics on $m\times m$1 modes? The construction of García-Escartín and Chamorro-Posada studies the homomorphism from the mode-space scattering matrix $m\times m$2 to the induced unitary $m\times m$3, and gives a criterion based on the adjoint map between the corresponding Lie algebras, together with an inverse construction $m\times m$4 when the target evolution is realizable (Garcia-Escartin et al., 2019). Demirel’s treatment of linear multi-particle, multi-port interferometers gives the same program a concrete lattice form: beam splitters are generated by anti-symmetric generalized Gell-Mann matrices, phase shifters by diagonal generators, and the whole interferometer is written as an ordered product of beam-splitter and phase-shifter matrices for arbitrary port number $m\times m$5 and particle number $m\times m$6 (Demirel, 2019).

2. Operator and transfer-matrix formulations of propagation

A second strand of matrix quantum optics treats propagation itself as a matrix representation of an operator. In paraxial ray optics, the standard $m\times m$7 ABCD matrix for free propagation,

$m\times m$8

is only the first-order truncation of a more general Taylor-expansion formalism. The generalized ABCD matrix acts on the infinite-dimensional vector of Taylor coefficients or derivatives of a ray trajectory,

$m\times m$9

and satisfies $U$0, where $U$1 is the matrix representation of $U$2 in the polynomial basis $U$3. This yields a one-to-one correspondence between generalized ABCD propagation and the quantum mechanical translation operator (Ornigotti et al., 2012).

In this representation, free-space propagation is a translation on the space of Taylor coefficients, and the semigroup property

$U$4

is the direct analogue of composing translation operators. This establishes an operator-theoretic bridge between classical matrix optics and quantum mechanics: the standard ABCD matrix becomes the lowest-order principal submatrix of an exact infinite-dimensional propagator (Ornigotti et al., 2012).

A related, but more field-theoretic, construction is the eight-dimensional matrix representation of Maxwell’s equations based on the Riemann–Silberstein–Weber vectors. For a linear homogeneous medium, the Maxwell operator is written as an $U$5 matrix built from Pauli matrices, is anti-Hermitian, and reduces under unitary transformation to a direct sum of four $U$6 Pauli blocks. In an inhomogeneous medium the operator decomposes into an unperturbed part and a perturbation containing gradients of refractive index and impedance, making perturbative and Foldy–Wouthuysen-like methods available. The same formalism rederives the Mukunda–Simon–Sudarshan substitution rule that lifts scalar paraxial optics to vector Maxwell optics by the replacement $U$7 (Khan et al., 2022).

3. Gaussian, symplectic, and density-matrix formulations

In continuous-variable quantum optics, matrix quantum optics takes the form of symplectic linear algebra. For $U$8 bosonic modes, the quadratures are collected into $U$9, with commutators $$\hat a_j^{(\mathrm{out})}=\sum_{k=1}^m U_{jk}\hat a_k^{(\mathrm{in})}$$0 and symplectic form

$$\hat a_j^{(\mathrm{out})}=\sum_{k=1}^m U_{jk}\hat a_k^{(\mathrm{in})}$$1

Any Gaussian unitary acts as $$\hat a_j^{(\mathrm{out})}=\sum_{k=1}^m U_{jk}\hat a_k^{(\mathrm{in})}$$2, where $$\hat a_j^{(\mathrm{out})}=\sum_{k=1}^m U_{jk}\hat a_k^{(\mathrm{in})}$$3, and covariance matrices transform as $$\hat a_j^{(\mathrm{out})}=\sum_{k=1}^m U_{jk}\hat a_k^{(\mathrm{in})}$$4. In this setting, Gaussian states are covariance matrices and Gaussian unitaries are symplectic matrices (Houde et al., 2024).

The main structured decompositions are then matrix decompositions with direct optical meaning. The Takagi or Autonne decomposition diagonalizes complex symmetric pair-creation matrices, the Bloch–Messiah or Euler decomposition writes any symplectic matrix as $$\hat a_j^{(\mathrm{out})}=\sum_{k=1}^m U_{jk}\hat a_k^{(\mathrm{in})}$$5 with passive interferometers $$\hat a_j^{(\mathrm{out})}=\sum_{k=1}^m U_{jk}\hat a_k^{(\mathrm{in})}$$6 and diagonal squeezing $$\hat a_j^{(\mathrm{out})}=\sum_{k=1}^m U_{jk}\hat a_k^{(\mathrm{in})}$$7, the Iwasawa decomposition gives a unique factorization $$\hat a_j^{(\mathrm{out})}=\sum_{k=1}^m U_{jk}\hat a_k^{(\mathrm{in})}$$8, and Williamson’s theorem diagonalizes any real symmetric positive definite covariance matrix as $$\hat a_j^{(\mathrm{out})}=\sum_{k=1}^m U_{jk}\hat a_k^{(\mathrm{in})}$$9. In optical terms, these are the normal-mode, interferometer–squeezer–interferometer, shear–squeeze–passive, and thermal normal-form decompositions of Gaussian quantum optics (Houde et al., 2024).

A finite-dimensional density-matrix version of matrix quantum optics appears in the unified treatment of coherence and polarization. There the central object is a $n$0 density matrix on $n$1, with basis $n$2. Spatial coherence is extracted from path off-diagonal elements such as $n$3 and $n$4, while polarization at each slit is encoded in $n$5 sub-blocks such as

$n$6

The same matrix supports Stokes parameters, degrees of polarization $n$7, the degree of coherence $n$8, and Kraus-map descriptions of path-only or polarization-dependent dephasing (Bernardo, 2016).

4. Optical realization and emulation of quantum operations

Another usage of matrix quantum optics is operational: classical optical hardware is arranged so that matrix multiplication physically implements a quantum-information transformation. In the structured-light scheme of optical matrix multiplication, the transverse $n$9- and $m$0-coordinates define a discrete Hilbert space $m$1, localized Gaussian spots on an $m$2 lattice provide a basis, one spatial light modulator encodes the input vector $m$3, a second encodes the matrix $m$4, and a cylindrical lens performs the column-wise summation at zero spatial frequency. In this language the optical field plays the role of a state vector, the SLM patterns are operators, and the system realizes $m$5 as a quantum-like map $m$6 (Koni et al., 2024).

Within that framework, Hadamard gates and the Deutsch–Jozsa algorithm were implemented by matrix–vector multiplication on the lattice. The reported fidelities are $m$7 for $m$8, $m$9 for $M=\binom{m+n-1}{n}$0, $M=\binom{m+n-1}{n}$1 for $M=\binom{m+n-1}{n}$2, and $M=\binom{m+n-1}{n}$3 for $M=\binom{m+n-1}{n}$4, with degradation attributed primarily to diffraction-driven overlap of neighboring Gaussian channels (Koni et al., 2024). The same work emphasizes that such optical emulation is an analog simulator for moderate $M=\binom{m+n-1}{n}$5, not a scalable quantum computer.

A simulation-oriented realization of the same idea appears in the Virtual Quantum Optics Laboratory. There, each propagating beam in a time bin is represented by a polarization Jones vector, Gaussian states are sampled as complex Gaussian random vectors whose statistics reproduce the Wigner functions of vacuum, thermal, coherent, or entangled Gaussian states, and passive optical components are implemented as $M=\binom{m+n-1}{n}$6 or $M=\binom{m+n-1}{n}$7 matrices acting on these vectors. Beam splitters, wave plates, polarizers, and depolarizers are therefore literal matrix transformations, while non-Gaussianity enters only through a nonlinear threshold-detector model that clicks when the amplitude exceeds a fixed threshold (Cour et al., 2021).

Imperfections are likewise matrix-encoded. The most general passive lossy two-port device is described by a non-unitary $M=\binom{m+n-1}{n}$8 scattering matrix supplemented by noise operators,

$M=\binom{m+n-1}{n}$9

with $U\in U(M)$0. In the lossy asymmetric beam-splitter analysis, loss relaxes the phase constraint present in a unitary beam splitter and thereby permits tunable two-photon interference, including programmable Hong–Ou–Mandel dips and peaks that are not available in the lossless case (Uppu et al., 2016).

5. Tensor-network and high-dimensional computational formulations

A major recent development in matrix quantum optics is the use of tensor networks to represent optical states and operators directly as matrix product states (MPS) and matrix product operators (MPO). In the time-domain treatment of ultrafast quantum nonlinear optics, a broadband pulse in a nonlinear waveguide is discretized into spatial bins, mapped to a Bose–Hubbard-type chain for $U\in U(M)$1 or a coupled FH–SH chain for $U\in U(M)$2, and then evolved with TEBD. The quantum state is encoded as an MPS, local observables are MPOs, and an additional “demultiplexing” algorithm constructs a linear-optical unitary that maps arbitrary temporal supermodes onto local bins, making it possible to extract reduced density matrices and Wigner functions for selected pulse modes (Yanagimoto et al., 2021).

That framework was used to simulate Kerr solitons and $U\in U(M)$3 simultons in the few-photon regime. For the Kerr case, the solitonic supermode develops Wigner-function negativity and purity loss because the nominal soliton mode becomes entangled with orthogonal temporal modes. For the simulton case, a two-mode FH–SH reduced density matrix reveals a basis in which one hybrid mode carries most of the non-Gaussianity while the other remains close to coherent, as diagnosed by entanglement negativity and single-mode Wigner functions (Yanagimoto et al., 2021).

A continuous-variable tensor-network formulation beyond the Fock basis pushes this program further by discretizing quadrature space on a grid $U\in U(M)$4, encoding the wavefunction as an MPS in the binary digits of the grid index, and representing $U\in U(M)$5, $U\in U(M)$6, $U\in U(M)$7, $U\in U(M)$8, and the SPDC Hamiltonian as MPOs. Time stepping is performed through an implicit Euler scheme solved variationally in MPS form. For degenerate SPDC with high-intensity pump fields $U\in U(M)$9, the method reports compression ratios above $n$0 while reproducing energy conservation, pump depletion, and quadrature squeezing (Kapridov et al., 19 Nov 2025).

For passive interferometers, an operator-basis MPS formalism addresses Boson-Sampling amplitudes directly from the interferometer matrix $n$1. Instead of evolving a Fock-space state, it encodes Heisenberg-evolved creation operators as an MPS-like object in an operator basis. The resulting graphical contraction rules compute permanents with complexity $n$2, matching Ryser’s algorithm, and extend naturally to partial distinguishability and photon loss through direct-sum tensor constructions and an additional loss mode (Cilluffo et al., 3 Feb 2025).

6. Scope, extensions, and distinct usages

The phrase “matrix quantum optics” is not restricted to passive interferometers or Gaussian optics. In a distinct large-$n$3 usage, it denotes driven-dissipative matrix quantum mechanical models whose planar limit reproduces phenomena familiar from Kerr resonators. The basic degrees of freedom are adjoint matrix oscillators $n$4, $n$5, the Lindblad jumps are $n$6, and the Hamiltonians are single-trace analogues of coherently or parametrically driven Kerr cavities. In this setting, the planar ungauged one-matrix models show evidence of nonequilibrium phase transitions analogous to those reported for driven-dissipative Kerr resonators, including transitions between unique and multiple steady states and $n$7-breaking behavior in the two-quantum pumped model (Cho, 6 Aug 2025).

That large-$n$8 usage also clarifies a conceptual boundary. In gauged matrix quantum mechanics, Lindblad dissipation is absent at large $n$9 because large-$E=L^{(\mathrm{out})} U L^{(\mathrm{in})}$0 factorization cancels the dissipative contributions for singlet jump operators; nontrivial planar open-system dynamics therefore require ungauged models with non-singlet jumps. This is a different construction from mode-space linear optics, but it preserves the central matrix-quantum-optics motif: optical nonequilibrium physics is reformulated in terms of matrix-valued dynamical objects, here at the level of matrix quantum mechanics rather than finite-mode photonics (Cho, 6 Aug 2025).

A neighboring semiclassical strand extends matrix optics into nonlinear polaritonic photonics. The transfer-matrix method for stratified media is generalized to inhomogeneous and nonlinear systems by representing the in-plane structure in a plane-wave basis and treating the nonlinear layer through a convolution operator in $E=L^{(\mathrm{out})} U L^{(\mathrm{in})}$1-space. For exciton-polariton microcavities this yields a nonlinear transfer matrix for a quantum well with intensity-dependent detuning, and the implementation is reported to scale as $E=L^{(\mathrm{out})} U L^{(\mathrm{in})}$2, in contrast with $E=L^{(\mathrm{out})} U L^{(\mathrm{in})}$3 behavior for PLaSK and related approaches (Sajnok et al., 10 Feb 2025). Because the treatment is monochromatic and mean-field, this is not a quantized-field formalism, but it is closely aligned with matrix-quantum-optics design methods for quantum photonic technologies.

In the cited literature, therefore, matrix quantum optics denotes several closely related but nonidentical constructions. What unifies them is not a single formal definition, but a common strategy: optical states, propagators, devices, or observables are organized as structured matrices whose algebra reflects the physics of translation, interference, Gaussian transformation, tomography, nonlinear propagation, tensor-network compression, or nonequilibrium steady-state dynamics (Ornigotti et al., 2012, Laing et al., 2012, Houde et al., 2024, Cilluffo et al., 3 Feb 2025, Cho, 6 Aug 2025).