Coarse grained simulation approaches provide powerful tools for the prediction of the equilibrium properties of polymeric systems. Recent efforts have sought to develop coarse-graining strategies capable of predicting the non-equilibrium behavior of entangled polymeric materials. Slip-link and slip-spring models, in particular, have been shown to be capable of reproducing several key aspects of the linear response and rheology of polymer melts. In this work, we extend a previously proposed multi-chain slip-spring model in a way that correctly incorporates the effects of the fluctuating environment in which polymer segments are immersed. The model is used to obtain the equation of state associated with the slip-springs, and the results are compared to those of related numerical approaches and an approximate analytical expression. The model is also used to examine a polymer melt confined into a thin film, where an inhomogeneous distribution of polymer segments is observed, and the corresponding inhomogeneities associated with density fluctuations are reflected on the spatial slip-spring distribution.

## I. INTRODUCTION

A central challenge in development of coarse-grained modeling approaches for soft matter is that of capturing the wide spectrum of spatial and temporal scales associated with the underlying collective behavior.^{1–4} In general, coarse-graining reduces the number of degrees of freedom and finer details are only kept in an average sense. Effective interactions between surviving interacting sites become softer, and the interaction potential between two interacting sites, separated by a distance *r*, remains finite in the limit $r\u21920$. Therefore, by ignoring the finer details regarding the local non-crossability between molecules, one can achieve faster equilibration of materials over mesoscopic length scales. There are, however, systems for which microscopic details related to the non-crossability between molecules have consequences that manifest themselves at long time scales. Examples include high molecular weight polymer melts and concentrated solutions. For such materials, topological constraints play a central role on their mechanical and rheological properties.^{3,5} Several phenomenological top-down approaches have therefore been proposed to reincorporate the effect of topological interactions on coarse grained polymer models. In this work, we focus on a particular class of models that offers the potential to describe the dynamics of entangled polymer melts subject to structural or compositional inhomogeneities. These models are generally referred to as slip-spring models.^{6–11}

In recent work, we have presented a simulation approach that is able to describe the linear and non-linear rheology of entangled polymer melts.^{6} This approach relies on a multi-chain representation of the macromolecular fluid, and effective fluctuating interactions, mediated by the so-called slip-springs (SS),^{7,10} between neighboring pairs of macromolecules are introduced to represent the topological effects that arise from the non-crossability constraint. In general, there is not a correspondence between a physical entanglement and a slip-spring, and the mapping between these two concepts should be elucidated by more detailed comparisons between atomistic and coarse-grained simulations. Depending on the implementation, the total number of slip-springs can be kept constant, or it can be allowed to fluctuate through a chemical potential that determines the average molecular weight between slip-springs. In both cases, our results have been found to be in quantitative agreement with experimental data and suggest that the proposed formalism may also be used to describe the dynamics of inhomogeneous systems, such as composites, block copolymers, or polymers near surfaces. In previous work, we also made a connection between our proposed multi-chain model and the existing single-chain mean-field formalisms^{6,12} and suggest that it may be possible to develop systematic coarse-graining formalisms that start with atomistic, multi-chain representations of a material and end up with predictive constitutive models for thermodynamics and rheology.

As encouraging as our past results and those of others have been, in efforts to compare the results of our multi-chain model with those of a single-chain theory, and to do so at a modest computational expense, several assumptions were made regarding the process of creation and destruction of slip-springs between a given pair of molecules. In particular, the local environment around a specific polymer segment was considered to be constant (i.e., non-fluctuating). That assumption is inherent to single-chain descriptions, but it need not be made in particle-based simulation approaches. Importantly, for materials that exhibit inhomogeneities (e.g., copolymers, composites, or polymers near a surface), it is necessary to include the effects of a fluctuating local environment in order to capture the equation of state of the material and its coupling to dynamics. With that goal in mind, in this work, we revisit our previous model for entangled polymers and incorporate the effects of the fluctuating environment in which a polymer segment is immersed. We use this new version to obtain the equation of state associated with the slip-springs and compare it to results from an analytical expression proposed by Chappa and co-workers.^{7} We then apply the model to explore the spatial behavior of slip-springs in inhomogeneous systems, where confinement and phase segregation may play an important role on the dynamics.

## II. MODEL AND SIMULATION APPROACH

