Percus–Yevick: mathematics and assumptions#

Percus–Yevick (PY) is a closure approximation for the correlations of a homogeneous, isotropic equilibrium hard-sphere mixture. It produces functions of separation and wavenumber rather than particle centres or a deposition history. See Analytical structure factors for the Python workflow, Metropolis Monte Carlo: mathematics and assumptions for explicit equilibrium sampling, and RSA Monte Carlo: mathematics and assumptions for irreversible deposition.

Physical inputs and correlation functions#

Species \(i\) has radius \(a_i>0\), number fraction \(x_i\), and number density \(\rho_i\), with \(\sum_i x_i=1\). A number fraction counts particles; a volume fraction weights their occupied volume. For total volume fraction \(\eta\),

\[\rho=\frac{\eta}{\sum_i x_i(4\pi a_i^3/3)},\qquad \rho_i=x_i\rho,\qquad \eta_i=\rho_i\frac{4\pi a_i^3}{3},\qquad \sum_i\eta_i=\eta.\]

Additive spheres have contact distance \(\sigma_{ij}=a_i+a_j\). Their potential is infinite below contact and zero outside. The equilibrium fluid assumptions exclude walls, directional order, attractions, and irreversible insertion. Mathematical validity of an input volume fraction does not establish a stable fluid phase at that density.

The partial pair correlation \(g_{ij}(r)\) gives the relative likelihood of finding a species \(j\) centre at separation \(r\) from species \(i\):

\[\mathrm{d}N_j=\rho_jg_{ij}(r)4\pi r^2\mathrm{d}r,\qquad h_{ij}(r)=g_{ij}(r)-1.\]

\(h\) is the total correlation: it includes both direct and indirect correlation through other particles. In a homogeneous fluid, \(g\) tends to one far away; the hard core requires \(g=0\) below contact.

Ornstein–Zernike equation and the PY postulate#

The direct correlation \(c_{ij}\) separates the direct part from chains mediated by other species. The mixture Ornstein–Zernike equation is

\[h_{ij}(\mathbf{r})=c_{ij}(\mathbf{r})+ \sum_\ell\rho_\ell\int_{\mathbb{R}^3} c_{i\ell}(\mathbf{r}-\mathbf{r}')h_{\ell j}(\mathbf{r}')\, \mathrm{d}^3\mathbf{r}'.\]

This relation alone contains two unknown functions. A closure supplies another relation; PY postulates

\[c_{ij}(r)=\bigl(1+h_{ij}(r)-c_{ij}(r)\bigr) \bigl(\exp[-\beta u_{ij}(r)]-1\bigr).\]

For hard spheres this is equivalent to

\[h_{ij}(r)=-1\quad(r<\sigma_{ij}),\qquad c_{ij}(r)=0\quad(r>\sigma_{ij}).\]

The first condition enforces impenetrability; the second is the approximation. The total correlation outside contact generally remains nonzero because the Ornstein–Zernike convolution propagates indirect correlations. A closure is not a deposition or particle-motion rule.

Fourier transforms and PackLab’s matrix normalization#

Wavenumber \(k\) has units of inverse length. For an isotropic function, PackLab’s three-dimensional radial transform convention is

\[\widetilde{h}_{ij}(k)=4\pi\int_0^\infty r^2h_{ij}(r)\operatorname{sinc}(kr)\,\mathrm{d}r,\qquad \operatorname{sinc}(z)=\frac{\sin z}{z},\quad\operatorname{sinc}(0)=1.\]

Here \(\widetilde h\) has units of volume. PackLab returns the density-normalized, dimensionless matrices

\[C_{ij}(k)=\sqrt{\rho_i\rho_j}\,\widetilde c_{ij}(k),\qquad H_{ij}(k)=\sqrt{\rho_i\rho_j}\,\widetilde h_{ij}(k).\]

The result field Cpy is \(C\), and H is \(H\), not the unnormalized transform \(\widetilde h\). Multiplying the Fourier-space Ornstein–Zernike equation by the density factors gives

\[H=C+CH,\qquad (I-C)H=C,\qquad S=I+H=(I-C)^{-1}.\]

\(I\) is the identity matrix and \(S\) is the corresponding partial structure-factor matrix. It describes density fluctuations in reciprocal space. H and Cpy have shape (species, species, wavenumber); g and h use distance as their last axis. Form \(S\) from H by adding the species identity at each wavenumber. An experimentally weighted mixture scattering intensity also needs species amplitudes and their weights; it is not simply every element of \(S\) summed together.

The inverse transform for positive-density species is

\[h_{ij}(r)=\frac{1}{2\pi^2\sqrt{\rho_i\rho_j}} \int_0^\infty k^2H_{ij}(k)\operatorname{sinc}(kr)\,\mathrm{d}k.\]

PackLab evaluates this numerically on a finite wavenumber grid. It explicitly enforces the hard core in the reported g; h retains the numerical inverse transform, so g = 1 + h can fail inside the core. For pairs involving a zero-density species, the inverse transform sets h to zero and g is zero inside contact and one outside; this convention does not compute infinite-dilution tracer correlations.

Resolution, periodic terminology, and interpretation#

Use wavenumber="auto" for a first calculation. For a chosen real-space resolution \(\delta r\) and largest requested distance \(r_{\max}\), a useful grid design is

\[k_{\max}=\frac{\pi}{\delta r},\qquad \Delta k\lesssim\frac{2\pi}{m r_{\max}},\]

where \(m\) is the number of samples per oscillation of the inverse kernel (default 12). Increasing \(k_{\max}\) resolves sharper features; decreasing \(\Delta k\) resolves oscillations over longer distances. These are independent requirements. Finite cutoff causes ringing near the contact discontinuity; coarse spacing causes integration errors. Retain PackLab’s inadequate-sampling warning and test convergence by increasing the cutoff and refining the spacing separately. A smooth-looking curve alone is not an accuracy test. Compare distances away from contact as well as the contact region, where limits from either side differ.

Periodic images and minimum-image distances apply to the explicit simulation box (see RSA Monte Carlo: mathematics and assumptions). The PY solver assumes an infinite bulk mixture; its wavenumber grid is an integration grid, not a list of periodic box modes. It has no wrapped centres, winding numbers, or finite particle count. Match radii, densities, and normalization when comparing with Metropolis, and check equilibration and box size on the simulation side. RSA need not agree with PY because its deposition process is irreversible. Grid refinement reduces numerical error within PY; it does not remove the physical error of the closure approximation.

Percus–Yevick: mathematics and assumptions#

Physical inputs and correlation functions#

Ornstein–Zernike equation and the PY postulate#

Fourier transforms and PackLab’s matrix normalization#

Resolution, periodic terminology, and interpretation#

Further reading#