Single-Particle Density Matrix Approach

• The single-particle density matrix approach is a formulation that represents many-body quantum systems via one-body reduced density matrices, capturing key observables like densities and occupation spectra.
• Natural orbital decomposition within this framework distinguishes between condensed and fragmented states, offering insights into quantum coherence and many-body correlations.
• Advanced formulations link SPDM to time-dependent kinetic equations, variational methods, and modern numerical techniques including stochastic estimation and machine learning.

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)$ (Sakmann et al., 2015). 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 (Blöchl et al., 2013). 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)$ (Gerasimenko, 2020). 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')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$0

where $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$1 are the natural occupations and $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$2 the natural orbitals (Sakmann et al., 2015).

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)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$3 of order $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$4, whereas a fragmented BEC has more than one eigenvalue $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$5 of order $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$6 (Sakmann et al., 2015). 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^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$7

with natural orbitals $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$8 and occupations $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$9. For fermions, the $\rho(\mathbf r)=\rho^{(1)}(\mathbf r|\mathbf r)$0-representability condition requires that $\rho(\mathbf r)=\rho^{(1)}(\mathbf r|\mathbf r)$1 be Hermitian and that its eigenvalues lie in $\rho(\mathbf r)=\rho^{(1)}(\mathbf r|\mathbf r)$2 (Blöchl et al., 2013). In spin-unpolarized closed-shell settings used in single-particle-exact density functional constructions, the effective one-body operator satisfies $\rho(\mathbf r)=\rho^{(1)}(\mathbf r|\mathbf r)$3, $\rho(\mathbf r)=\rho^{(1)}(\mathbf r|\mathbf r)$4, and $\rho(\mathbf r)=\rho^{(1)}(\mathbf r|\mathbf r)$5 (Trappe et al., 2023).

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 $\rho(\mathbf r)=\rho^{(1)}(\mathbf r|\mathbf r)$6-body probability structure.

2. From many-particle density operators to kinetic equations

In Gerasimenko’s operator-hierarchy formulation, observables are represented by sequences $\rho(\mathbf r)=\rho^{(1)}(\mathbf r|\mathbf r)$7, states by sequences of density operators $\rho(\mathbf r)=\rho^{(1)}(\mathbf r|\mathbf r)$8, and expectation values by

$\rho(\mathbf r)=\rho^{(1)}(\mathbf r|\mathbf r)$9

up to the normalization factor $$\rho_{a,b}:=\sum_j P_j\langle\Phi_j|\hat c_b^\dagger \hat c_a|\Phi_j\rangle,$$0 (Gerasimenko, 2020). The exact $$\rho_{a,b}:=\sum_j P_j\langle\Phi_j|\hat c_b^\dagger \hat c_a|\Phi_j\rangle,$$1-particle density operators evolve by the von Neumann equation,

$$\rho_{a,b}:=\sum_j P_j\langle\Phi_j|\hat c_b^\dagger \hat c_a|\Phi_j\rangle,$$2

with generator

$$\rho_{a,b}:=\sum_j P_j\langle\Phi_j|\hat c_b^\dagger \hat c_a|\Phi_j\rangle,$$3

The same framework rewrites the density operator sequence in terms of correlation operators $$\rho_{a,b}:=\sum_j P_j\langle\Phi_j|\hat c_b^\dagger \hat c_a|\Phi_j\rangle,$$4 through cluster expansions. Reduced density operators are then obtained by partial tracing,

$$\rho_{a,b}:=\sum_j P_j\langle\Phi_j|\hat c_b^\dagger \hat c_a|\Phi_j\rangle,$$5

and satisfy the BBGKY hierarchy

$$\rho_{a,b}:=\sum_j P_j\langle\Phi_j|\hat c_b^\dagger \hat c_a|\Phi_j\rangle,$$6

The specifically single-particle closure appears under the initial “chaos” state

$$\rho_{a,b}:=\sum_j P_j\langle\Phi_j|\hat c_b^\dagger \hat c_a|\Phi_j\rangle,$$7

for which the entire state evolution can be expressed functionally in terms of the one-particle reduced correlation operator $$\rho_{a,b}:=\sum_j P_j\langle\Phi_j|\hat c_b^\dagger \hat c_a|\Phi_j\rangle,$$8. Higher reduced correlations become functionals