General details pertaining to the “Theoretically-Informed Entanglement Polymer Simulations” or TIEPOS approach have been presented in the literature. Here we provide a short summary of the model and methods, and highlight important differences between past implementations and the formalism employed here. Without loss of generality, we consider an AB copolymer melt consisting of *N*_{c} macromolecules in volume *V* at temperature *T*. The polymers are represented by the discretized Gaussian chain model, and each macromolecule is composed of $N\alpha $ beads, $\alpha \u2208[1,Nc]$. The system’s configuration is specified by the coordinates of all polymer beads {**r**_{i}(*s*)} at different times, where **r**_{i}(*s*) is the position of the *s*th bead in the *i*th chain. Intra-molecular interactions are given by $HbkBT=32\u2211i=1Nc\u2211s=1Ni\u22121[\mathbf{r}i(s+1)\u2212\mathbf{r}i(s)]2/b2$, where *b*^{2} is the mean-squared bond length of an ideal chain and *k*_{B} is Boltzmann’s constant. Generalizations to more complex situations, including blends of polymers with different bond lengths or chain stiffness, are straightforward. The inter-molecular interactions are expressed as a functional of the local densities $\varphi A,B(\mathbf{r})$,^{13–17}

where $\rho $ is the bead number density. The repulsion between unlike monomers, quantified by the Flory-Huggins parameter $\chi $, is represented by the first term in the functional. A finite melt compressibility is prescribed by the second term, which corresponds to Helfand’s quadratic approximation,^{18} where the parameter $\kappa $ is related to the inverse isothermal compressibility. By regularizing the scalar density fields, $\varphi \alpha (\mathbf{r})$, employing a density cloud $w$(**r**) associated with each interacting segment of the polymer chains, it is possible to transform the field-based description in Eq. (1) into a particle-based representation with a pairwise interaction potential, $HnbkBT=\u2211i<jU(\mathbf{r}ij)$, where the sum runs over all pairs of polymer segments in the system. The interaction range, *r*_{int}, should be large enough such that individual beads interact with a sufficiently large number of particles, in the simulations reported in this work, we have chosen $rint\u22480.12Re$. These intra- and inter-molecular interactions are soft and fail to capture entangled dynamics. One must therefore resort to slip-springs (or other forms of effectively describing non-crossability constraints) to represent such topological interactions.

As mentioned above, entanglement effects are represented in our work by fluctuating interactions, mediated by slip-springs (SS). A SS consist of two anchoring units that move along different chain contours, i.e., self-entanglements are excluded in the present approach, and connect the respective polymer segments via a harmonic spring: $usskBT=32[\mathbf{R}k\u2212\mathbf{R}nk]2/bss2$, where **R**_{k} denotes the position of the *k*th anchoring unit and *n*_{k} denotes its partner. The strength of the spring is quantified by *b*_{ss}, which is chosen here as *b*_{ss} = *b*. In the implementation adopted here, the SS can only hop between neighboring polymer beads, and at most one SS occupies a given segment. Here we note that previous implementations have also considered “dynamic” springs that evolve according to their own equation of motion;^{10,11,19} however, the short-time dynamics on the length scale below a coarse-grained segment has no particular physical significance.^{20} In the grand-canonical implementation considered here, the total number of SS fluctuates around an average controlled by a chemical potential $\mu $. The average number of slip-springs per chain is $\u27e8Z\u27e9$. To initialize the SS configuration, for each SS pair, a polymer bead is randomly chosen and the first anchoring unit is attached to that segment. Then, a nearby segment is chosen to place the partner anchoring unit, in such a way as to generate a Boltzmann distribution of SS lengths.

The time evolution of the system is performed by relying on a hybrid scheme, where polymer configurations are updated via dynamic simulations over a short time interval, $\tau ss$, and SS configurations are updated via Monte Carlo moves at the end of such a time interval. The slip-springs can be updated in a variety of ways; one could use a stochastic variable, $\tau ^ss$, with an arbitrary probability distribution $\Psi (\tau ^ss)$, and use this function to generate a set of $\tau ^ss$’s and update slip-springs after these different waiting times. In the interest of simplicity, here we have opted for a scheme with $\Psi (\tau ^ss)=\delta (\tau ^ss\u2212\tau ss)$. This choice provides good agreement with experiments. In this work, also for simplicity, we adopt a Brownian dynamics approach to evolve polymer segments. Thus, the time evolution of configuration {**r**_{j}} is governed by

