---
abstract: 'We present a theory to predict the structure and thermodynamics of mixtures of colloids of different diameters, building on our earlier work \[J. Chem. Phys. 145, 074904 (2016)\] that considered mixtures with all particles constrained to have the same size. The patchy, solvent particles have short-range directional interactions, while the solute particles have short-range isotropic interactions. The hard-sphere mixture without any association site forms the reference fluid. An important ingredient within the multi-body association theory is the description of clustering of the reference solvent around the reference solute. Here we account for the physical, multi-body clusters of the reference solvent around the reference solute in terms of occupancy statistics in a defined observation volume. These occupancy probabilities are obtained from enhanced sampling simulations, but we also present statistical mechanical models to estimate these probabilities with limited simulation data. Relative to an approach that describes only up to three-body correlations in the reference, incorporating the complete reference information better predicts the bonding state and thermodynamics of the physical solute for a wide range of system conditions. Importantly, analysis of the residual chemical potential of the infinitely dilute solute from molecular simulation and theory shows that whereas the chemical potential is somewhat insensitive to the description of the structure of the reference fluid the energetic and entropic contributions are not, with the results from the complete reference approach being in better agreement with particle simulations.'
author:
- Artee Bansal
- Arjun Valiya Parambathu
- 'D. Asthagiri'
- 'Kenneth R. Cox'
- 'Walter G. Chapman'
title: 'Thermodynamics of mixtures of patchy and spherical colloids of different sizes: a multi-body association theory with complete reference fluid information'
---

[^1]

Introduction
============

In thermodynamic perturbation theory of association involving short range interactions between molecules, the properties of the reference fluid plays a central role. In the typical situation when the reference is a hard-sphere fluid, perturbation theories usually use information about two body, and at times three body, correlations in the reference fluid to describe the physical (associating) system. For example, Wertheim’s theory [@wertheim_fluids_1984; @wertheim_fluids_1984-1] and its extensions based on the statistical associating fluid theory (SAFT) [@chapman_new_1990] use pair correlation information at contact to estimate extent of association between pairs of molecules. In SAFT, for the hard-sphere reference either the Carnahan-Starling [@carnahan_equation_1969] equation for a single component fluid or the Boublik-Mansoori-Carnahan-Starling-Leland [@boublik_hardsphere_1970; @mansoori_equilibrium_1971] equation for a mixture are used to describe the pair-correlation information at contact. The structure of Wertheim’s theory or SAFT is such that one can obtain accurate extent of association and thermodynamics even for systems with strong inter-particle interactions provided the representation of the reference is adequate. However, as the complexity of the interaction increases in the physical system, such as may result from multiple bonding and size asymmetries, information about two or three body correlations in the reference no longer suffices.

In our previous work [@bansal_structure_2016], we studied the multi-body correlation functions of a symmetric hard sphere reference fluid in terms of the probabilities of observing $n$ molecules in the bonding region. These occupancy probabilities were obtained from enhanced sampling Monte Carlo simulations for the hard sphere fluid. We developed a procedure to use this information within the Marshall-Chapman formalism [@marshall_molecular_2013; @marshall_thermodynamic_2013] to describe multiple association of solvent molecules to a solute molecule. This *complete reference* approach proved successful in predicting the bonding state and thermodynamics of a colloidal solute in a patchy solvent for a wide range of system conditions [@bansal_structure_2016].

Here we study mixtures where the solute diameter is as small as half to as large as twice the diameter of the solvent. The solvent particles are spheres with directional interaction sites and the solute particles are spheres with isotropic interactions, and the solute is capable of bonding with multiple solvent particles. The structure and thermodynamics of mixture of hard spheres with different diameters has been studied in detail before [@torquato_microstructure_1986; @reiss_statistical_1959; @mayer_integral_1947; @torquato_microstructure_1982; @torquato_microstructure_1983; @torquato_microstructure_1985], but a compact form for the correlations beyond the contact value is still unavailable. Further, for systems with large size asymmetries even the pair-correlation information obtained using the Boublik-Mansoori-Carnahan-Starling-Leland equation is inadequate [@feng_contact_2011]. Our approach of including multi-body correlations rests on using the occupancy statistics [@reiss_upper_1981; @pratt_quasichemical_2001; @pratt_selfconsistent_2003; @bansal_structure_2016] of the hard-sphere solvent around the hard-sphere solute. We find that representing multi-body correlation functions in terms of occupancy statistics in physically reasonable observation volumes accurately captures the multi-bonding effects in asymmetric mixtures. These occupancy statistics are obtained from particle simulations. Importantly, we also present a physically transparent, statistical mechanical model to describe the occupancy probabilities in symmetric and asymmetric hard sphere fluids for different packing fractions. This model corrects multi-body effects obtained for isolated clusters by incorporating the role of the cluster-bulk interface and the bulk medium effects. We also investigate the energy-entropy decomposition of the chemical potential of the solute in a model system with only solute-solvent interactions to better appreciate the role of the reference fluid. Throughout, theoretical results are validated versus molecular simulations

The rest of the paper is organized in the following way. In Section \[sc:bentheory\] we discuss the association potential of the system and describe how packing effects are important for the given potential. The Marshall-Chapman [@marshall_thermodynamic_2013] theory is briefly introduced to show the multi-density representation of the free energy, and based on our previous work [@bansal_structure_2016], an improved representation of multi-body correlations(*complete reference*) is presented. In Section \[sc:HStheory\], we examine hard sphere packing around a reference particle and develop models based on statistical mechanics [@reiss_upper_1981] and hard sphere simulation data for different densities. We apply our complete reference approach for different asymmetric mixtures of solute and solvent and present results in Section \[sc:res\_asso\]. In Section \[sc:res\_hs\] we present results for the hard sphere reference system (symmetric and asymmetric mixtures) based on the correlation developed in section \[sc:HStheory\]. We also provide simulation results for isolated cluster probabilities in asymmetric hard sphere mixtures in the appendix (Section \[sc:appen\]).

Theory
======

