---
title: Single-Particle Density Matrix Approach
url: https://www.emergentmind.com/topics/single-particle-density-matrix-approach
type: topic
---

# Single-Particle Density Matrix Approach

Searching arXiv for recent and foundational papers on the single-particle density matrix approach.
arXiv search query: "single-particle density matrix approach"
The single-particle density matrix approach denotes a class of many-body formulations in which a quantum system is represented, propagated, measured, or approximated through its one-body reduced density matrix rather than through the full \(N\)-particle wave function or density operator. In a bosonic field-theoretic representation, the first-order reduced density matrix is
\[
\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,
\]
and its diagonal gives the single-particle density \(\rho(\mathbf r)=\rho^{(1)}(\mathbf r|\mathbf r)\) [1501.03224]. In a grand-canonical fermionic formulation, the one-particle reduced density matrix is
\[
\rho_{a,b}:=\sum_j P_j\langle\Phi_j|\hat c_b^\dagger \hat c_a|\Phi_j\rangle,
\]
which can be diagonalized into natural orbitals and occupations [1309.2406]. In operator-hierarchy approaches, the full many-particle density operator can be reduced, via cluster expansions and marginalization, to a typical-particle description governed by a generalized quantum kinetic equation for a one-particle operator \(G_1(t)\) [2001.01180]. Taken together, these constructions place the single-particle density matrix at the interface between microscopic unitary dynamics, reduced kinetic theory, observable one-body coherence, and computational electronic structure.

## 1. Fundamental object, natural orbitals, and representability

For an \(N\)-boson state expanded as
\[
|\Psi\rangle = \sum_{\vec n} C_{\vec n}(t)\,|\vec n\rangle,
\]
with bosonic field operator
\[
\hat\Psi(\mathbf r)=\sum_j \hat b_j\,\phi_j(\mathbf r),
\]
the first-order reduced density matrix may be written in the orbital basis as
\[
\rho^{(1)}(\mathbf r|\mathbf r')=\sum_{i,j}\rho_{ij}\,\phi_i^*(\mathbf r')\,\phi_j(\mathbf r), \qquad \rho_{ij}=\langle \Psi|\hat b_i^\dagger \hat b_j|\Psi\rangle.
\]
Diagonalizing \(\rho_{ij}\) yields the natural-orbital decomposition
\[
\rho^{(1)}(\mathbf r|\mathbf r')=\sum_i \rho_i\,\phi_i^{NO}(\mathbf r)\phi_i^{NO*}(\mathbf r'),
\]
where \(\rho_i\) are the natural occupations and \(\phi_i^{NO}\) the natural orbitals [1501.03224].

This decomposition is the standard spectral organization of one-body coherence. In the bosonic examples emphasized in the literature, it is also the diagnostic for condensation and fragmentation: a condensed BEC has exactly one eigenvalue \(\rho_1\) of order \(N\), whereas a fragmented BEC has more than one eigenvalue \(\rho_i\) of order \(N\) [1501.03224]. The same work explicitly links fragmentation to correlation: correlated states are at least partially fragmented, and fragmentation is visible through the occupation spectrum of the first-order reduced density matrix.

In fermionic reduced-density-matrix theory, diagonalization gives
\[
\rho=\sum_n f_n |\psi_n\rangle\langle\psi_n|,
\]
with natural orbitals \(|\psi_n\rangle\) and occupations \(f_n\). For fermions, the \(N\)-representability condition requires that \(\rho\) be Hermitian and that its eigenvalues lie in \([0,1]\) [1309.2406]. In spin-unpolarized closed-shell settings used in single-particle-exact density functional constructions, the effective one-body operator satisfies \(\rho\ge 0\), \(\operatorname{tr}\rho=N\), and \(\rho\le 2\) [2305.03233].

The central significance of the single-particle density matrix is therefore twofold. First, it is the reduced object from which one-body densities, currents, and occupation spectra are obtained directly. Second, its spectral data encode the extent to which a many-body state departs from a single-determinant or fully condensed structure. A recurrent limitation, stressed repeatedly in the literature, is that this object captures average one-body properties but not the full \(N\)-body probability structure.

## 2. From many-particle density operators to kinetic equations