where $\mathbf{f}jb$ and $\mathbf{f}jnb$ are the intra- and intermolecular forces acting on the polymer segment *j*, and $\mathbf{f}jss$ is the force associated with the slip-springs. Here $\bm{\xi}j$ represents the stochastic influence of the coarse-grained microscopic degrees of freedom and satisfies the fluctuation-dissipation relation: $\u27e8\xi i\alpha (t1)\xi j\beta (t2)\u27e9=2\gamma kBT\delta ij\delta \alpha \beta \delta (t1\u2212t2)$, where $\alpha ,\beta =x,y,z$, and $\gamma $ is the friction coefficient.

During the SS configuration update, one performs $\u27e8Nss\u27e9+2Nc$ attempts to update SS where *N** _{ss}*(

*t*) are the number of SS present at that time

*t*, as well as to create new SS at chain ends (see the supplementary material for details). For each attempt, an existing SS is selected with probability

*P*

_{SS}, or a chain-end is chosen randomly with probability

*P*

_{ch}= 1 −

*P*

_{SS}. In our implementation, we have chosen $PSS=\u27e8Z\u27e92+\u27e8Z\u27e9$. Then, one of the following MC moves is attempted:

A

*hopping*move, which involves the transfer of one end of a slip-spring from one polymer segment to a neighboring one, along the same polymer. We propose to move in one or the other direction along the chain with equal probability. This move is accepted with probability $p=min(1,e\u2212\Delta uss/kBT)$, where $\Delta uss$ is the difference on the slip-spring energy between the trial and old configuration.A

*destruction*move involves the deletion of an existing slip-spring pair if one of the rings reaches the chain end and moves out of the contour with probability 1/2. This move is accepted with probability $p$ = $min(1,PchPSS2NssNc\beta eW\u22121)$. Here, $W$ = $e\u2212uss(r1)/kBT+\u2211k=2Me\u2212uss(rk)/kBT$, where*r*_{1}is the bond length associated with the slip-spring to be destroyed.A

*creation*move creates a new slip-spring pair by adding a new slip-spring end to a randomly selected unoccupied chain end. The location of the ring’s partner is selected from all*M*unoccupied polymer beads within a radius*R*_{c}with probability $e\u2212uss(rj)/kBT/W$, where $W=\u2211k=1Me\u2212uss(rk)/kBT$. This move is accepted with probability $p$ = $min(1,PSSPchNc2(Nss+1)\beta e\u22121W)$.

where $\beta e=e\u2212\mu /kBT$. Thus, the value of the chemical potential dictates the average number of slip-springs per chain. Note that, as in previous implementations, the main relaxation mechanisms of entangled polymers are incorporated naturally in this approach: reptation, contour length fluctuations, and constraint release.

The important point to highlight here is the following: in this multi-chain model, at each instant of time there is a fluctuating discrete number of possible neighbors to anchor the other end of a slip-spring. Also the microstates associated with those neighboring beads can change from one time to another. All these microscopic details are encoded in $W$, which is, therefore, a fluctuating quantity. Note that in the absence of fluctuations in the local environment, as is the case of single-chain models, then $W$ is constant, let us say $Wo$, and could be adsorbed into $\beta e$, thus renormalizing the SS chemical potential, i.e., $\mu sc=\mu +kBTlnWo$. That assumption was made in a previous report.^{6} As discussed below, while its consequences are minor for homogeneous systems, e.g., homopolymer melts, it could lead to incorrect predictions for inhomogeneous materials, including block copolymers, thin films, or solutions. The implementation of the model outlined above correctly incorporates such microscopic fluctuations, and thus, $W$ is a fluctuating quantity, varying both in space and time.

## III. RESULTS

### A. Homogeneous systems

In this section, we explore the effects of incorporating a fluctuating $W$ into a multi-chain description of homopolymers, as indicated in Section II. Note that $W$ depends on particle density $\rho $; the higher the number density, the smaller become the relative fluctuations of $W$. If $W$ is assumed to be a constant, it plays no role on the acceptance probabilities and, therefore, the average number of slip-spring per chain $\u27e8Z\u27e9$ is no longer sensitive to changes in particle density. In the acceptance rules outlined above, $W$ appears explicitly, and number density therefore affects the resulting $\u27e8Z\u27e9$.