$$\rho_{a,b}:=\sum_j P_j\langle\Phi_j|\hat c_b^\dagger \hat c_a|\Phi_j\rangle,$$9

and $G_1(t)$0 satisfies the generalized quantum kinetic equation

$G_1(t)$1

with $G_1(t)$2 (Gerasimenko, 2020).

In the mean-field scaling limit, higher-order reduced correlation operators vanish,

$G_1(t)$3

which is the propagation of chaos. The one-particle operator converges to $G_1(t)$4 satisfying the quantum Vlasov equation

$G_1(t)$5

For pure states, this reduces to the Hartree equation and, in related representations, to nonlinear Schrödinger or Gross–Pitaevskii equations (Gerasimenko, 2020).

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 $G_1(t)$6-particle probability density

$G_1(t)$7

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 (Sakmann et al., 2015).

The same work gives a first-principles sampling procedure for single shots by factorizing the $G_1(t)$8-body probability into conditional probabilities,

$G_1(t)$9

and updating a reduced many-body state via repeated application of the field operator. In the MCTDHB setting,

$N$0

with

$N$1

the method generates not only single-shot configurations but also full counting distributions and correlation functions of any order (Sakmann et al., 2015).

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$2 bosons. Upon delocalising the superposition in momentum space, the condensate fraction reaches maximum values larger than $N$3 in the range of particles studied (Ayet et al., 2015).

In optical lattices, the off-diagonal SPDM element between distant sites can be reconstructed from site occupations after two quenches. The effective Hamiltonians are

$N$4

and in the one-particle sector the coherence is

$N$5

Because the $N$6 and $N$7 sectors contribute zero to $N$8, the full SPDM element is

$N$9

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 (Ardila et al., 2018).

For interacting two-band fermionic systems, the SPDM can also encode topology. In the spin-$$|\Psi\rangle = \sum_{\vec n} C_{\vec n}(t)\,|\vec n\rangle,$$0 Haldane model with repulsive on-site interaction, the topological Hamiltonian

$$|\Psi\rangle = \sum_{\vec n} C_{\vec n}(t)\,|\vec n\rangle,$$1

has, under small quasiparticle and quasihole linewidths, the same eigenvectors as $$|\Psi\rangle = \sum_{\vec n} C_{\vec n}(t)\,|\vec n\rangle,$$2, 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 $$|\Psi\rangle = \sum_{\vec n} C_{\vec n}(t)\,|\vec n\rangle,$$3 component at the Dirac point $$|\Psi\rangle = \sum_{\vec n} C_{\vec n}(t)\,|\vec n\rangle,$$4 (Zheng et al., 2018).

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 $$|\Psi\rangle = \sum_{\vec n} C_{\vec n}(t)\,|\vec n\rangle,$$5-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 $$|\Psi\rangle = \sum_{\vec n} C_{\vec n}(t)\,|\vec n\rangle,$$6 sits within the BBGKY hierarchy and depends explicitly on the two-body reduced density matrix $$|\Psi\rangle = \sum_{\vec n} C_{\vec n}(t)\,|\vec n\rangle,$$7: $$|\Psi\rangle = \sum_{\vec n} C_{\vec n}(t)\,|\vec n\rangle,$$8 The practical obstacle is the interaction term. The reviewed time-dependent density-matrix propagation literature identifies energy conservation, positivity, natural-orbital occupation bounds, and $$|\Psi\rangle = \sum_{\vec n} C_{\vec n}(t)\,|\vec n\rangle,$$9-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 $$\hat\Psi(\mathbf r)=\sum_j \hat b_j\,\phi_j(\mathbf r),$$0-representability and can violate positivity conditions (Elliott et al., 2016).

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 $$\hat\Psi(\mathbf r)=\sum_j \hat b_j\,\phi_j(\mathbf r),$$1 make the generator depend on $$\hat\Psi(\mathbf r)=\sum_j \hat b_j\,\phi_j(\mathbf r),$$2, yet the paper proves that the eigenvalues of the single-particle density matrix remain in $$\hat\Psi(\mathbf r)=\sum_j \hat b_j\,\phi_j(\mathbf r),$$3. 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 (Rosati et al., 2014).

In finite-temperature correlated dynamics, density-matrix coupled-cluster theory rewrites the time-dependent density matrix in thermo-field form as