Asymmetric mixtures with different association geometries  {#sc:bentheory}
----------------------------------------------------------

The focus of our study is asymmetric mixtures containing molecules with short range attractive interactions. The short range association potential is the same as that in previous work [@bansal_structure_2016]: the solute molecule can associate with multiple solvent molecules isotropically and the patchy solvent has directional interactions. The total potential is a sum of hard sphere and association contributions $$u{(r)}=u_{HS}{(r)}+u_{AS}{(r)}
\label{eq:potT}$$

The association potential for patchy-patchy $(p,p)$ and spherical-patchy $(s,p)$ particles is: $$u_{AB}^{(p,p)}{(r)}=
    \begin{cases}
        -\epsilon_{AB}^{(p,p)}, r<r_c \,\text{and}\, \theta_A\leq \theta_c^{(A)}\,\text{and}\,\theta_B \leq \theta_c^{(B)}
        \\
        0   \text{ \ \ \ \ \ otherwise}
        \\    
    \end{cases}
    \label{eq:potential1}$$

$$u_A^{(s,p)}{(r)}=
    \begin{cases}
        -\epsilon_A^{(s,p)}, r<r_c\, \text{and}\, \theta_A\leq \theta_c^{(A)}
        \\
        0   \text{ \ \ \ \ \ otherwise}
        \\    
    \end{cases}
    \label{eq:potential2}$$

where the subscripts $A$ and $B$ represent the type of site and $\epsilon$ is the association energy; $r$ is the distance between the particles; and $\theta_A$ is the angle between the vector connecting the centers of two molecules and the vector connecting association site $A$ to the center of that molecule (Fig. \[fig:1\]). The critical distance beyond which particles do not interact is $r_c$ and $\theta_c$ is the solid angle beyond which sites cannot bond. Fig. \[fig:1\] shows examples of solute-solvent and solvent-solvent short range interaction geometries for different sizes of solute particles.

![Association between solute and solvent (a) and solvent molecules (b). Different Cases with solute larger (middle) and smaller (left) than solvent molecules are studied. $r$ is the center-to-center distance and $\theta_A$ and $\theta_B$ are the orientation of the attractive patches $A$ and $B$ relative to line connecting the centers. The solute (colored red) can only interact with patch $A$ (colored red).[]{data-label="fig:1"}](figure1)

Since the solute can associate with multiple solvent molecules (Eq. \[eq:potential2\]), it is important to study the multi-body correlations that determine the packing of solvent particles around the solute in the reference fluid [@bansal_structure_2016]. The difficulty in determining these interactions arises due to the limited knowledge in describing multi-body correlation functions for $n\ge 3$. But the volume integral of the multi-body correlation has a clear physical meaning in terms of average number of $n$-solvent clusters ($F^{(n)}$, Fig. \[fig:Fn\]). In particular, for the distinguished solute, $$\begin{aligned}
{F^{(n)}} & = &\frac{{{\rho_p^n}}}{{n!}}\int\limits_{v} {d{{\vec r}_1} \cdots \int\limits_{v} d{{\vec r}_n}{g_{HS}}\left( {{{\vec r}_1} \cdots {{\vec r}_n}|0} \right)} \nonumber \\
& = & {\sum\limits_{m = n}^{{n^{\max }}} {C^m_n p_m}} \, ,
\label{eq:Fn}\end{aligned}$$ where $\rho_p$ is the density of solvent particles, $p_n$ is the probability of observing exactly $n$ solvent particles in the observation volume of the solute ($v$) defined by the spherical region of radius $r_c$, $C^m_n$ ($=m! / (m-n)!\cdot n!$), and $g_{HS}({\vec r}_1 \cdots {\vec r}_n|0)$ is the distribution function of the $n$-solvent particles around the solute at the center of the observation volume, indicated by $(\ldots | 0)$. $n^{max}$ is the maximum number of solvent molecules that can occupy the observation volume around the reference solute.

![image](figure2)

In Wertheim’s multi-density formalism [@wertheim_fluids_1986; @wertheim_fluids_1986-1], the free energy due to association ($A^{AS}$) is expressed as $$\frac{{A^{AS}}}{{V{k_{\rm B}}T}} = \sum {\left( {{\rho_k}\ln \frac{{\rho ^{(0)}_{k }}}{{{\rho_{ k }}}} + {Q^{\left( k \right)}} + {\rho _{ k }}} \right)} -   \frac{  \Delta c^{(0)}}{V}
%  - {{\Delta {c^{(0)}}} \mathord{\left/
%       {\vphantom {{\Delta {c^{(0)}}} V}} \right.
%       \kern-\nulldelimiterspace} V}
  \label{eq:3}$$

where $k_{\rm B}$ is the Boltzmann constant, $T$ is the temperature, the summation is over the species ($k={s,p}$), $\rho$ is the number density, $\rho^{(0)}$ is the monomer density, $Q^{(k)}$ is obtained from the Marshall-Chapman development [@marshall_thermodynamic_2013] and $\Delta{c^{0}}$ is the contribution to the graph sum due to association between the solvent-solvent $(p,p)$ and solute-solvent $(s,p)$ molecules, i.e. $$\Delta c_{}^{\left( 0 \right)} = \Delta c_{pp}^{\left( 0 \right)} + \Delta c_{sp}^{\left( 0 \right)}
  \label{eq:4}$$

Marshall and Chapman [@marshall_molecular_2013; @marshall_thermodynamic_2013] extended Wertheim’s theory beyond the single bonding condition to incorporate multi-body effects in a solution consisting of an isotropic solute and solvent with directional interactions. The contribution to free energy due to association between solute and solvent molecules was obtained as $$\frac {\Delta c_{sp}^{(0)}}{V} = \sum\limits_{n = 1}^{{n^{\max }}} {\frac{\Delta c_n^{( 0)}}{V} }
    \label{eq:5}$$ where the sum is over different coordination states of the solute and $\Delta c_n^{(0)}$ is given by:

$$\begin{aligned}
        \Delta c_n^{(0)} & = & \frac{\rho ^{(0)}_s {( {\rho_p X_A^{(p)}})}^n} {\Omega^{n + 1} n!} \int d(1)\cdots d(n + 1) \,
        g_{HS}( 1 \cdots n + 1) \cdot \prod\limits_{k = 2}^{n + 1} {( f_{A}^{(s,p)} ( 1,k))}  \, .
        \label{eq:14}
    \end{aligned}$$

In Eq. \[eq:14\], $\rho_p=\rho\cdot x^{(p)}$ is the density of solvent molecules obtained from the mole fraction of solvent($x^{(p)}$) and the total density($\rho$), $X_A^{(p)}$ is the fraction of solvent molecules not bonded at site A, $\Omega =4\pi$ is the total number of orientations, $f_{A}^{( {s,p})}(1,k) = (\exp (\varepsilon_A^{(s,p)}/k_BT) - 1)$ is the Mayer function for association between $p$ and $s$ molecules corresponding to potential in Eq. \[eq:potential2\] and the integral is over all the orientations and positions of the $n+1$ particles. By taking the average association strength and acceptable orientations out of the integral and fixing the solute at the origin, the above integral can be rewritten as

$$\begin{aligned}
        \frac{  \Delta c_n^{(0)}}{V} & = & \frac{\rho _s^{(0)} {( {{\rho}_{p}X_A^{(p)}}f_A^{(s,p)}\sqrt {\kappa _{AA}})}^n} {n!} \int_{v} d\vec r_1 \cdots  \int_{v} d\vec r_n \,
        g_{HS}(\vec r_1 \cdots \vec r_n |0)  \, .
        \label{eq:3}
    \end{aligned}$$

Marshall and Chapman[@marshall_molecular_2013; @marshall_thermodynamic_2013] approximated the integral in Eq. \[eq:3\] as $$\begin{aligned}
  \int_{v} d\vec r_1 \cdots \int_{v} d\vec r_n\, g_{HS}(\vec r_1 \cdots \vec r_n |0)  \approx y_{HS}^n( \sigma) \delta ^{(n)} \Xi ^{(n)} \, ,
  \label{eq:MCA}
  \end{aligned}$$ where $\Xi^{(n)}$ is the partition function for an isolated cluster of $n$ solvent hard-spheres around a solute hard-sphere, $y_{HS}( \sigma)$ is (pair) cavity correlation function at contact, and $\delta^{\left(n\right)}$ corrects the superposition of cavity correlation functions for three body interactions. We will hereafter refer to Eq. \[eq:MCA\] as the Marshall-Chapman approximation (MCA).

As shown earlier [@bansal_structure_2016], MCA fails for high densities and high association energies, conditions where multi-body interactions are important. But recognizing that the integral in Eq. \[eq:3\] is related to $F^{(n)}$ (Eq. \[eq:Fn\]) we have [@bansal_structure_2016] $$\begin{aligned}
\frac{  \Delta c_n^{(0)}}{V}  = {\rho_s^{(0)} {( {x ^{(p)}X_A^{(p)}}f_A^{(s,p)}\sqrt {\kappa _{AA}})}^n} F^{(n)}  \, .
\label{eq:Cn_new}\end{aligned}$$ It can be observed that all the multi-body correlation information is subsumed in $F^{(n)}$ which is obtained from the occupancy distribution $\{p_n\}$. We follow our earlier work [@bansal_structure_2016] to estimate this distribution. Importantly, since $\{p_n\}$ forms the basis of our *complete reference* approach, we also develop an analytical model to describe these distribution functions.

Finally, with the above information, and based on the Marshall-Chapman theory [@marshall_thermodynamic_2013], the fraction of solute associated with $n$ solvent molecules is $$X_n^{\left( s \right)} = \frac{ {( {x ^{(p)}X_A^{(p)}}f_A^{(s,p)}\sqrt {\kappa _{AA}})}^n F^{(n)} }{{1 + \sum\limits_{n = 1}^{{n^{\max }}} {( {x ^{(p)}X_A^{(p)}}f_A^{(s,p)}\sqrt {\kappa _{AA}})}^n} F^{(n)} } \, ,
\label{eq:301}$$ and the fraction of solute not bonded to any solvent molecule is $$X_0^{\left( s \right)} = \frac{1}{{1 + \sum\limits_{n = 1}^{{n^{\max }}} {( {x ^{(p)}X_A^{(p)}}f_A^{(s,p)}\sqrt {\kappa _{AA}})}^n} F^{(n)} } \, .
\label{eq:300}$$ Using these distributions for associating mixture, the average number of solvent associated with the solute is given by: $$n_{avg} = \sum\limits_n {n\cdot {X^{(s)}_n}}  \, ,
\label{eq:81}$$ The fraction of solvent not bonded at site $A$ and site $B$ can be obtained by simultaneous solution of the following equations: $$X_A^{\left( p \right)} = \frac{1}{{1 + \xi {\kappa _{AB}}f_{AB}^{\left( {p,p} \right)}{\rho_p}X_B^{(p)} + \frac{{{\rho_s}}}{{{\rho_p}}}\frac{{ n_{avg} }}{{X_A^{(p)}}}}}  \, ,$$ $$X_B^{\left( p \right)} = \frac{1}{{1 + \xi {\kappa _{AB}}f_{AB}^{\left( {p,p} \right)}{\rho_p}X_A^{(p)}}} \, .
\label{eq:82}$$ where $$\begin{aligned}
    \xi  & = & 4\pi {\sigma^2}\left( {{r_c} - \sigma} \right){y_{HS}}(\sigma) \\ 
    \kappa_{AB} & = &\left[1-cos(\theta_c)\right]^2/{4} \\
    f_{AB}^{({p,p})} & = & \exp ( \varepsilon _{AB}^{({p,p})}/k_{\rm B}T)-1 \, .\end{aligned}$$

Occupancy distribution $\{p_n\}$ for the hard-sphere fluid {#sc:HStheory}
----------------------------------------------------------

Consider a hard sphere fluid with one solute and $N$ solvent particles in a volume $V$ and temperature $T$. We are interested in the occupancy statistics $\{p_n\}$ of the solvent in the coordination volume around the solute. To this end consider the reaction $$S{P_{n = 0}} + {P_n} \rightleftharpoons S{P_n} \, ,
\label{eq:20}$$ with the equilibrium constant $${K_n} = \frac{{{\rho _{S{P_n}}}}}{{{\rho _{S{P_{n = 0}}}}\rho _p^n}} \, ,
\label{eq:21}$$ where $\rho_{SP_n}$ is the density of species $SP_n$ and $\rho_p$ is the density of the solvent. Clearly, we have [@pratt_quasichemical_2001; @lrp:book; @lrp:cpms] $$\frac{p_n}{p_0} = K_n \rho_p^n \, ,
\label{eq:pnp0}$$ where $p_0$ is the probability that the coordination volume is empty of solvent particles. Following earlier work in studying clusters with quasichemical theory [@pratt_quasichemical_2001; @pratt_selfconsistent_2003; @merchant_thermodynamically_2009; @merchant:jcp11b], we can show that $$K_n = \frac{(e^{ \beta\mu^{\rm ex}_{p} })^n}{n!}\int\limits_v d{\vec r}_1\ldots\int\limits_v d{\vec r}_n e^{-\beta U_{SP_n}(R^n)} e^{-\beta \phi(R^n;\beta)} \, , 
\label{eq:kn1}$$ where $U_{SP_n}(R^n)$ is the potential energy of the $n$-solvent cluster (with the solute $S$ fixed at the center of the cluster), $\beta = (k_{\rm B}T)^{-1}$, $\phi(R^n;\beta)$ is the free energy of interaction of the cluster with the rest of the bulk medium for a given configuration $R^n$ of the cluster, and $\beta\mu^{\rm ex}_{p} $ is the excess chemical potential of the solvent particle. (For completeness, in appendix A we derive the above expression for $K_n$.) $\phi(R^n;\beta)$ can also be thought of as a field imposed by the bulk solvent medium [@pratt_selfconsistent_2003; @merchant:jcp11b] on the solute-solvent cluster in the observation volume. Earlier Pratt and Ashbaugh [@pratt_selfconsistent_2003] modeled this field using a self-consistent approach. Here we take a different approach.

First note that without the field term, the cluster integral presents a simpler $n$-body problem (where $n$ is small, typically less than 20 for systems of interest here). The field is thus an interfacial term that couples the local cluster with the bulk medium. To make this explicit, we can rewrite Eq. \[eq:kn1\] as $$K_n = \frac{(e^{ \beta\mu^{\rm ex}_{p} })^n}{n!} \langle e^{-\beta \phi(R^n;\beta)}\rangle_0 \int\limits_vd{\vec r}_1\ldots\int\limits_v d{\vec r}_n e^{-\beta U_{SP_n}(R^n)} \, .$$ Here $\langle \ldots\rangle_0$ indicates averaging with over the normalized probability density for cluster conformations $R^n$ in the absence of interactions with the rest of the medium, i.e. over the density $e^{-\beta U_{SP_n}(R^n)}  / (n! K_n^{(0)})$, where $$n! K_n^{(0)} =  \int\limits_vd{\vec r}_1\ldots\int\limits_v d{\vec r}_n e^{-\beta U_{SP_n}(R^n)} \, , 
\label{eq:kn0}$$ and the interfacial contribution is $\beta\Omega_n = -\ln \langle e^{-\beta \phi(R^n;\beta)}\rangle_0$.

From analysis of simulation data for different densities, we find that $\Omega_n$ can be described by a two parameter equation as $$\beta \Omega_n =-\zeta_1 \cdot  n^2+\zeta_2\cdot n \, .$$ This model of the interfacial term was anticipated in our previous work [@bansal_structure_2016]. Thus we finally obtain $${p_n} = \frac{\exp(\beta \cdot n \cdot \mu^{ex}_{p})\rho_p^n [{\exp{(\zeta_1 \cdot  n^2-\zeta_2\cdot n)}}]{K_n}^{(0)}}{{1 + \sum\limits_{j \ge 1} {\exp(\beta \cdot j \cdot \mu^{ex}_{p})\rho_p^j [{\exp{(\zeta_1 \cdot j^2-\zeta_2 \cdot j)}}]{K_j}^{(0)}} }} \ 
\label{eq:pn_2par}$$ The above equation can also be derived using a two moment maximum entropy approach, with the mean and variance of the occupancy as constraints and $K_n^{(0)}$ as the default (see appendix A).

Drawing upon the work of Reiss and Merry [@reiss_upper_1981], we can model the interfacial term in terms of surface sites (of the cluster) that are available to interact with the bulk fluid. On the basis of such a mean field approach and guided by Monte Carlo (MC) simulation data for different densities of hard sphere systems, we find that $$\beta \Omega_n =(-0.0109 \cdot  n^2+1.0109 \cdot n)\cdot \zeta
\label{eq:surf_int}$$ $${p_n} = \frac{\exp(\beta \cdot n \cdot \mu^{ex}_{p})\rho_p^n [{\exp{(-(-0.0109 \cdot  n^2+1.0109 \cdot n)\cdot \zeta)}}]{K_n}^{(0)}}{{1 + \sum\limits_{j \ge 1} {\exp(\beta \cdot j \cdot \mu^{ex}_{p})\rho_p^j [{\exp{(-(-0.0109 \cdot  j^2+1.0109 \cdot j)\cdot \zeta)}}]{K_j}^{(0)}} }} \ 
\label{eq:pn_1par}$$

Eq. \[eq:pn\_2par\] and Eq. \[eq:pn\_1par\] are the 2-parameter and 1-parameter models, respectively, for $p_n$ on the basis of which we obtain $F^{(n)}$ (Eq. \[eq:Fn\]) to describe multi-body correlations in the reference fluid. The parameter values for different densities are given in Table \[table: param1\]. These parameters were obtained based on hard spheres mixtures with all particles of the same size. Since the information about size-asymmetry is already contained in the isolated cluster partition function, we will use these parameters to study asymmetric mixtures and will discuss limitations for cases with extreme size ratio.

METHODS
=======

To compare the theory results, we perform Monte-Carlo (MC) simulations for a range of systems. This section presents the details of the MC simulations for different associating and hard sphere systems.The Marshall-Chapman approximation (MCA) and the models developed for hard sphere distribution functions require isolated cluster probabilities; these were also computed for different size ratios.

Monte Carlo Simulations
-----------------------

MC simulations were carried out for reference hard sphere systems and associating systems to compare the results of Marshall-Chapman theory using MCA and the complete reference approach. The associating mixture contains the patchy solvent particles and the isotropically interacting solute defined by the potentials given by Eq. \[eq:potential1\] and Eq. \[eq:potential2\], respectively. The solute diameter is $\sigma_s$ and solvent diameter is $\sigma_p$.

The observation volume is defined by a critical radius $r_c = 1.1\bar{\sigma}$, where $\bar{\sigma} = (\sigma_s + \sigma_p)/2$ is the closest distance of approach. For cases where $\sigma_s/\sigma_p \geq 1.5$, a cutoff of 1.1$\bar{\sigma}$ can include some of the second-shell solvent. To avoid this and focus attention to the first observation shell, we set $r_c = \bar{\sigma} + 0.1\sigma_p$ for these cases. In the associating system (Fig. \[fig:1\] and Eqs. \[eq:potential1\] and \[eq:potential2\]), the critical angles for interaction are $\theta_c^{(A)} = \theta_c^{(B)} = 27\degree$.

For hard sphere mixtures 255 solvent particles and 1 solute particle were studied in a given simulation cell. Ensemble reweighting technique was used to reveal low probability states [@merchant_water_2011]. The system was equilibrated for 1 million steps with translational factors chosen to yield an acceptance rate of 0.3, and data was collected every 100 sweeps. Analysis was carried out for different densities and size ratios $\sigma_s/\sigma_p$.

For associating mixtures bonding distributions and average bonding numbers were studied for mixtures with different sizes and different association strengths for solute-solvent and solvent-solvent interactions. The excess chemical potential of the coupling of the colloid with solvent was also calculated using thermodynamic integration of average binding energy of solute with solvent as a function of solute-solvent interactions, using the three-point Gauss Legendre quadrature technique [@Hummer:jcp96]. For a symmetric mixture with no solvent-solvent interactions, energetic and entropic contributions for solute chemical potential were also studied at constant volume and temperature.

Concentration effects were also computed considering a total of 864 particles, with varied number of solute particles. Due to the difference in size of the solute and solvent, the computations were performed keeping the packing fraction constant. Hence, the density of the system changed with respect to the concentration. Also the maximum angle for which the patch can form single bond is computed using the law of cosines, and hence the critical angle needs to be altered when the solute size is smaller than the solvent. For a size ratio of $\sigma_s/\sigma_p=0.8$, critical angle $\theta_c^{(A)} = \theta_c^{(B)} = 20\degree$ was used to ensure single bonding condition for the $A$ patch on solvent molecules which can associate with solute molecules.

Isolated cluster probabilities
------------------------------

For asymmetric mixtures, we also study isolated cluster probabilities for different size ratios of solute and solvent molecules. The observation volume around the isolated solute is the same as defined in the previous section. For different size ratios we use the spherical code (appendix B) to estimate the maximum number of solvent molecules which can be inserted in the observation volume of the isolated solute. Successive insertion probabilities are calculated as in previous works [@marshall_molecular_2013; @bansal_structure_2016], where $10^8$ to $10^9$ insertions were carried out in the observation volume around the solute and the cases with no overlap with remaining $n-1$ particles studied. The data for isolated cluster probabilities for different sizes is given in the appendix B.

Results and discussions
=======================

Hard sphere $\{p_n\}$ distribution {#sc:res_hs}
----------------------------------

Recall that Eq. \[eq:pn\_1par\] and Eq. \[eq:pn\_2par\] are the 1-parameter and 2-parameter models for the occupancy distribution $\{p_n\}$. To obtain the parameters for hard spheres all of the same size, for both the models we use the average occupancy ($n^{HS}_{avg}= \sum\limits_n n\cdot {p_n}$) as a fitting constraint. For the 2-parameter model we additionally use the exclusion probability ($p_0$) — the probability when no hard sphere solvent particle is present in the observation volume — as a constraint. For the 1-parameter model we study the surface interactions based on the mean field approach developed by Reiss and Merry [@reiss_upper_1981]. By analyzing the distribution functions $\{p_n\}$ for different densities, we obtain geometric effects (density independent) that describe the mutual interference of different surface sites (Eq. \[eq:surf\_int\]). The density (or packing fraction) dependent parameters for these two models are given in Table  \[table: param1\].

   $\rho \sigma^3$   $\eta$   $\zeta_1$   $\zeta_2$   $\zeta$
  ----------------- -------- ----------- ----------- ---------
         0.2         0.105     0.0175       0.773      0.75
         0.3         0.157      0.015       1.267      1.251
         0.4         0.209     0.0256       1.979      1.947
         0.5         0.262      0.023       2.88       2.876
         0.6         0.314     0.0361       4.179      4.172
         0.7         0.367     0.0457       5.947      5.985
         0.8         0.419     0.0609       8.432      8.562
         0.9         0.471     0.0829      12.088     12.414

  : Parameters for Eq. \[eq:pn\_2par\] ($\zeta_1$,$\zeta_2$) and Eq. \[eq:pn\_1par\] ($\zeta$) for different densities $\rho \sigma^3$ and packing fractions ($\eta$)

\[table: param1\]

Fig. \[fig:HS\_sym\] presents the results corresponding to the average occupancy ($n^{HS}_{avg}$) and exclusion probability ($\ln p_0$) based on the models. We compare these results with the MC simulation values and also include results from literature [@chang_real_1994; @torquato_microstructure_1985].

![image](figure3)

Fig. \[fig:HS\_sym\] makes it clear that 2-parameter model can simultaneously capture both $n_{avg}$ and $\ln p_0$ in excellent agreement with simulation. Importantly, even the 1-parameter model is able to capture most of the details, affirming the physical ideas underlying the models (Eqs. \[eq:pn\_1par\] and  \[eq:pn\_2par\]).

Next, using the parameters obtained above for a fluid where both solute and solvent hard-spheres are the same size (a symmetric mixture), we describe the occupancy in a fluid where the solute and solvent are of different sizes (an asymmetric mixture). Our *ansatz* is that the information about size asymmetry is adequately captured by the isolated cluster partition function $K_n^{(0)}$ (Eq. \[eq:kn0\]). For an infinitely dilute system comprising one solute in a solvent bath, Fig. \[fig:pn1\_asym\] shows the predictions of $n_{avg}^{HS}$ and $\ln p_0$ based on the 1-parameter model (Eq. \[eq:pn\_1par\]) for different size ratios and different reduced densities. The level of agreement is encouraging, but perhaps not surprising since the size ratios are not much different from a symmetric case. Thus the geometric effects describing the mutual interference of different surface sites in the packing around the solute will be similar to what is observed in the symmetric mixture. For extreme size ratios this should break down, as we will discuss in the context of average bonding numbers in associating mixtures.

![image](figure4)

For asymmetric mixtures, the packing fraction is a better measure of packing in the fluid. Fig. \[fig:pn1\_asym\_eta\]

![image](figure5)

presents the variation of average occupancy for different packing fractions for three different size ratios for an infinitely dilute solution. The results show that the 1-parameter model is able to predict the average occupancy quite well. We note that as the concentration of solute is changed in asymmetric mixtures, the packing fraction changes (for a given density) and parameters corresponding to the resulting packing fraction should be used from Table \[table: param1\].

Associating mixture {#sc:res_asso}
-------------------

We next consider associating fluids and investigate both size and concentration effects.

### Infinite Dilution

We first study an infinitely dilute solution and vary the size of solute with respect to a fixed size of solvent particles. In our complete reference theory, the reference fluid $\{p_n\}$ distribution is either computed directly from simulations (‘$p_n$-Simulation’ in Fig. \[fig:n\_avg\]) or from the 1-parameter model discussed in Sec. A above (‘$p_n$-Model1’ in Fig. \[fig:n\_avg\]). Using $\{p_n\}$ we compute $F^{(n)}$ (Eq. \[eq:Fn\]), and on that basis, the average number of bonds in the associating mixture using Eq. \[eq:Cn\_new\], Eq. \[eq:301\] and Eq. \[eq:81\]. Fig. \[fig:n\_avg\] shows the variation of average bonding numbers with size ratio of solute and solvent molecules for a density of 0.8 and association strength of 7 $k_{\rm B}T$.

![image](figure6)

As the size of the solute increases with respect to the size of the solvent, more solvent molecules can associate with the solute. With accurate information about the reference fluid $\{p_n\}$ (and hence $F^{(n)}$) from direct simulations, the complete reference theory is able to capture this increase in bonding numbers quite accurately. Interestingly, even $\{p_n\}$ obtained using the 1-parameter model (Eq. \[eq:pn\_1par\]) suffices. But note that the Marshall-Chapman approximation with only up to 3-body effects incorporated in the theory underestimates the average bonding numbers. These results emphasize the importance of multi-body interactions in describing the association correctly.

As noted in Sec.  \[sc:HStheory\], the amount of surface exposure (or surface sites) is an important factor in determining the packing effects in the models. For the size ratios considered in Fig. \[fig:n\_avg\], the maximum number of surface sites and hence the geometric effects in surface interactions are expected to be similar to the symmetric case, and not surprisingly, the agreement of bonding numbers between simulation and with the 1-parameter model for $p_n$ is very good. For extreme size ratios (Table. \[table:extreme\]), a high error with the 1-parameter model is expected. This results because of the disparity in surface sites for these ratios relative to the symmetric mixture within which the density independent geometric effects were obtained (Eq. \[eq:surf\_int\]). However, even for these extreme size ratios, with information for reference hard sphere distribution functions from simulation, the complete reference theory is able to capture the average bonding numbers for these extreme size ratios quite well.

   $\sigma_s / \sigma_p $    MC     $p_n$-Simulation   $p_n$-Model1
  ------------------------ ------- ------------------ --------------
            0.5             3.61          3.55             4.50
             2              15.84        15.04            17.66

  : Comparison of average bonding numbers ($n_{avg}$) of solute for extreme size ratios.

\[table:extreme\]

### Varying association strengths

To understand the effect of varying association strengths between solute-solvent and solvent solvent particles, we studied a case with a fixed size ratio of $\sigma_s / \sigma_p = 0.8 $ and for varying association strengths (fig. \[fig:Asso\_asym\]).

![image](figure7)

As the strength of solute-solvent association is increased as compared to solvent-solvent interactions, multi-body effects become important. It was observed that higher deviations are observed with TPT2-based Marshall-Chapman approximation, especially for increasing strength of solute-solvent association. Importantly, excellent agreement with the complete reference theory is observed for all cases noted in the figure.

### Chemical potential and energy entropy contributions

Fig. \[fig:Asso\_mu\] shows the chemical potential of the solute for two limiting cases: one where the solvent-solvent association is present and one where it is absent, for two size ratios.

![image](figure8)

The results indicate that despite the large deviations in $n_{avg}$ (Fig. \[fig:n\_avg\]) noted above, the deviations using the Marshall-Chapman approximation are not as high relative to the complete reference approach.

To better understand this result, we consider a symmetric system where there is no association between solvent particles. For this case, the partial molar energy is the energy of the system with solute-solvent association and there is no contribution for change in solvent-solvent association due to inclusion of an ideal solute at infinite dilution. We decompose the residual chemical potential of the solute ($\mu^{Asso}_{s}$) into its energetic and entropic contributions [@yu_thermodynamic_1988; @yu_solvation_1990]:

$$\beta \mu_s^{Asso} =  \beta E_s^{Asso}-T\cdot \beta S_{s,V}^{Asso}$$

where $ E_s^{Asso} =  ( \frac{\partial \mu^{Asso}_s/T}{\partial 1/T})_{{\rho _s},{\rho _p}}$ is the partial molar energy and $ S_{s,V}^{Asso}=  {( {\frac{{\partial \mu _s^{Asso}}}{{\partial T}}})_{{\rho _s},{\rho _p}}}$ is the partial molar entropy contribution. The entropy and energy values can be obtained with MCA and complete reference approach based on the corresponding temperature derivatives of the chemical potential of solute.

![image](figure9)

Fig. \[fig:en\_entr\] clearly shows that energy and entropy contributions are not captured accurately by MCA. But the apparent reasonable prediction of the chemical potential (using MCA) results because of cancellation of the errors between the entropy and energy contributions. The above result shows that in comparing perturbation theories, it could be useful and prudent to study chemical potential and also its energy and entropy contributions.

### Concentrated systems

We also study the variation of average bonding number of the solute for different concentration of the solute in the asymmetric mixtures. As the concentration of the solute increases, the system becomes limited in the number of solvent molecules that can bond to the solute molecules and hence TPT2 approximation used in MCA becomes more accurate. For low concentrations, multi-body correlations are important and hence deviations are observed with MCA. Fig. \[fig:n\_avg\_conc\] shows the results for two different size ratios. The theory with the complete reference is able to capture the average bonding numbers for the whole concentration range.

![image](figure10)

Concluding discussions
======================

We have studied asymmetric mixtures having strong short-range association between differently sized solute and solvent molecules. The solute molecules have isotropic interactions and the solvent molecules have directional interactions. Such systems are archetypes of colloidal mixtures that are being actively studied in designing materials from the nanoscale level. These systems can also describe the short range ion-solvation and ion-pairing effects in electrolyte systems which is another focus of our research.

The isotropic interactions of the solute can allow multiple solvent molecules to associate and hence multi-body effects become important for these systems. Previously [@bansal_structure_2016], we discussed the development of an accurate perturbation theory for these systems is hindered by the difficulty in obtaining the multi-body correlations in the reference system (typically hard sphere). We discussed the limitation of an approach based on obtaining multi-body correlations for the hard sphere system in the gas phase and approximating bulk solvent effect with a linear superposition of pair correlation (together with term to account for three body correction). It was observed that this second order perturbation method fails at high association strengths and high densities. We introduced an approach to represent the multi-body clusters in terms of occupancy distribution to accurately describe the packing in the hard sphere system. Excellent agreement with MC simulation for a range of conditions of association and concentrations were obtained with this *complete reference* approach for symmetric mixtures.

Here, we have built upon the earlier work and study systems with size asymmetry. Our study shows that the multi-body correlations for asymmetric hard sphere mixtures can be accurately studied in terms of occupancy distributions. With these accurate packing effects, our approach gives excellent agreement with MC simulation for the asymmetric associating mixtures. These occupancy distributions were obtained by particle simulations. Based on ideas borrowed from quasichemical theory, we have also developed parametric models to describe occupancy distributions in the hard sphere systems of different densities and different asymmetries. These distributions were obtained by describing the effects of clustering, medium and surface interactions simultaneously in the hard sphere packing around a solute and can be incorporated in perturbation theories (for eg. statistical associating fluid theory) without having to perform particle simulations for the reference fluid. We validate this complete reference theory (using parameterized models of the reference distribution) against several simulations. A critical test was in analyzing the energy and entropy contributions to the chemical potential of the solute. For a system where solvent-solvent association is not present, such a decomposition of chemical potential in energy and entropy contributions with [complete reference]{} theory showed excellent agreement with MC simulations. The apparently reasonable agreement of the second order perturbation approach is shown to arise from the balancing of errors in the energy and entropy contributions. This important finding suggests the need to study different properties while validating perturbative theories for fluids. This present framework can prove useful in modeling real solutions where concentration of solute is low and its size is different from the solvent molecules.

Acknowledgment
==============

We acknowledge RPSEA / DOE 10121-4204-01 and the Robert A. Welch Foundation (C-1241) for financial support.

Appendix {#sc:appen}
========

Hard sphere distribution
------------------------

### Expression for $K_n$ {#sc:Kn}

From Eq. \[eq:pnp0\], we find that obtaining an expression for $K_n$ reduces to evaluating the ratio of $p_n / p_0$. The total potential energy of the system when $n$ solvent particles are coordinated with the solute and the remaining $N-n$ solvent particles are outside the observation volume can be formally written as $U = U_{SP_n} + U_{N-n|SP_n} + U_{N-n}$. $U_{SP_n}$ is the potential energy of the solute-$n$-solvent cluster. $U_{N-n|SP_n}$ is the interaction energy of the cluster with the rest of the solvent; specifically $U_{N|SP_0}$ is the interaction energy of the solute with the fluid outside the observation volume. In the particular case of hard-spheres, $U_{N|SP_0} = 0$. Finally, $U_{N-n}$ is the potential energy of the solvent constituting the bulk.

Since $p_n \propto Q(n,N-n,V,T)$, where $Q(n,N-n,V,T)$ is the canonical partition function of the system with $n$ solvent in the observation volume around the solute and $N-n$ in the bulk, we immediately have $$\frac{p_n}{p_0} = \frac{N!}{(N-n)!n!} \frac {\int\limits_v d{\vec r}_1\ldots \int\limits_v d{\vec r}_n e^{-\beta U_{SP_n}} \cdot \int\limits_{V-v} d{\vec r}_{n+1} \ldots \int\limits_{V-v} d{\vec r}_{N-n} e^{-\beta U_{N-n|SP_n}} e^{-\beta U_{N-n}}}{\int\limits_{V-v} d{\vec r}_1 \ldots \int\limits_{V-v} d{\vec r}_N e^{-\beta U_{N}}} \, 
\label{eq:pnporatio}$$ where we have implicitly moved the center of the coordinates to the center of the solute and thus canceled a common factor of $V$ from both the numerator and denominator. Further, since $U_{SP_0} = 0$ and $U_{N|SP_0} = 0$, the denominator simply depends on the potential energy of the solvent in the bulk. (Of course for a general solute, this restriction is easily removed [@merchant_thermodynamically_2009; @merchant:jcp11b].)

Next consider the ratio $$\frac{Q(0,N,V-v)}{Q(0,N-n,V-v)} = \frac{Q(0,N-n+1,V-v)}{Q(0,N-n,V-v)}\cdot \frac{Q(0,N-n+2,V-v)}{Q(0,N-n+1,V-v)}\ldots \frac{Q(0,N,V-v)}{Q(0,N-1,V-v)}
\label{eq:pdt}$$ where we suppress $T$ for conciseness and the 0 indicates that there is no solute in the system (or as is the case here, $U_{SP_0} = U_{N|SP_0} = 0$). In the thermodynamic limit of large $V >> v$ and $N >> n$, from the standard potential distribution relation [@lrp:book; @lrp:cpms; @widom:jpc82], each of the above factor on the right is simply $e^{-\beta \mu_p^{\rm ex}} / \Lambda_p^3 \rho_p$, where $\Lambda_p$ is the thermal de Broglie wavelength of the solvent sphere and $\mu^{\rm ex}_p$ is its excess chemical potential, and $\rho_p$ the density of solvent.

Since, $$Q(0,N-n,V-v,T) = \frac{1}{\Lambda_p^{3(N-n)} (N-n)!} \int\limits_{V-v}d{\vec r}_1\ldots \int\limits_{V-v} d{\vec r}_{N-n} e^{-\beta U_{N-n}} \,$$ we multiply and divide Eq. \[eq:pnporatio\] by the factor $\int\limits_{V-v}d{\vec r}_1\ldots \int\limits_{V-v} d{\vec r}_{N-n} e^{-\beta U_{N-n}}$. Rearranging the resulting equation using Eq. \[eq:pdt\] in the large $V$ and large $N$ limit, and noting that the momentum partition functions (for both solute and solvent) cancel exactly, we obtain Eq. 20 (main text), where $$\begin{aligned}
e^{-\beta \phi(R^n; \beta)} & = &  \frac{ \int\limits_{V-v} d{\vec r}_{n+1} \ldots \int\limits_{V-v} d{\vec r}_{N-n} e^{-\beta U_{N-n|SP_n}} e^{-\beta U_{N-n}}}{\int\limits_{V-v}d{\vec r}_1\ldots \int\limits_{V-v} d{\vec r}_{N-n} e^{-\beta U_{N-n}}} \nonumber \\
& = & \langle e^{-\beta U_{N-n|SP_n}} \rangle_{N-n} \,\end{aligned}$$ where $ \langle \ldots \rangle_{N-n}$ denotes averaging over the configurations of the $N-n$ solvent particles in the volume outside the observation shell.

### MaxEnt model for $\{p_n\}$ {#sc:MaxEnt}

Here we present an alternative derivation of the two parameter model (Eq. \[eq:pn\_2par\]) on the basis of information theoretic modeling of $\{p_n\}$ [@sivia; @lrp:jpcb98; @asthagiri:pre03]. On the basis of the isolated cluster partition function, we have the distribution of occupancy probabilities $\{p_n^{(0)}\}$ as $$p_n^{(0)} = \frac{K_n^{(0)} \rho_p^{n}}{1+\sum\limits_{m\ge 1} K_m^{(0)} \rho_p^{m}}$$

With $\{p_n^{(0)}\}$ as the default model and accepting the availability of the mean (first moment) and variance (second moment) of the distribution $\{p_n\}$ from simulation data, by standard maximum entropy arguments, we have $$\label{it04}
    \frac{p_n}{p^0_n} = e^{-C} e^{-\lambda_1 \cdot n} e^{-\lambda_2 \cdot n^2} \,$$ where the Lagrange multipliers $C$, $\lambda_1$, and $\lambda_2$ are, respectively, obtained from enforcing the following constraints $$\label{it05}
    \sum_n p_n = 1$$ $$\label{it06}
  \sum_n p_n \cdot n = n_{avg}$$ $$\label{it07}
    \sum_n p_n \cdot n^2= -\overline{\sigma^2}$$

We thus obtain $$\label{it11}
    p_n = \frac{\big [ e^{-\lambda_1} e^{-\lambda_2 \cdot n} \big ]^{n} K_n^{(0)} \rho_p^{n}}{1+ \sum \limits_{m\ge 1} \big [ e^{-\lambda_1} e^{-\lambda_2 \cdot m} \big ]^{m} K_m^{(0)} \rho_p^{m}}$$ Eq. \[it11\] is the same form as obtained in section \[sc:HStheory\] for a two parameter correlation. By using $\lambda' = e^{-\lambda_1}$, we can also represent Eq. \[it11\] by

$$\label{it13}
    p_n = \frac{\big [ \lambda' e^{-\lambda_2 \cdot n} \big ]^{n} K_n^{(0)} \rho_p^{n}}{1+ \sum\limits_{m\ge 1} \big [\lambda' e^{-\lambda_2 \cdot m} \big ]^{m} K_m^{(0)} \rho_p^{m}}$$

which is the same form introduced in our previous work [@bansal_structure_2016]. Based on this derivation, we can see that $ \lambda_2 $ is the term corresponding to surface interactions discussed in section  \[sc:HStheory\] and constraints the variance of the $\{p_n\}$ distribution.

Isolated cluster probabilities
------------------------------

For asymmetric mixtures we study the isolated cluster probabilities and to find the maximum number of solvent molecules that can occupy the observation volume, we use spherical code [@sphcode; @sph_code2].

### Spherical Code

The spherical code provides information to place $n$ points optimally across a sphere. It gives the optimal angle between the points, considering the center of the sphere as the origin. This is extended to place solvent particles across the solute surface, point of contact between the sphere serving as the optimizing points. The angle between the two contact points, $\theta$, is determined as given in the fig. \[fig-sph\_pack\]. The number of solvent molecule that can be tightly packed is obtained from the data [@sphcode]. The sphere onto which the points are optimally placed is an imaginary sphere which includes the critical radius as shown in the dashed lines. It is to be noted that this is still a theoretical estimate for contact packing on the imaginary larger sphere, and due to higher freedom of packing in our case, coordination states can be marginally higher for very large size ratios, $\sigma_s/\sigma_p \geq 5$. Table in fig. \[fig-sph\_pack\] also gives the ma ximum an gle for which single bonding condition holds for a given size ratio ($\theta_{c,max}$) and specified critical distance ($r_c$).

![image](figure11){width="55.00000%"}

\[fig-sph\_pack\]

Once the $n^{max}$ is defined, we find $P^{( n )}$ as the probability that there is no hard sphere overlap for randomly generated solvent molecules in the observation volume (or inner-shell) of solute molecules. As was discussed previously[@marshall_molecular_2013] that a hit-or-miss Monte Carlo [@hammersley; @pratt_quasichemical_2001] approach to calculate $P^{(n)}$ proves inaccurate for large values of $n$ ($n>8$). But since $${P^{( n)}} = P_{insert}^{( n )}{P^{( {n - 1} )}} \, , 
 \label{eq:12}$$ where $P_{insert}^{(n)}$ is the probability of inserting a [*[single]{}*]{} particle given $n-1$ particles are already in the bonding volume, an iterative procedure was be used to build the higher-order partition function from lower order one [@marshall_molecular_2013]. The one-particle insertion probability $P_{insert}^{(n)}$ is easily evaluated using hit-or-miss Monte Carlo. Following is the table for isolated cluster probabilities for different size ratios.

\[table:IC\]

[37]{}ifxundefined \[1\][ ifx[\#1]{} ]{}ifnum \[1\][ \#1firstoftwo secondoftwo ]{}ifx \[1\][ \#1firstoftwo secondoftwo ]{}““\#1””@noop \[0\][secondoftwo]{}sanitize@url \[0\][‘\
12‘\$12 ‘&12‘\#12‘12‘\_12‘%12]{}@startlink\[1\]@endlink\[0\]@bib@innerbibempty @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} [****, ()](\doibase 10.1063/1.1746468) @noop [****,  ()]{} [****,  ()](http://scitation.aip.org/content/aip/journal/jcp/78/6/10.1063/1.445245) [****,  ()](http://scitation.aip.org/content/aip/journal/jcp/82/2/10.1063/1.448475) @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [**]{} (, , ) in @noop [**]{}, , Vol. ,  (, , ) Chap. , pp.  @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [**]{} (, ) @noop [****,  ()]{} @noop [****,  ()]{} @noop [“,” ]{} () @noop [“,” ]{} () @noop [**]{} (, , )

[^1]: wgchap@rice.edu