We fix the value of the chemical potential in both schemes and vary the polymer bead density (the results reported in this section use *N* = 32, and the length is measured in units of *R*_{e}). It is important to remark that the chemical potential in both schemes does not have the same numerical value and does not depend on $\rho $ either. Figure 1 shows $\u27e8Z\u27e9$ as a function of density $\rho $. As expected, for constant $W$, the number of slip-springs per chain does not depend on the bead density for large enough values (green dashed curve). In contrast, for a fluctuating $W$, the new set of acceptance rules gives rise to a monotonic dependence on $\rho $. In fact, for relatively small densities, a linear dependence is found, which agrees with the prediction of Ref. 7 for a model similar to that considered here: This behavior can be explained by invoking an analytical prediction for the SS equation of state, originally derived in Ref. 7, which reads^{34}:

where $g(\mathbf{r})$ is the slip-spring pair correlation function, which only includes pairs of segments that can be connected by a slip-spring and $\beta e=1/\zeta =exp(\u2212\mu /kBT)$ denotes the inverse fugacity of SS. Note, in particular, that if the melt’s microscopic structure did not depend on density, i.e., if *g*(*r*) were not affected by changes of $\rho $, then for a fixed chemical potential $\u27e8Z\u27e9$ would linearly grow with $\rho $. However, if by changing $\rho $ there is a change of the slip-spring pair correlation function, then a non-linear dependence should be expected.

One key difference between the model in Ref. 7 and ours is that in the model of Chappa *et al.*^{7} a single segment can be occupied by any number of SS, only SS-connections to the identical bead or nearest neighbors along the same chain molecule were excluded. In our case, only one SS can be attached to a given polymer segment and the two ends of a SS have to reside on different polymers. We stress, however, here that at this point it is still unclear whether one or multiple SS per bead provide a better description of entangled polymers.

Interestingly, the curve $\u27e8Z\u27e9$ vs $\rho $ departs from a linear relationship for large densities. We computed the slip-spring pair correlation function for different particle densities. In Figure 2, we show *g*(*r*) for two extreme cases: $\rho $ = $400/Re3$ and $\rho $ = $8000/Re3$. As can be seen in the figure, the local microscopic structure has a strong dependence on the particle density, thereby affecting the number of SS per chain. Since an increase in density gives rise to a larger *g*(*r*) at short distances, the number of potential SS-partners of a segment increases faster than $\rho $, and Eq. (3) predicts an enhanced increase of $\u27e8Z\u27e9$ with density, contrary to what is observed in Fig. 1. Thus, the deviations from the linear behavior in Figure 1 indicate a saturation in the number of SS due to the additional constraints of single-segment occupancy of the intermolecular SS. Also, we note again that the fact that we are excluding more than a single SS per segment induces effective contact interactions (and associated correlations) between SS that were not considered in the analysis leading to Eq. (3).

If the number density is fixed at a given value ($\rho $ = $1798/Re3$), we can explore how $\u27e8Z\u27e9$ is affected by the chemical potential. As reported in Figure 3, a monotonic dependence of the chemical potential is obtained. Again the growth of $\u27e8Z\u27e9$ with SS fugacity $\zeta \u2261exp(\mu /kBT)$ agrees with that expected from Eq. (3). However, in this case, a linear relationship was reported in Ref. 7. The origin of this discrepancy between both approaches relies on one key difference: first, as mentioned above, in the approach considered here each polymer bead can be occupied by only one slip-spring and the two ends of a SS are connected to different polymers, whereas in Ref. 7 there are no such restrictions on the SS connectivity. At a fixed density, the incorporation of slip-springs has no effect on the microscopic structure of the melt, even for large values of the SS chemical potential (see Figure 4), i.e., for this density the compensating potential that eliminates the effect of SS on equilibrium properties^{7} is not necessary. Therefore, as the chemical potential is increased, the number of available unoccupied beads decreases, giving place to a slower growth of $\u27e8Z\u27e9$ with $\mu $. In the model introduced in Ref. 7, any number of SS on a single bead is permitted, and increasing $\zeta \u2261exp(\mu /kBT)$ leads to a linear increase of the number of SS that can be accommodated between any neighboring pairs of segments except nearest neighbors along a chain.