In Gerasimenko’s operator-hierarchy formulation, observables are represented by sequences \(A=(I,A_1,A_2,\dots)\), states by sequences of density operators \(D=(I,D_1,D_2,\dots)\), and expectation values by
\[
(A,D)=\sum_{n=0}^{\infty}\frac{1}{n!}\operatorname{Tr}_{1,\dots,n} A_n D_n,
\]
up to the normalization factor \((I,D)\) [2001.01180]. The exact \(n\)-particle density operators evolve by the von Neumann equation,
\[
D_n(t)=G_n^*(t)D_n^0=e^{-itH_n}D_n^0\,e^{itH_n},
\]
with generator
\[
N_n^* f_n = -i(H_n f_n-f_n H_n).
\]

The same framework rewrites the density operator sequence in terms of correlation operators \(g(t)=(I,g_1(t),g_2(t),\dots)\) through cluster expansions. Reduced density operators are then obtained by partial tracing,
\[
F_s(t,1,\dots,s)=\frac{1}{(I,D)}\sum_{n=0}^{\infty}\frac{1}{n!}\operatorname{Tr}_{s+1,\dots,s+n}D_{s+n}(t,1,\dots,s+n),
\]
and satisfy the BBGKY hierarchy
\[
\frac{\partial}{\partial t}F_s(t,1,\dots,s) = N_s^*F_s(t,1,\dots,s) + \sum_{j=1}^s \operatorname{Tr}_{s+1}N^*_{\mathrm{int}(j,s+1)}\,F_{s+1}(t,1,\dots,s,s+1).
\]

The specifically single-particle closure appears under the initial “chaos” state
\[
G^{(c)}=(0,G_1,0,\dots),
\]
for which the entire state evolution can be expressed functionally in terms of the one-particle reduced correlation operator \(G_1(t)\). Higher reduced correlations become functionals
\[
G_s(t,1,\dots,s\,|\,G_1(t)),\qquad s\ge 2,
\]
and \(G_1(t)\) satisfies the generalized quantum kinetic equation
\[
\frac{\partial}{\partial t}G_1(t,1) = N^*(1)G_1(t,1) + \operatorname{Tr}_2\,N^*_{\mathrm{int}(1,2)}\,G_1(t,1)G_1(t,2) + \operatorname{Tr}_2\,N^*_{\mathrm{int}(1,2)}\,G_2(t,1,2\,|\,G_1(t)),
\]
with \(G_1(0)=G_1^0\) [2001.01180].

In the mean-field scaling limit, higher-order reduced correlation operators vanish,
\[
\lim_{\varepsilon\to 0}\varepsilon^s G_s(t)=0,\qquad s\ge 2,
\]
which is the propagation of chaos. The one-particle operator converges to \(g_1(t)\) satisfying the quantum Vlasov equation
\[
\frac{\partial}{\partial t}g_1(t,1) = N^*(1)g_1(t,1) + \operatorname{Tr}_2\,N^*_{\mathrm{int}(1,2)}\,g_1(t,1)g_1(t,2).
\]
For pure states, this reduces to the Hartree equation and, in related representations, to nonlinear Schrödinger or Gross–Pitaevskii equations [2001.01180].

This line of work gives the single-particle density matrix approach a rigorous hierarchy-theoretic interpretation. The one-particle object is not introduced as an ad hoc closure; it is obtained as the exact reduced descriptor after cluster expansion, marginalization, and, where appropriate, mean-field scaling.

## 3. Coherence, single-shot physics, and experimental reconstruction

A central distinction in ultracold-boson applications is the difference between the reduced one-body description and actual measurement outcomes. The single-particle density gives the probability of finding a particle at a given position after averaging over many realizations, whereas a single experimental shot is governed by the full \(N\)-particle probability density
\[
P(\mathbf r_1,\dots,\mathbf r_N)=|\Psi(\mathbf r_1,\dots,\mathbf r_N)|^2.
\]
Accordingly, interference fringes, fluctuating vortices, and broad center-of-mass distributions can appear in single shots even when they are absent or weak in the average single-particle density [1501.03224].

