Depletion force between disordered linear macromolecules

TL;DR

This paper investigates the depletion force between disordered linear macromolecules by combining analytical calculations, molecular dynamics simulations, and theoretical approaches to understand how disorder affects entropic interactions.

Contribution

It introduces a method to calculate the depletion force for disordered macromolecules modeled as random hard sphere arrangements, combining analytical and simulation techniques.

Findings

Depletion force depends on the position along the macromolecule.

Disorder influences the magnitude and profile of the depletion force.

Comparison shows differences between stiff and flexible disordered polymers.

Abstract

When two macromolecules come very near in a fluid, the surrounding molecules, having finite volume, are less likely to get in between. This leads to a pressure difference manifesting as an entropic attraction, called depletion force. Here we calculate the density profile of liquid molecules surrounding a disordered rigid macromolecules modelled as a random arrangement of hard spheres on a linear backbone. We analytically determine the position dependence of the depletion force between two such disordered molecules by calculating the free energy of the system. We then use molecular dynamics simulations to obtain the depletion force between stiff disordered polymers as well as flexible ones and compare the two against each other. We also show how the disorder averaging can be handled starting from the inhomogenous RISM equations.

Equations

$$\gamma (r,\kappa )=\{\begin{array}{ll}\kappa \sqrt{1-{\left(r/\kappa \right)}^2}& r\le \kappa \\ 0& r>\kappa \end{array}.$$

$${\Gamma }_N$$

$${\tilde{\gamma }}_{\theta }(r,\kappa )=\{\begin{array}{ll}0& r\ge \kappa \\ \gamma (r,\kappa )& \kappa >r\ge r_0\\ (r_0-r)\tan \theta +\gamma (r_0,\kappa )& r_0>r\ge 0\end{array}$$

$$\displaystyle Pr(r|x_{p}\in[x_{k},x_{k+1}])=\frac{\Xi\big{(}x_{k}-x_{k+1}-2R_{1}\xi(r)\big{)}}{x_{k+1}-x_{k}}$$

$$Pr(r|\vec{x})$$

$$\displaystyle=\frac{1}{L}\sum_{k=1}^{N}\Xi\big{(}x_{k}-x_{k-1}-2R_{1}\xi(r)\big{)}$$

$$p_{\varphi }(r,N)$$

$$\frac{1}{L}\sum_{k=1}^{N}\Xi\big{(}x_{k}-x_{k-1}-2R_{1}\xi(r)\big{)}=1-\phi\,\xi(r)$$

$$p_{\varphi }(r,N,\theta )=1-\varphi \xi (r)\text{mbox}for\xi (r)\le 1.$$

$$p_{\varphi }($$

$$\{\begin{array}{ll}1-\varphi \xi (r)& \xi (r)\le 1\\ \left(1-\varphi \right){\left(1-\varphi \frac{\xi (r)-1}{N(1-\varphi )}\right)}^N& 1<\xi (r)\le (\frac{N}{\varphi }-N)+1\\ 0& N(\frac{1}{\varphi }-1)+1<\xi (r)\end{array}$$

$$p_{\varphi }(r,\theta ,\kappa )=\{\begin{array}{ll}1-\varphi \xi (r)& \xi (r)\le 1\\ (1-\varphi )\exp \left[\frac{\varphi \left(1-\xi (r)\right)}{1-\varphi }\right]& \xi (r)>1\end{array}$$

$$P(x)=(1-\varphi )\exp [\varphi (1-2x)/(1-\varphi )].$$

$$V_E(x)=V-v_{\text{out}}-Lx+L\lambda (x).$$

$$\lambda (x)$$

$$={\int }_0^x{p}_{\varphi }(r)dr+{\int }_x^2{p}_{\varphi }(r)p_{\phi }(r^{\prime })dr+{\int }_2^{2+x}{p}_{\phi }(r^{\prime })dr$$

$$F(x)=F_0-N{k}_{B}T\log \left[1-\frac{v_{\text{out}}}{V}-\frac{L}{V}\left(x-\lambda (x)\right)\right]$$

$$F(x)\simeq F_0^{\prime }-\rho k_{B}T\left[x-\lambda (x)\right]L.$$

$$\mathcal{P}(x)$$

$$u(x)$$

$$g_{ll}(x)=\exp \left(-\rho k_{B}T\left[x-2-\lambda (x)+\lambda (2)\right]L\right).$$

$$h({\mathbf{r}}_{01})=\int d{\mathbf{r}}_{2}d{\mathbf{r}}_3\omega ({\mathbf{r}}_{02})C({\mathbf{r}}_{23})\left[\omega ({\mathbf{r}}_{31})+\rho h({\mathbf{r}}_{31})\right]$$

$${\omega }_{m,k}^{(\alpha )}(\mathbf{r})=\delta (r_y-h_k^{(m)}+p_y^{(m)})\times \delta (r_x-p_x^{(m)}).$$

$$h_{l{l}^{\prime }}^{(\alpha )}(\Delta p_x)=\rho \sum _{j,k=1}^N{\int }_{{\mathbb{R}}^2}d{\mathbf{r}}_s{C}_{ls,j}^{(\alpha )}\left((x_s-p_x^{(l)})\hat{x}+(y_s-h_j^{(l)})\hat{y}\right)$$

$$\times C_{s{l}^{\prime },k}^{(\alpha )}\left((p_x^{(l^{\prime })}-x_s)\hat{x}+(h_k^{(l^{\prime })}-y_s)\hat{y}\right).$$

$$-\beta u(x)=⟨\log \left(1+h^{(\alpha )}(x)\right)⟩=⟨h^{(\alpha )}(x)⟩+\mathcal{O}({\rho }^2)$$

$$=\rho \sum _{j,k=1}^N{\int }_{{\mathbb{R}}^2}d{\mathbf{r}}_s⟨C_{ls,j}^{(\alpha )}({\mathbf{r}}_{l,j})C_{s{l}^{\prime },k}^{(\alpha )}({\mathbf{r}}_{m,k})⟩$$

$$=\rho {\int }_0^{h+2}d{x}_s{\int }_0^{L}d{y}_s⟨\sum _{j,k}{C}_{ls,j}^{(\alpha )}({\mathbf{r}}_{l,j})C_{s{l}^{\prime },k}^{(\alpha )}({\mathbf{r}}_{m,k})⟩.$$

$$\chi (x)$$

$$=L{c}_{ls}(x-p_x^{(l)})c_{ls}(p_x^{(l^{\prime })}-x)$$

$$c_{ls}(x-p_x^{(l)})$$

$$\displaystyle C^{(\alpha)}_{ls,i}(\mathbf{r})=\begin{cases}0&\|\mathbf{r}-\mathbf{p}_{i}\|\geq 2R\\ -\Big{[}\sum_{j}\Theta\left(2R-\|\mathbf{r}-\mathbf{p}_{j}\|\right)\Big{]}^{-1}&\|\mathbf{r}-\mathbf{p}_{i}\|<2R\end{cases}.$$

$$\chi (x)$$

$$-\beta u(x)$$

$$\frac{u(x)}{\rho k_{B}TL}=-x+{\int }_0^{x+2}{p}_{\varphi }(r)p_{\phi }(x+2-r)dr$$

Full text

Objects immersed or dissolved in a liquid will experience an emergent attractive force, as the liquid molecules, having finite volume, cannot squeeze between them Likos (2001). Put another way, it is entropically favorable for the objects to be close, since each have surrounding volumes unavailable to the liquid molecules; and when objects approach to the extent that these volumes overlap, the molecules have more volume to explore Asakura and Oosawa (1954, 1958). Thus objects are more likely to be near each other, as if they attract. This entropic force is called “depletion.”

Since the seminal papers of Asakura and Oosawa, depletion forces has had far reaching implications from molecular physics Minton (1981); Gholami et al. (2006); Mravlak (2008); Maghrebi et al. (2011); Sapir and Harries (2015) and biochemistry Hall and Minton (2003); Hanlumyuang et al. (2014); Braun et al. (2016), to high energy physics Visser (2011); Basilakos and Sola (2014); Typel (2016); Feng et al. (2016); Fadafan and Tabatabaei (2016). It has even been suggested that gravity Verlinde (2011); Plastino and Rocca (2018); Plastino et al. (2018); Bhattacharya et al. (2018) and the Coulomb force Wang (2010) might be entropic forces. So far, depletion forces between plates immersed in rods Mao et al. (1997a); Chen and Schweizer (2002) or spherocylinders Mao et al. (1997b), forces between colloids Crocker et al. (1999), semiflexible chains Castelnovo and Gelbart (2004), spherocylinders Li and Ma (2005), and ellipsoids König et al. (2006); Miao et al. (2014) immersed in colloids, and forces between colloids immersed in polymer Striolo et al. (2004) have been established. Forces mediated by mixtures of two types of particles have also been studied Triplett and Fichthorn (2010).

Depletion interactions are particularly relevant in polymer physics Fuchs and Schweizer (2001, 2002); Tuinier et al. (2015); Likos (2016); Pérez-Ramírez et al. (2017). The Ornstein-Zernike (OZ) equation Ornstein and Zernike (1914); Schöll-Paschinger and Kahl (2003); Brader and Schmidt (2013), Reference Interaction Site Model (RISM) Misin et al. (2015); Tormey (2016); Johnson et al. (2016), and Polymer RISM (PRISM) Schweizer and Curro (1994); Schweizer (1993); Perry and Sing (2015) yields correlation functions between flexible polymers, rods, and colloids. Generally, additional assumptions are needed to close OZ and RISM-type equations, since they each contain multiple functions whose forms are not specified by the defining equation Baxter (1968); Yethiraj and Schweizer (1992); Martynov et al. (1999); Saumon et al. (2012). Closure equations specify properties of the functions involved in PRISM, and include the Percus-Yevick approximation Percus and Yevick (1958) and the hypernetted-chain equation Fries and Patey (1985). These models and closure relations are then typically solved numerically, though some analytical methods also exist Wertheim (1963, 1964); Yasutomi and Ginoza (2000); Adda-Bedia et al. (2008); Rohrmann and Santos (2011).

Depletion forces are also important in biophysics, where they play a crucial role in the organization of cells, particularly in chromosomal organization and compaction Stavans and Oppenheim (2006); Jeon et al. (2016a, b); Pereira et al. (2017); Shendruk et al. (2015); Kumar et al. (2019), the behavior of DNA inside cells Odijk (1998), crowding-induced chromatin compaction Oh et al. (2018), and polymer chain looping Bian et al. (2019).

While depletion forces always originate from disordered arrangements of a solvent, the objects experiencing the force themselves have always been chosen by authors to be orderly geometric shapes, such as planes, cylinders, squares, and spheres. In polymer physics, while RISM allows for the incorporation of intramolecular structure, the structure function can only represent specific objects rather than an ensemble of disordered objects.

Here, we study the depletion force between two disordered objects. Specifically, we consider rigid and flexible disordered polymers (a random arrangement of hard spheres) immersed in a solvent (also hard spheres). These objects can be seen as a model of disordered block copolymer, where one monomer is very small and the other is large Stepanow et al. (1996); Westfahl Jr and Schmalian (2005); Kim et al. (2018). We first find the probability that a sphere incident towards a straight polymer will approach by $r$ before any collision. This is the disorder averaged polymer-liquid correlation function in the low density limit. We then use this to evaluate the average depletion force between two straight polymer chains. Molecular dynamics simulations are used to verify our theoretical results. We then numerically calculate the depletion force for the more general case of flexible polymer chains. Finally, we show how RISM can be extended to treat the straight disordered polymers.

Liquid profile near a disordered chain. We model a rigid disordered polymer as a line of length $L$ with $N$ non-overlapping, but otherwise randomly-placed spheres of radius $R_{1}$, and the surrounding liquid molecules as spheres of radius $R_{2}$.

Consider a liquid molecule colliding with one of the monomers as shown in Fig. 1A, with a dotted empty circle on the right, and a filled gray circle respectively. We define our coordinate system such that the polymer lies on the $x$-axes horizontally, and that the monomer and liquid molecule collide at a point where the vertical coordinate of the liquid molecule is $y=r$. At the point of contact, the difference in the $x$ coordinates between the monomer and the liquid molecule is

[TABLE]

where $\kappa=1+R_{2}/R_{1}$. The “stopping function,” $\gamma(r,\kappa)$, is useful because it tells us that a liquid molecule $r$ away from the backbone cannot have a horizontal distance less than $\gamma(r,\kappa)$ to a monomer. This region is excluded since the two disks would have to overlap.

Here and throughout $\gamma(r)$ and $r$ are defined in units of $R_{1}$. For the rest of the paper we pick $R_{1}=R_{2}$ for cosmetic and pedagogical reasons, however the $R_{1}\neq R_{2}$ case is a straightforward generalization. Also note that much of our analysis can be further generalized to non-spherical shapes by replacing $\gamma$ with another appropriate stopping function.

We also define a second related quantity $\tilde{\gamma}_{\theta}(r,\kappa)$. We consider a liquid molecule incident at angle $\theta$ that passes by the monomer tangentially, as shown in Fig. 1A with a dotted empty circle on the left. When $r$ away from the backbone, its horizontal distance to the monomer is defined as $\tilde{\gamma}_{\theta}$. Note that $\tilde{\gamma}_{\theta}\geq\gamma$ due to the “shadow” of the disk on the line. We will evaluate $\tilde{\gamma}_{\theta}$ in a moment.

To estimate the liquid density profile around the chain, we calculate the probability that a molecule incident towards the chain at an angle $\theta$ can approach to a distance $r$ without contact. Given a configuration of spheres at coordinates $\{x_{0},...,x_{N-1}\}$, and a incidence angle $\theta$, we will count up the fraction of the interval $[0,L]$ in which the incident sphere could start to make it to a distance at least $r$ above the central axis (cf. Fig 1). We will then integrate over all valid configurations of spheres.

The measure of configurations of spheres on a line is equivalent to the configurations of hard lines on a line, known as a Tonks gas Tonks (1936); Giaquinta (2008). We write it as

[TABLE]

where $\phi\equiv 2NR_{1}/L$ is the (linear) volume fraction of the disks on the line, i.e. the amount of the line that is covered by disks.

As the incident sphere approaches the chain, it will at some point be at a distance $r$ to the backbone. The $x$ coordinate of the incoming sphere when it is $r$ above the backbone will be labeled $x_{p}$. We will approach the computation by conditioning on the interval that $x_{p}$ falls into.

Each sphere in the chain excludes a distance $R_{1}\gamma(r)$ to its right, and $R_{1}\tilde{\gamma}_{\theta}(r,\kappa)$ to its left, where

[TABLE]

where $r_{0}=r_{0}(\theta)\equiv\kappa\sin\theta$. The case for $\kappa\tan\theta>r$ emerges because of the “shadow” the excluded volume that a sphere casts on the line.

Let $x_{p}\in[x_{k},x_{k+1}]$ (note that these $x$’s are have dimensionality of length). Then probability that it is able to approach to a distance $r$ above the line is proportional to the length $[x_{k},x_{k+1}]$ that is not excluded by $\gamma$ or $\tilde{\gamma}$

[TABLE]

where $\xi(r)\equiv[\gamma(r)+\tilde{\gamma}(r)]/2$, and the clamp function $\Xi(x)$ is $\Xi(x)=x$ for $x>0$ and 0 otherwise. The clamp function is necessary, because the excluded areas of $\gamma$ and $\tilde{\gamma}$ may overlap. Since the particle’s intersection position, $x_{p}$ will be in $[x_{0},x_{1}]$ or in $[x_{1},x_{2}]$, etc., and the particle is initialized uniformly at random,

[TABLE]

Therefore, to obtain $p_{\phi}(r,N)=Pr(r)$ we just need to integrate over all the arrangements of $x_{k}$,

[TABLE]

where the integration measure is the same as (2), and $1/\Gamma_{N}$ comes from the normalization of the measure, since we are averaging over all valid configurations of spheres.

If $r>2$, then $p_{\phi}(r,N)=1$. If $r$ is large enough that $\xi(r)\leq 1$, then every $\Xi$ function is non-zero, and the sum in (6) telescopes and can be evaluated easily,

[TABLE]

which is independent of any of the positions of the spheres on the line. This term then pulls out of the integral, which cancels with the factor of $1/\Gamma_{N}$, leaving us with

[TABLE]

If $\xi(r)>1$, integrating (6) is more involved, since the clamp functions can be zero (see Appendix A),

[TABLE]

In the thermodynamic limit $N,L\to\infty$, $N/L=\phi/(2R)$,

[TABLE]

Note that $\xi(r)$ is actually a function of $\kappa$, $\theta$, and $r$, and therefore dimensionless.

To confirm this result, we carry out 2D simulations where we launch 50,000 disks from random positions towards a line of randomly arranged disks, and recording the $r$ value at which they hit. We do this at both $\theta=0$ (Fig. 2 A) and at various nonzero $\theta$’s (Fig. 2 B) and find excellent agreement between analytic formulas and simulations.

The approach probability (9) can be viewed in several additional ways. One way of looking at (9) is that its limiting case is the probability that there is a gap of size $x=\xi(r)/2$ in a Tonks gas. Indeed, the probability that there is a gap of size $x$ at a generic point in a Tonks gas in the thermodynamic limit is Torquato et al. (1990); Elkoshi et al. (1985); Giaquinta (2008)

[TABLE]

Another view, which we use to find the entropic force, is that $p_{\phi}(r)\equiv p_{\phi}(r,\theta=0)$ is the expected free volume $r$ away from the chain, see Fig. 3 A. Indeed, our simulations confirm that even for moderately flexible polymers,” $p_{\phi}(r)$ is a good approximation of the pair correlation between the polymer and the liquid (Fig. 4). Figs. 3A and 4A, B, were generated by molecular dynamics simulations of flexible or inflexible polymers immersed in a fluid. See Appendix D for more details on how the simulations were performed.

Note that the fluid density profile predicted by (9), $p_{\phi}(r)$, is different than that of an ordered chain, where the monomers are evenly spaced. For an ordered chain, the approach probability is $p_{\phi}^{(\text{u})}(r)=\Xi(1-\phi\,\gamma(r))$, which is significantly different for $r<\sqrt{3}$. Interestingly, we see that the approach probability for an ordered chain is always less than or equal to that of a disordered chain with the same average density (Fig. 3B).

Of course, the interpretation of (9) as free volume or pair correlation function only holds exactly when the liquid density $\phi_{L}=Nv/V$ is small. When the liquid density is large, packing effects and sphere-sphere exclusion can come into play, and the density profile of the liquid will deviate from $p_{\phi}(r)$ (Fig. 4B).

Depletion force between rigid disordered chains. Now that we know the free volume near a disordered polymer, we can calculate the depletion force between two rigid polymer chains. While we present the derivation for polymers in 2D, the idea is the same for 3D (see Appendix C). Consider two parallel disordered lines of length $L$ in two dimensions with sphere densities $\phi$ and $\varphi$ a distance $X=(2+x)R$ away from one another. Suppose that the lines are in a hard sphere fluid, also with radius $R$, and with number density $\rho=N_{s}/V$, total system volume (area) $V$, and temperature $T$. If the lines are closer than $x=2$, the excluded volumes of the spheres on separate lines can overlap, resulting in an entropic force. The entropic force will depend the arrangement of spheres on each line, which are random variables, but we can calculate the average entropic force between the lines (which becomes exact as the line length $L\to\infty$) using (9). From here on, by $p_{\phi}(x)$, we mean $p_{\phi}(x,0,0)$, i.e. $\theta=0,\kappa=2$.

If the lines themselves cannot interact ($x>0$), then the arrangements of spheres on each line is independent. The expression $p_{\phi}(x)$, (9), tells us the probability that a point a distance $x$ from a line can be occupied by a sphere. Letting $rR$ denote the distance of a point from the left line, $r^{\prime}R=(2+x-r)R$ is the distance of that point from the right line. For points at distances $0<r<x$, a sphere can only be excluded by spheres on the left line, the probability that a sphere can occupy the volume is $p_{\phi}(r)$. For points at distances $x<r<2$, a sphere can be excluded by spheres on either line, so the probability that a sphere can occupy the volume is $p_{\phi}(r)p_{\varphi}(x+2-r)$. Finally, for points at distances $2<r<2+x$, a sphere can only be excluded by the right line, so the probability that a sphere can occupy the volume is $p_{\varphi}(x+2-r)$.

Depletion Force. Each of these expressions involving $p$’s is the expected free volume per length at points between the lines. From this, we can calculate the free energy of the system, and then the entropic force. The partition function is $Z(x)=V_{E}(x)^{N}/(N!\Lambda^{2N})$, where $\Lambda=h/\sqrt{2\pi mk_{B}T}$ is the kinetic part (the $h$ in $\Lambda$ is Planck’s constant), and $V_{E}(x)$ is the free volume of the system. The expected free volume of the system is

[TABLE]

The term $v_{out}$ corresponds to the reduced expected volume to the left of the left chain and right of the right chain, which can also be expressed in terms of $p$, but does not depend on $x$, so it will not matter to the entropic force calculation. In the third and fourth term, we subtract the volume between the lines, and add back the expected volume between the lines. From the considerations in the last section, we know that

[TABLE]

where $r^{\prime}\equiv x+2-r$. The free energy $F=-k_{B}T\log Z$ is

[TABLE]

where $F_{0}$ is independent of $x$. In the thermodynamic limit $V\to\infty$, we use $\log(1+\epsilon)\simeq\epsilon$ to get

[TABLE]

From this, we find that the force per unit length $\mathcal{P}(x)=-L^{-1}\partial F/\partial x$, and disorder averaged potential of the force $u(x)$ are

[TABLE]

which is our second main result. Note that the expression is the same if $\phi$ and $\varphi$ are exchanged (and integrating by parts), and so $\mathcal{P}$ can be written in a symmetric form if that is preferred. Furthermore, $\mathcal{P}$ is always negative, so the entropic force is attractive, and $\mathcal{P}$ is zero for $x>2$.

The correlation function for two chains can be easily found from (12),

[TABLE]

To verify our result, (11), we ran molecular dynamics simulations 111Using the GFlow molecular dynamics package, https://github.com/nrupprecht/GFlow for $\phi=\varphi$ measuring the force on disordered lines, and binning this by line separation. The lines were arranged to be rigid, and constrained to move only in the horizontal direction in order to measure the attractive force accurately (Fig. 5).

As molecular dynamics requires forces, we simulate the hard spheres as “very stiff elastic spheres”, which turns out to be very good approximation. See Appendix D for more simulation details.

Depletion force between flexible disordered chains. We evaluate the entropic force between flexible chains using molecular dynamics simulations (Fig.5A,B). Since we construct flexible chains of macromolecules instead of rigid linear configurations of macromolecules, the horizontal distance between chains is now a stochastic quantity both in position and time.

To get around this problem, we calculate the force on a per particle, i.e. the depletion pressure, which is the time averaged pressure $F/(2R)$ on a particle as a function of distance, $x$, from the closest point of the other chain. This way, we can also meaningfully compare depletion pressure of a flexible polymer with the net force on the rigid polymer per $2R$ length (cf. Fig.5).

While the entropic pressure between flexible polymers is similar to the entropic pressure between rigid linear polymers, several interesting differences are apparent (Fig. 6). One is that the force between flexible chains are always less than that between the rigid chains. A second interesting feature is that while the force between rigid chains is zero beyond $x=2$, for flexible chains it becomes repulsive.

The repulsive part of the force is due to the fact that a flexible chain has many configurations it can be in, unlike a rigid chain. A snapshot of the instantaneous configuration of the two chains can be seen in Fig. 6 B. When two chains are near each other, they block each other from folding in as many ways as they could if they were far away from one another. This leads to an entropic repulsion term, which is what causes the pressure to be smaller overall than that between rigid polymers. Since the entropic repulsion depends on the number of allowable configurations for the entire chain, an entropic force can be felt by particles that are not close enough to feel the fluid entropic force that is due to fluid arrangements (i.e. particles with $x>2$). The force is entirely due to the chain entropic force.

Connection to RISM. Our method can be seen as an alternate route to solving the Reference Interaction Site Model (RISM) equation, a technique for finding correlation functions and depletion forces that is common in the polymer physics literature. The inhomogeneous RISM equation Ishizuka et al. (2008) reads

[TABLE]

where $h$ is the matrix of total correlation function, $\omega$ is the matrix of intramolecular correlation functions, and $c$ is the matrix of direct correlation function. We use the inhomogeneous version of RISM since the fixed vertical lines breaks radial symmetry. The pair correlation function is $g(\mathbf{r})=h(\mathbf{r})+1$.

In our system, the correlation functions are random variables that depend on the arrangement of both chains. Suppose the chains, $l$, $l^{\prime}$, are parametrized so their sites have positions $\{p^{(m)}_{x},p^{(m)}_{y}\}$, $m\in\{l,l^{\prime}\}$. Since we want to find the correlation between just the $x$ coordinates of the lines, we will use the molecular correlation between $\{p^{(m)}_{x},0\}$ and the sites in the line as $\omega_{m,i}$, as opposed to the typical $\omega_{m,ij}$ correlation between sites $i$ and $j$. For two random chain configurations $\alpha=\{h^{(l)}_{1},\dots h^{(l)}_{N};h^{(l^{\prime})}_{1},\dots h^{(l^{\prime})}_{N}\}$ where $l,l^{\prime}$ index the two chains, this is

[TABLE]

In our model, where the chains do not directly interact in our range of interest of $h$, $C_{ll^{\prime},ij}(\mathbf{r})=0$. Inserting these functions in (14) and keeping only terms of order $\rho$, the total correlation function between $p_{x}^{(l)},p_{x}^{(l^{\prime})}$ is

[TABLE]

Hereafter, we abbreviate the vector positions in the $C$ functions as $\mathbf{r}_{l,j}$, $\mathbf{r}_{m,k}$. Note that $C_{ls}^{\alpha}$, $C_{sl^{\prime}}^{\alpha}$ are $N$ component vectors, one entry for each interaction site on the corresponding line.

Often, this is where the analytical part of the RISM procedure stops, the equations for $h$ and $C$ are evaluated using Picard iteration or some similar technique, leading to the radial correlation function. But our object of study is $u(r)\equiv\langle U^{(\alpha)}(r)\rangle_{\alpha}$, not the radial correlation function or potential of mean force for any specific random polymer arrangement. However, the potential of mean force for long chains should converge to the disorder averaged potential of mean force in the thermodynamic limit.

By definition, $-\beta U^{(\alpha)}(r)=\log\left(1+h^{(\alpha)}(r)\right)$. Recalling that $\rho$ is small, and averaging over the disorder,

[TABLE]

Note that the smallness of $h$ is an essential element in the simplification of this problem, since for most probability distributions, $\langle\log f(x,\omega)\rangle\neq\log\langle f(x,\omega)\rangle$.

Concentrating on the inner integral, which ignores edge effects, and noting that the two correlation functions depend on different and independent chain arrangements,

[TABLE]

Going from the first to second line is a consequence of the fact that if $L$ is large, or when we are using periodic boundary conditions, the $C$s will on average be independent of $y$, and the disorder of each chain can be averaged over independently, so the disorder averaged sum of $C$s is not a function of $y$.

We now have to determine the vector or direct correlation functions, $C^{(\alpha)}_{ls}(x)$. For hard spheres at low density, the direct correlation function is known to be $\Theta(2R-r)$ Tejero and Haro (2007), which makes physical sense given the nature of the hard sphere interaction. Due to the fact that our disordered chain sites’ can be close to each other and strongly correlated, the correlation function will differ from this. Physically speaking, $\sum_{i}C_{i}(\mathbf{r})$ should sum to $1$ in regions of space into which the hard sphere fluid cannot penetrate, and 0 where it can. Therefore, a physically realistic choice when the liquid is dilute is

[TABLE]

This evenly splits the responsibility of excluding liquid particles among all sites that exclude that region of space.

Since $\langle\sum_{j}C^{\alpha}_{ns}(x)\rangle$ just counts the probability that a random configuration excludes a hard sphere from being at a distance $x$ from the chain, we get

[TABLE]

Expanding the terms in the integral, we get that

[TABLE]

Note that the term in parenthesis is a constant (for fixed $\phi,\varphi$), so it does not matter in the averaged potential of mean force and we will omit it. Some simplification arise due to the fact that $p_{\phi,\varphi}(r)=1$ for $2\leq r$. We have, therefore, rederived our expression via RISM, which matches our expression for free energy (up to constants that do not depend on $h$), $u(x)=-\rho k_{B}T[x-\lambda(x)]L$.

Conclusion. Our two main theoretical results are the fluid density profile surrounding a disordered straight polymer, and depletion force between two such chains (eqn. (9, 11)). Additionally, we have numerically determined the depletion force between flexible polymers. While there are many studies on the interactions between polymers with definite structure functions, here we focused on disordered polymers, with random structure functions. We evaluated the average correlation between a disordered chain and a hard sphere fluid (i.e. the approach probability, $p_{\phi}$), and used this to find the depletion force between disordered chains in an approach similar to that of Asakura-Oosawa equations. We also showed how RISM can be extended, in the low density limit, to derive the average mean potential for this random ensemble.

It would be interesting to extend the random RISM method presented here to take into account fluids in the high density limit. There has been success in spin glass theory in dealing with averages over quenched disorder, e.g. the replica approach, which is also used in replica density functional theory Reich and Schmidt (2004); Schmidt (2005) and replica Ornstein-Zernike theory Pizio and Sokolowski (1997) to study fluids in random porous materials. A similar approach might be fruitful here.

Appendix A Calculation details

To evaluate (6) when $\xi(r)>1$ we exchange the order of the sum and integrals, and look at each term in the sum,

[TABLE]

But since our choice of which sphere was 0 is arbitrary, and we are using periodic boundary conditions (so there is circular symmetry), we must have that $T_{j}(r)=T_{k}(r)$ for all pairs $j,k$. In other words, all the $T_{j}$ are equal, call this common value $T(r)$. We will evaluate the one that requires the least work, $T(r)=T_{1}(r)$, since we can integrate all the integrals after the $\Xi$ function inductively.

[TABLE]

We can deal with the $\Xi$ function by adjusting the bounds of integration. Making the necessary adjustments,

[TABLE]

if $\xi(r)-1<N(1/\phi-1)$, and is [math] otherwise.

It is safe to exchange the lower bound of integration in (16), $x_{0}+2R$, with $x_{0}+2R\xi(y)$ since we are treating the case that $\xi(r)>1$. The second case comes from the fact that $\xi(r)-1\geq N(1/\phi-1)$ means that the lower bound of integration would have to be above the upper bound. The $T$ integral can be performed as follows, call $\alpha\equiv 2R\xi(r)$ and $\beta\equiv L-2R(N-1)$. Then

[TABLE]

where we made the substitution $u=\beta-x$. The integral is then just two polynomial integrals which can be performed with ease.

[TABLE]

if $\xi(r)-1<N(1/\phi-1)$, and is zero otherwise. This can be put into (15) and gives us our approach probability for $\xi(r)>1$, and can be combined with (7) to yield our complete expression for approach probability.

Appendix B Calculation without periodic boundary conditions

Here, we give some details about the approach probability for a chain, without periodic boundary conditions. In the thermodynamic limit, the probability is (9), but for finite $N,L$, there are corrections of order $\mathcal{O}(1/N)$, and the problem is more complicated. For simplicity, we assume that the spheres are of equal size (so $\kappa=2$), and that $\theta=0$.

Let $\Gamma^{(l)}_{N,R,L}$ be the measure of configurations of $N$ spheres of radius $R$ on a line of length $L$ such that no sphere extends “beyond” the line - that is the centers of the spheres must fall in the range $[R,L-R]$. We suppose that the line is along the x axis, from 0 to $L$. Using the same methods used to solve the analogous problem with periodic boundary conditions, we get $\Gamma^{(l)}=\frac{L^{N}}{N!}(1-\phi)^{N}$.

A sphere launched towards the chain can collide if it has x coordinate in the range $[-R,L+R]$, so we consider our launched spheres to be chosen uniformly at random from this range. A sphere with x coordinate less than $x_{1}$ can only interact with the first sphere, likewise, a sphere with x coordiate greater than $x_{N}$ can only interact with the last sphere.

[TABLE]

In the case where $\gamma<1$, the $L(r)$, $R(r)$ integrals can be evaluated (they have the same value, by symmetry)

[TABLE]

and the sums telescopes and so can be evaluated,

[TABLE]

For the $\gamma>1$ case, all the $S_{j}$ are equal, like before, and $L(r)=R(r)$ by symmetry. By adjusting the integration bounds in the integral, just as in (18), we can evaluate $S(r)$.

[TABLE]

Putting all this together, we find that the approach probability to a chain without periodic boundary conditions is

[TABLE]

While this equation is much more complex than its periodic boundary condition counterpart, they both have the same limiting value in the thermodynamic limit.

Appendix C Free energy in three dimensions

For a pair of disordered chains in three dimensions, the procedure for finding the free volume is similar to the case of disordered parallel chains. Suppose the fluid radii, and the radii of the monomers are both $R$. Let $V_{A}$ be the volume within $2R$ of the first chain, $V_{B}$ be the volume within $2R$ of the second chain, and $V_{\text{int}}\equiv V_{A}\cap V_{B}$. The expected free volume is now a complicated function of the positions of the centers of the chains, $\mathbf{r}_{1},\mathbf{r}_{2}$, and normal vectors describing the orientations of the chains, $\hat{\mathbf{n}}_{1}$, $\hat{\mathbf{n}}_{2}$.

[TABLE]

The function $p_{\phi}$ is the familiar probability of approach from before. Since the partition function is $Z(\{\mathbf{r},\hat{\mathbf{n}}\})=V_{E}(\{\mathbf{r},\hat{\mathbf{n}}\})^{N}/(N!\Lambda^{N})$, and $F=-k_{B}T\log Z$, we take the log and expand in terms of $1/V$,

[TABLE]

where $F_{0}$ is a term that does not depend on the configuration of the cylinders.

Unfortunately, after this point, further simplification, becomes impossible. To find the depletion force on the chains, we evaluate the gradient of the free energy,

[TABLE]

The torque on the chains could be evaluated similarly. Note that the volumes over which we are integrating depend on the orientation of the particles. Note that $F$ and the magnitude of $\mathcal{F}$ only actually depend on the relative locations and orientations of the two chains.

As a final word of caution, this formula is only valid when the minimum distance between the lines is $\geq 2R$, which is necessary to ensure that the arrangement of spheres on the lines are independent from one another. Furthermore, edge effects are not included, so this formula assumes that the parts of the lines that are close to one another are away from the ends of the lines.

Appendix D Simulation details

We give some further details concerning the molecular dynamics simulations used to obtain data for Fig.s 3, 4, 5. Straight polymers, the objects that we also treat theoretically in this paper, are generated by first specifying the length of the “linear backbone” of the polymer (which is not explicitly represented by objects in the simulation, but helps us keep track of where the monomers should go), $L$, and a target linear volume fraction, $\phi=2NR/L$. Only certain volume fractions are possible, the ones that correspond to integral numbers of disks on the backbone. The total number of spheres, $N$, is calculated from $\phi,\,R=0.05$.

The positions of the disks on the backbone are chosen by picking $N$ numbers uniform at random from $[0,L-2RN]$ and sorting them so $x_{1}<...<x_{N}$. The positions of the disks are set to $r_{1}\equiv x_{1}+R$, $r_{2}\equiv r_{1}+x_{2}-x_{1}+2R$, … , $r_{N}=r_{N-1}+x_{N}-x_{N-1}+2R$. This method samples uniformly at random from the set of all configurations of disks on the backbone with the constraint that no disks overlap.

To create a flexible polymer, as we do to obtain data for Fig. 3 and 6, we modify the above algorithm slightly. We pick a small monomer radius, $R_{c}=0.02$, for the disks that act as chain links in the polymer. The chain disks do not interact with any other particles except via harmonic bond and angle forces with their (at most) two adjacent neighbors in the chain. Again, based on the desired chain length and linear volume fraction, the total number of (large) disks is calculated. Positions for the large disks are calculated as before, but since we are filling the space between large disks with chain monomers, we calculate the expected number of chain particles between the two large disks, $n_{x}=(r_{k}-r_{k-1})/(2R_{c})+1$ and place $n_{c}=\text{floor}(n_{x})+\hat{C}$ chain particles between adjacent large disks. The random variable $\hat{C}$ is $1$ with probability $(n_{x}-n_{c})\in[0,1]$ and 0 otherwise. This compensates on average for the fact that we always round $n_{x}$ down to get the deterministic part of $n_{c}$.

As mentioned before, the hard sphere interaction was approximated by using a harmonic repulsion force with a very large stiffness. This force was used for chain-chain, chain-solvent, and solvent-solvent forces,

[TABLE]

where $R_{1},\,R_{2}$ are the radii of the interacting spheres, $\kappa$ is the repulsion constant, and $\Theta$ is the Heaviside step function. In our simulations, the radii of all the interacting particles (the large monomers and the hard sphere fluid) are the same, $R=0.05$, each particle had density $1$ (and therefore mass $m=0.00785$), and the spring stiffness used was $\kappa=500$. The period of collision was therefore about $t_{c}=0.025$. In all simulations, periodic boundary conditions were used, and time evolution is carried out using Velocity-Verlet integration, modified to run a Nose-Hoover thermostat Evans and Holian (1985) to keep the system at constant temperature. To evaluate non-bonded forces, a linked cells structure is used to create Verlet list for force evaluation.

The straight polymers are each rigid disordered chains, created as described above, which we constrain to only move in the x direction, and remain vertical, all net torque and force in the $y$ direction are set to zero. The right chain was displaced a random amount in the $y$ direction, uniformly chosen to be between 0 and $2R$ above the start of the left chain. This insures that measurements are not effected by systematic correlations between the spheres on the two chains, which can occur at high $\phi$. At distances greater than $5R$, a harmonic force was activated to prevent the chains from departing far away thereby allowing us to collect larger amount of data.

For the straight polymer, each $\phi$ curve was composed using data from 18 simulations, 6 where the chain were in square boxes of solvent of side length $160R$, with 408 solvent particles, and 12 where the chains were in boxes of solvent of side length $320R$, with 1630 solvent particles. In both cases, the length of the lines were $32R$, and the volume density of the solvent was $\phi=\frac{\pi R^{2}}{V}=0.05$.

The pressure for the straight chains is calculated by binning the total force on each chain, projected in the $\pm\hat{x}$ direction, towards the other chain. Each bin is averaged, and the force is divided by the length of the line, $L$.

As alluded to in the main text, calculating the pressure for the flexible chains is more complicated since the orientation of each chain segment is different, and the chains are no longer parallel. For each monomer (disk) in each chain, we calculate the nearest point on the flexible backbone of the other chain, which may lie between two monomers. If this is the case, we assume that the chain is straight between monomers. We then bin the force on the monomer projected in this direction, $\hat{n}$, divided by $2R$ (this is the pressure in the $\hat{n}$ direction on the particle), binning by the distance to that nearest point on the other chain.

Since the flexible polymers are not constrained to be near one another, have many more internal degrees of freedom, and because we use smaller simulations than in the linear polymer case, we average the system much more that in the straight polymer case to obtain accurate data. Each $\phi$ curve was composed of 150 simulations run for 30,000 seconds each. Each simulation took place in a box of solvent with side length $100R$ with 319 solvent particles. The length of the chain was $40R$, and the density of the solvent was $\phi=0.1$.

The GFlow molecular dynamics package is freely available at https://github.com/nrupprecht/GFlow.

Appendix E Numerical approximation of entropic pressure

We have use Schmidt M. (2014) to find an approximate closed form expression for (11). For $\phi=\varphi$,

[TABLE]

in units of $\rho k_{B}T$. With $c_{1}=-0.158,c_{2}=0.031,c_{3}=3.667,c_{4}=0.899$ and $c_{5}=0.236$ the mean error is $0.5\%$.