We conclude this section by comparing the analytical equation of state obtained for the multi-chain system (Eq. (3)) with that predicted by the single-chain theory of Schieber *et al.*^{21–24} In the latter theory, for discrete chains, the equation of state reads

where $\beta esc=e\u2212\mu sc/kBT$, the superscript “sc” refers to the single-chain prediction, and, $\beta esc$ is chosen to be equal to the average number of beads between entanglements, *N*_{e}. Note that the theory of Schieber *et al.* excludes double SS occupancy by construction, whereas Eq. (3) captures the multi-chain correlations. Our model contains both ingredients, and the following considerations provide a reasonable approximation for dense systems that allows a comparison. For weakly entangled melts, the single-chain equation of state can be approximated by $\u27e8Z\u27e9\u2248N\beta esc$. On the other hand, ignoring that in Ref. 7 SS must not connect the same or nearest-neighbor beads along a chain, we can approximate Eq. (3) by

where

is the average of the Rosenbluth factor defined in the creation-destruction rules described above. In the case of homogeneous systems, $\u27e8W(\rho )\u27e9$ is independent of space and time. Thus, for weakly entangled homopolymers in bulk, we can identify $\beta esc$ = $\beta e\u27e8W(\rho )\u27e9\u22121$, i.e., the chemical potentials are related by $\mu sc$ = $\mu +kBTln\u27e8W(\rho )\u27e9$. To corroborate this approximate relationship, we first plot the data in Figure 3 in terms of $\beta e$ = $exp(\u2212\mu /kBT)$. This is shown in the main panel of Figure 5. As can be seen, a similar dependence of $\u27e8Z\u27e9$ on $\beta e$ is obtained (see Figure 2 in Ref. 1), but the latter has quite different numerical values. Then, we use the time average of $W$ sampled during the simulations, $\u27e8W\u27e9$ and plot $\u27e8Z\u27e9$ vs $\mu +kBTln\u27e8W\u27e9$. The resulting curve is presented in the inset of Figure 5 along with the equation of state predicted by the single-chain theory, Eq. (4). As can be seen, the simulations follow the analytical predictions, even for systems with small $\beta esc$. Therefore, for homogeneous systems, the results presented in Ref. 1 provide a reasonable approximation to the results obtained using the fluctuating-environment creation-destruction rules adopted here. Note, however, that for inhomogeneous systems the spatiotemporal character of $W(r,t;\mu )$ will affect the creation-destruction rules locally, possibly leading to inhomogeneous slip-spring distributions, as expected, for example, in phase segregated block copolymers. A specific example is considered in Section III B.

### B. Inhomogeneous systems: Block copolymers and thin films

In what follows, we consider two inhomogeneous a homopolymer melt confined into a thin film, and a phase segregated block copolymer thin film. In both cases, instead of evolving the polymer configurations by dynamical equations (Eq. (2)), we have chosen to perform local Monte Carlo moves to relax polymer chains.^{28} The reason for this has to do with the fact that hard boundary conditions are easier to handle with an MC approach, and we are only interested in the polymer and slip-spring density distributions, as opposed to dynamic properties. As above, polymer segments are relaxed by a certain number of MC steps; the SS configuration is updated using the rules described in Section II. The coarse grained model we use requires that we use a large chain discretization, *N*, in order to describe the narrow interface of a nearly incompressible polymer melt in the vicinity of a solid.^{25} That is computationally demanding, and we have therefore opted here for an effective bonded interaction which helps to reproduce the configurational statistical properties of Gaussian chains under confinement with a low value of *N*.^{25} Note that an alternative technique, particularly useful when the statistical segment length is smaller than the polymer-solid interface width, would be to use an external potential obtained via Self-Consistent Field Theory (SCFT)^{26} or use a higher-order discretization of the Gaussian paths.^{27} Explicitly, we add the extra term to the bonded contribution^{25}