The same work gives a first-principles sampling procedure for single shots by factorizing the \(N\)-body probability into conditional probabilities,
\[
P(\mathbf r_1,\dots,\mathbf r_N) = P(\mathbf r_1)\,P(\mathbf r_2|\mathbf r_1)\cdots P(\mathbf r_N|\mathbf r_{N-1},\dots,\mathbf r_1),
\]
and updating a reduced many-body state via repeated application of the field operator. In the MCTDHB setting,
\[
i\frac{\partial}{\partial t}|\Psi\rangle=\hat H|\Psi\rangle,
\]
with
\[
\hat H=\sum_{i=1}^N\left(-\frac{1}{2}\frac{\partial^2}{\partial \mathbf r_i^2}+V(\mathbf r_i)\right) +\lambda_0\sum_{i<j}\delta_\epsilon(\mathbf r_i-\mathbf r_j),
\]
the method generates not only single-shot configurations but also full counting distributions and correlation functions of any order [1501.03224].

Direct evaluation of off-diagonal SPDMs is also possible in integrable many-body settings. For quantum bright solitons constructed from superpositions of attractive Lieb–Liniger string states, a coordinate-space diagrammatic method computes the SPDM and its eigenvalues directly from Bethe-ansatz wave functions, with modest numerical resources, for systems up to \(N=10\) bosons. Upon delocalising the superposition in momentum space, the condensate fraction reaches maximum values larger than \(97\%\) in the range of particles studied [1510.04311].

In optical lattices, the off-diagonal SPDM element between distant sites can be reconstructed from site occupations after two quenches. The effective Hamiltonians are
\[
H_{ab}^z = B_z (n_a - n_b), \qquad H_{ab}^x = B_x(c_a^\dagger c_b + c_b^\dagger c_a),
\]
and in the one-particle sector the coherence is
\[
\chi_{ab}^{(1)} = \mathrm{tr}(\rho_{ab} \, c_a^\dagger c_b) = \frac{r_x + i r_y}{2}.
\]
Because the \(n=0\) and \(n=2\) sectors contribute zero to \(c_a^\dagger c_b\), the full SPDM element is
\[
\chi_{ab}=p_1\chi_{ab}^{(1)}.
\]
The scheme applies to fermions and hard-core bosons, relies on engineered distant tunneling and site-resolved occupation measurements, and does not generalize to soft-core bosons [1806.08171].

For interacting two-band fermionic systems, the SPDM can also encode topology. In the spin-\(\tfrac12\) Haldane model with repulsive on-site interaction, the topological Hamiltonian
\[
H_t({\bf k}) \equiv -G^{-1}({\bf k}, i\omega=0)
\]
has, under small quasiparticle and quasihole linewidths, the same eigenvectors as \(\rho_{\bf k}^T\), enabling reconstruction of Berry curvature and first Chern number from SPDM tomography. The transition point is identified by SPDM-gap closing and by the sign change of the \(a_{{\bf k},3}\) component at the Dirac point \(K\) [1812.01991].

These developments make explicit that the SPDM is experimentally accessible, but only through protocols matched to the observable of interest. A common misconception is that a density snapshot directly reveals the full SPDM; the optical-lattice and topological protocols show that off-diagonal coherence requires controlled rotations, quenches, or engineered couplings, while single-shot bosonic structure requires the full \(N\)-body probability distribution rather than the reduced one-body density alone.

## 4. Propagation, dissipation, and excitation structure