$$\hat\Psi(\mathbf r)=\sum_j \hat b_j\,\phi_j(\mathbf r),$$4

evolves it along the Keldysh contour, and preserves the trace exactly: $$\hat\Psi(\mathbf r)=\sum_j \hat b_j\,\phi_j(\mathbf r),$$5 In the single-impurity Anderson model, DMCC-S is exact for the noninteracting case $$\hat\Psi(\mathbf r)=\sum_j \hat b_j\,\phi_j(\mathbf r),$$6, whereas DMCC-SD substantially improves the interacting dynamics and becomes nearly indistinguishable from the renormalization-group benchmark (Shushkov et al., 2019).

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 $$\hat\Psi(\mathbf r)=\sum_j \hat b_j\,\phi_j(\mathbf r),$$7 and $$\hat\Psi(\mathbf r)=\sum_j \hat b_j\,\phi_j(\mathbf r),$$8, and the density matrix obeys a Lindblad master equation with a spontaneous-emission dissipator. The diagonal entries $$\hat\Psi(\mathbf r)=\sum_j \hat b_j\,\phi_j(\mathbf r),$$9 and $$\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.$$0 are populations, the off-diagonal element $$\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.$$1 is the bunching/coherence variable, and the model is reliable in the practical operating regime $$\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.$$2 (Fares et al., 2018).

For Bose-Einstein condensates, an extended matrix formalism beyond the usual Nambu $$\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.$$3 language unifies the single-particle propagator, two-particle response, and vertex structures. In the low-energy and low-momentum limit at $$\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.$$4, 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 $$\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.$$5 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 (Watabe, 2020).

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,

$$\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.$$6

where $$\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.$$7 is universal in the sense that it depends on the interaction $$\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.$$8 but not on the external one-particle Hamiltonian $$\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.$$9. 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 (Blöchl et al., 2013).

A broader second-quantized formulation generalizes the density variable itself. For any choice of single-particle modes $\rho_{ij}$0, the basic density is the occupation-number list

$\rho_{ij}$1

This includes configuration-space density $\rho_{ij}$2, momentum-space density $\rho_{ij}$3, and occupancies of the eigenstates of the one-body Hamiltonian. In the latter basis, the energy functional takes the explicitly single-particle-exact form

$\rho_{ij}$4

so that all one-body contributions are treated exactly and only the interaction functional remains to be approximated (Englert et al., 2022).

Single-particle-exact density functional theory develops this idea into a practical variational scheme. The single-particle Hamiltonian

$\rho_{ij}$5

defines the 1pEx basis

$\rho_{ij}$6

and the variational variables are the participation numbers

$\rho_{ij}$7

with, for an unpolarized even-$\rho_{ij}$8 system,

$\rho_{ij}$9

The exact one-particle contribution is then

$$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$00

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^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$01, 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 (Trappe et al., 2023).

These functionals establish the SPDM as a variational object rather than only a reduced observable. The formal bridge to Green’s functions (Blöchl et al., 2013), the basis-independent occupation-number formulation (Englert et al., 2022), and the explicit participation-number optimization strategy (Trappe et al., 2023) 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 $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$02, the defining constraints are

$$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$03

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 (Truflandier et al., 2020).

For large sparse Hamiltonians, stochastic estimation of the density matrix exploits the identity

$$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$04

where $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$05 is the single-particle density matrix and $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$06 the free-energy matrix function. Gradient-based probing estimates $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$07 stochastically and differentiates it, rather than probing $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$08 directly. In zero-temperature metals, the stochastic error for local density-matrix elements scales as

$$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$09

while the convergence becomes exponential for finite-temperature or insulating systems (Wang et al., 2017).

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 $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$10 a.u., the learned density matrices often reduce SCF convergence to one iteration for $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$11 and $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$12 and to two iterations for $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$13, and the resulting speedup is about $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$14 to $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$15 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 (Hazra et al., 2024).

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 $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$16-expansion. In the cited application to dipolar Hartree–Fock theory, the second-order correction produces a finite gradient term of the form

$$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$17

with negative and small coefficient $$\rho^{(1)}(\mathbf r|\mathbf r')=\langle \Psi|\hat\Psi^\dagger(\mathbf r')\hat\Psi(\mathbf r)|\Psi\rangle,$$18 (Bencheikh et al., 2016).

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.