where *z*_{1} and *z*_{2} are the perpendicular coordinates, with respect to the hard surface, of two connected segments along the backbone. This equation is obtained by modeling the hard walls as infinite step functions and solving the diffusion equation with Dirichlet boundary conditions on the impenetrable surface (see Ref. 25 for details). This term compensates for the lost of accessible configurations close to the hard interface. The films have a thickness *L*_{z} and the hard walls are located at *z* = 0 and *z* = *L*_{z}. In the directions parallel to the surfaces, we impose periodic boundary conditions. Using these effective interactions, it is possible to recover the density profiles predicted by the field theory with *N* as low as *N* = 32. In Figure 6, we present the density profile for polymer segments (black curve) of a melt confined into a thin film of thickness *L*_{z} = 1.25*R*_{e}. As can be seen in the figure, there is a depletion layer close to the surfaces, where most of the segments are chain-ends (red dashed curve), which are entropically favored to be localized near the walls. Thus, this breaking of translational symmetry induces changes to the polymer configurations and it becomes stronger as the thickness of the films gets smaller. As a consequence of the inhomogeneous distribution of polymer segments and chain configurations, the slip-springs also distribute inhomogeneously inside the film (green dashed curve); as can be appreciated in the figure, most of the entanglements are located in the middle of the film. We should note here that recent work^{29–31} addressing the effect of strong confinement on polymer configurations and entanglements in microscopic models has reported an entanglement spatial distribution that displays the same qualitative behavior as the one reported in here. Importantly, not only does the distribution of slip-springs become inhomogeneous, but their number decreases as the films become thinner (see Figure 7). This effect is consistent with experimental observations showing that entangled melts have a lower viscosity as they undergo nanoscale confinement,^{32,33} which is an indirect indication that chains become disentangled under such circumstances. As mentioned above, in our current approach, we have not considered self-entanglements, which are likely to arise in these thin film systems.

An analogous situation arises in lamellae-forming block copolymers. For these materials, we chose a discretization parameter *N* = 64, and the Flory-Huggins parameter is chosen as $\chi N=37$, with a lamellar periodicity of *L*_{o} = 1.79*R*_{e}. The melt is confined into a thin film (thickness *L*_{z} = *L*_{o}, and we have also made use of Eq. (7)) and the bottom substrate is patterned with stripes of width $w$ = *L*_{o}/2, and periodicity equal to *L*_{o}. The stripes are preferential for block A, and the space between the stripes is preferential for block B. Thus, the translational symmetry is broken in two dimensions: perpendicular to the substrate and perpendicular to the A-B interfaces. As expected, there is an inhomogeneous distribution of slip-springs in both directions. More importantly, for the parameters we have chosen, there is a minimum in the slip-spring spatial distribution located at the A-B interface, and a shallow minimum in the middle of the domains (Figure 8). Entanglements are in this case localized in the middle of the blocks, as opposed to the chain ends. The latter results are quite promising, as they suggest some correlations between molecular conformations and slip-springs, as expected to occur for real entanglements. We should note here that it could be possible that the local chain orientation would play a stronger role defining the local entanglement length; thus SS update rules could require a term associated with the local orientation. To address this issue, more comparative studies between microscopic models and the current approach are needed.

## IV. CONCLUSIONS

A previously reported multi-chain slip-spring formalism has been extended to incorporate the effects of a fluctuating local number of entanglements. The simulation algorithm was used to obtain the equation of state associated with the slip-springs, and its results were found to be in good agreement with predictions from an analytical expression.^{7,21} The formalism was also used to examine the effects of confinement and phase separation on the entanglement distribution. In both cases, these inhomogeneities are clearly reflected on the spatial distribution of slip-springs and are expected to play a major role in the dynamic response of the material. Results for the linear behavior and non-linear rheology of these materials will be reported in a future publication.

## SUPPLEMENTARY MATERIAL

See supplementary material for a detailed description of the algorithm to update SS configurations.

## ACKNOWLEDGMENTS

This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Science and Engineering Division. L.S. and M.M. thank the German Science Foundation for financial support (No. DFG Mu 1674/16-1). Marat Andreev gratefully acknowledges the support from NIST through a CHiMaD postdoctoral award. We gratefully acknowledge the computing resources provided on Blues, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory.

## REFERENCES

Assuming that there is not a constraint on how many SS a single bead can hold, the physical interpretation of Eq. (3) is as follows: the probability of a given pair of beads, *i* and *j*, to have a SS is $e\mu /kBT\u2212uss(r)/kBT$, where *r* is their separation. Now, the probability of having such separation is given by *g*(*r*)/*V*, thus the average number of SS between this pair is: $\u27e8Zij\u27e9=e\mu /kBTV\u222bd3rg(r)e\u2212uss(r)/kBT$. The total number of SS is obtained by summing over all possible pairs, and from that Eq. (3) follows (see Ref. 7 for details).