For many-electron dynamics, the exact equation of motion for the one-body reduced density matrix \(\rho_1\) sits within the BBGKY hierarchy and depends explicitly on the two-body reduced density matrix \(\rho_2\):
\[
i\dot{\rho_1}(\mathbf r',\mathbf r,t) = \Big(-\nabla^2/2 + v_{\rm ext}(\mathbf r,t) + \nabla'^2/2 - v_{\rm ext}(\mathbf r',t)\Big)\rho_1(\mathbf r',\mathbf r,t) + \int d^3r_2\, f_{ee}(\mathbf r,\mathbf r',\mathbf r_2)\, \rho_2(\mathbf r',\mathbf r_2,\mathbf r,\mathbf r_2,t).
\]
The practical obstacle is the interaction term. The reviewed time-dependent density-matrix propagation literature identifies energy conservation, positivity, natural-orbital occupation bounds, and \(N\)-representability as the key exact constraints. Adiabatic density-matrix approximations typically keep the occupation numbers fixed and miss double excitations, whereas the Frozen-Gaussian-assisted TDDMFG scheme captures changing occupation numbers and doubly excited structure but does not guarantee \(N\)-representability and can violate positivity conditions [1602.03723].

A distinct route begins from a many-body Lindblad master equation and derives a closed single-particle density-matrix equation by mean-field factorization. The resulting equation is nonlinear and generally non-Lindblad because the Pauli-blocking factors \((\delta-\rho)\) make the generator depend on \(\rho\), yet the paper proves that the eigenvalues of the single-particle density matrix remain in \([0,1]\). This positivity-preserving single-particle formulation is contrasted with conventional non-Lindblad Markov approaches, whose mean-field equations can lead to positivity violations and unphysical results. The construction is further extended to open quantum systems with spatial boundaries, yielding a Lindblad-like system-reservoir scattering superoperator [1408.1898].

In finite-temperature correlated dynamics, density-matrix coupled-cluster theory rewrites the time-dependent density matrix in thermo-field form as
\[
\ket{\rho_t} = e^{T_t}\ket{\rho_0},
\]
evolves it along the Keldysh contour, and preserves the trace exactly:
\[
\bra{1}\ket{\rho_t}=\bra{1}\ket{\rho_0}.
\]
In the single-impurity Anderson model, DMCC-S is exact for the noninteracting case \(U=0\), whereas DMCC-SD substantially improves the interacting dynamics and becomes nearly indistinguishable from the renormalization-group benchmark [1907.11962].

Single-particle density matrices also support compact open-system models. In the quantum free-electron laser, the electron is approximated as a two-level system with momentum states \(|0\rangle\) and \(|-1\rangle\), and the density matrix obeys a Lindblad master equation with a spontaneous-emission dissipator. The diagonal entries \(P_0\) and \(P_{-1}\) are populations, the off-diagonal element \(B\) is the bunching/coherence variable, and the model is reliable in the practical operating regime \(D\le 0.07\) [1803.01941].

For Bose-Einstein condensates, an extended matrix formalism beyond the usual Nambu \(2\times2\) language unifies the single-particle propagator, two-particle response, and vertex structures. In the low-energy and low-momentum limit at \(T=0\), the single-particle Green’s function and the density response share the same phonon pole, consistent with the Gavoret–Nozières correspondence, while the Nepomnyashchii identity \(\Sigma_{12}(0)=0\) controls the infrared structure. At nonzero temperature, random-phase-approximation calculations show that both the single-particle spectral function and the density response can exhibit satellite structures due to beyond-mean-field effects [2003.10813].

The technical lesson across these developments is that propagation of a single-particle density matrix is not merely an economical truncation. Its viability depends on preserving trace, positivity, and representability, and on handling the residual two-body information either through exact hierarchy functionals, controlled mean-field reductions, or explicitly constructed correlated closures.

## 5. Variational and functional formulations

Reduced-density-matrix functional theory casts equilibrium many-body physics as a constrained search over the one-particle density matrix. In grand-canonical form,
\[
\Omega_{\beta,\mu}(\hat h+\hat W) = \min_\rho \left\{ [\rho(h-\mu 1)] + F_\beta^{\hat W}[\rho] \right\},
\]
where \(F_\beta^{\hat W}[\rho]\) is universal in the sense that it depends on the interaction \(\hat W\) but not on the external one-particle Hamiltonian \(\hat h\). A central result is that the exact density-matrix functional can be derived from the Luttinger–Ward functional of the Green’s function and is convex, so the grand potential obeys a true minimum principle rather than a merely stationary one [1309.2406].

A broader second-quantized formulation generalizes the density variable itself. For any choice of single-particle modes \(|a\rangle\), the basic density is the occupation-number list
\[
n_a=\operatorname{tr}\bigl(\psi^\dagger(a)\psi(a)\rho\bigr), \qquad \sum_a n_a = N.
\]
This includes configuration-space density \(n(\mathbf r)\), momentum-space density \(n(\mathbf p)\), and occupancies of the eigenstates of the one-body Hamiltonian. In the latter basis, the energy functional takes the explicitly single-particle-exact form
\[
E[n]=\sum_{k=0}^\infty \varepsilon_k n_k + E_{\text{pair}}[n],
\]
so that all one-body contributions are treated exactly and only the interaction functional remains to be approximated [2206.10097].

Single-particle-exact density functional theory develops this idea into a practical variational scheme. The single-particle Hamiltonian
\[
H_{\mathrm{1p}}(\mathbf P,\mathbf R) = \frac{1}{2m}\mathbf P^2 + V_{\mathrm{ext}}(\mathbf R)
\]
defines the 1pEx basis
\[
H_{\mathrm{1p}}|a\rangle = E_a |a\rangle,
\]
and the variational variables are the participation numbers
\[
n_a = \langle a|\rho|a\rangle,
\]
with, for an unpolarized even-\(N\) system,
\[
0 \le n_a \le 2, \qquad \sum_a n_a = N.
\]
The exact one-particle contribution is then
\[
E_{\mathrm{1p}}=\sum_a n_a E_a.
\]
The interaction functional is approximated through density-matrix constructions based on two schemes: an iterative “matrix mixer” that enforces the closed-shell representability condition \(\rho^2=2\rho\), and a Thomas–Fermi-inspired phase-space ansatz. Proof-of-principle simulations on interacting Fermi gases and on atoms and ions, with and without relativistic corrections, are reported to be typically accurate at the one-percent level [2305.03233].

These functionals establish the SPDM as a variational object rather than only a reduced observable. The formal bridge to Green’s functions [1309.2406], the basis-independent occupation-number formulation [2206.10097], and the explicit participation-number optimization strategy [2305.03233] all pursue the same objective: to preserve the exact single-particle structure while parameterizing the genuinely many-body remainder through controlled functionals.

## 6. Numerical estimators, perturbative algorithms, and learned density matrices

Density matrix perturbation theory replaces the sum-over-states structure of Rayleigh–Schrödinger perturbation theory by direct equations for perturbed one-particle density matrices. For the unperturbed density matrix \(D\), the defining constraints are
\[
[D,F]=0,\qquad D=D^\dagger,\qquad D^2=D,\qquad \mathrm{tr}(D)=N.
\]
The benchmarked variants comprise SOS-McWeeny DMPT, Sylvester-DMPT, and recursive purification DMPT extended to hole-particle canonical purification. HPCP-DMPT is reported to have stable convergence profiles but, for a given perturbation order, to require roughly three times as many matrix multiplications as TC2 [2007.04739].

For large sparse Hamiltonians, stochastic estimation of the density matrix exploits the identity
\[
f(H)=\frac{d\Omega}{dH^T},
\]
where \(f(H)\) is the single-particle density matrix and \(\Omega=g(H)\) the free-energy matrix function. Gradient-based probing estimates \(g(H)\) stochastically and differentiates it, rather than probing \(f(H)\) directly. In zero-temperature metals, the stochastic error for local density-matrix elements scales as
\[
\Delta f \sim S^{-(d+2)/2d},
\]
while the convergence becomes exponential for finite-temperature or insulating systems [1711.10570].

Machine learning has been used to predict the ground-state one-particle density matrix of Kohn–Sham DFT directly from atomic positions. In the reported models, the test error is of order \(\mathrm{MAE}\sim 2\times 10^{-4}\) a.u., the learned density matrices often reduce SCF convergence to one iteration for \(\mathrm{H_2O}\) and \(\mathrm{S_2O}\) and to two iterations for \([\mathrm{Fe(H_2O)_6}]^{2+}\), and the resulting speedup is about \(3\times\) to \(5\times\) compared with standard guesses. The predicted density matrices are also accurate enough for force evaluation and accelerated ab initio molecular dynamics with little or no self-consistent iteration [2401.06533].

Semiclassical constructions provide a different numerical route. For a two-dimensional Fermi gas, the Grammaticos–Voros Wigner-transform expansion yields a one-particle density matrix approximation that preserves Hermiticity and idempotency to all orders in the \(\hbar\)-expansion. In the cited application to dipolar Hartree–Fock theory, the second-order correction produces a finite gradient term of the form
\[
\Delta E_{dd}^{(1)} = \varepsilon \int d^2R\,\frac{(\nabla\rho)^2}{\sqrt{\rho}},
\]
with negative and small coefficient \(\varepsilon\) [1605.08014].

These algorithmic developments clarify the present computational status of the single-particle density matrix approach. It is simultaneously a target of perturbative solvers, a sparse object for stochastic estimation, a learned surrogate for self-consistent-field initialization, and a semiclassical carrier of exact algebraic constraints. The literature also delineates the associated caveats: adiabatic closures can freeze occupations, stochastic gains depend on spatial decay, and data-driven density matrices must still respect electron-number and basis constraints.

Source: https://www.emergentmind.com/topics/single-particle-density-matrix-approach