An analytic Percus–Yevick approach to the RISM model of monodisperse, homopolymer melts

The reference interaction site model is often used as a nonperturbative description of the equilibrium structure of monodisperse polymer melts on nanoscopical scales. Assuming that the single polymer chains in the melt have ideal (Gaussian) conformations, we derive for the melt a set of generalized Percus–Yevick (PY) equations which can be exactly reduced to a simple algebraic systems of equations. While in the case of the simple atomic sphere systems this involves cubic polynomials, in the analogous hard-sphere polymer system this leads to quintic polynomials. In the limit of infinitely long chains the exact, generalized PY static structure factors, the compressibility and the direct correlation functions are found as functions of the monomer density.


Introduction
The behaviour of a polymer melt differs strongly from polymer solutions. Whereas in recent years significant progress has been made in the understanding of polymer solutions, especially by using scaling arguments and renormalization group techniques [1], the microscopic description of polymer melts is still less clear. The reason for this difference are the relevant length scales, which in polymer solutions are of the order of the radius of gyration, whereas in dense polymer systems the short-range interactions remain effective, with a characteristic length scale of the order of the polymer segment size.
These intermolecular interactions are very strong and cannot be treated using the usual perturbation theories for polymer solutions. Numerical Monte Carlo and molecular dynamics simulations of both lattice and continuum models [2]- [6] have been limited to small systems (i.e. to systems containing relatively short chains). Furthermore, the numerical investigation of such systems have to be performed on a time scale longer than the longest relaxation time, which determines the evolution of a given initial state to thermodynamic equilibrium. This relaxation time increases rapidly with increasing chain length and density [7].
A successful alternative approach is given by a microscopic equilibrium theory of polymer melts based on non-perturbative statistical theories. Such theories are characterized by nonlinear, self-consistent equations of the Ornstein-Zernike type, such as the Percus-Yevick (PY) equation [9] and the hypernetted chains equation [10], which have previously been used for simple atomic fluids. The generalization of the PY equation to one-component molecular liquids [11,12] leads to the reference interaction site model (RISM). First analytical and numerical approaches to the equilibrium structure of dense polymer systems using the RISM theory have been given by Schweizer and Curro [13]- [15]; extensions to polyelectrolytes are given by Hoffman et al [16]. A very nice theoretical solution of the RISM problem was given recently by Fuchs [17]. Since most of these models collect the action of chemically connected atoms into some monomers with an idealized hard spherical feature, the expected results may be reasonable at sufficiently large scales in comparison with the typical atomic diameter. Thus, the RISM theory for polymer melts is mainly a nanoscopical theory.
It should be remarked that the RISM theory is a diagrammatically improper theory, because the distinction of different topological classes of diagrams is neglected. Thus, this theory has been corrected by additional contributions [18,19]. The analysis of the RISM theory is based on a closed system of nonlinear equations which are somewhat simpler and probably of similar accuracy as such proper theories [20]. Thus, the RISM theory can be used as an approximation which allows a reasonable description of the static structure of molecular systems. Therefore, the RISM theory can be very helpful for a discussion of dense polymer systems down to nanoscopic scales; in particular, the behaviour of the static structure factor for long wavelengths and its dependence on the chain lengths and the monomer density are still unsolved problems. Quantitative statements are very important for an understanding of the glass transition of a polymer melt in terms of the mode coupling theory [21], where long-range density fluctuations have a considerable effect on the glass kinetics.
The purpose of the present paper is to determine an analytical approximation (i.e. a system of simple algebraic equations, which can be solved by standard methods) of the RISM theory for a dense system of identical polymer chains with a well-defined chain length L of connected hard core monomers and randomly choosen interaction constants.

Static structure factor and the RISM theory of chains
The correlation function g(r − r ) is given by the equilibrium ensemble average of δ(r − r i ) δ(r − r j ), performed over all position configurations (1) with the particle density ρ (i.e. in the present case, ρ represents the number of monomers per unit volume). On the other hand, the static structure factor is defined by where N is the monomer particle number. Hence, we have the well-known relation between the correlation function and the static structure factor. It is useful to split the index i in (1) and (2) into two subindices in the case of polymer chains: each molecule (index I) contains L monomers (or effective monomer segments), determined by the index α, i.e. a unit of the system is given by i = (I, α), with I ∈ [1, M] (M is the total number of chains in the system) and α ∈ [1, L] (N = ML). Therefore, the static structure factor becomes S(k) = ω(k) + ρh(k).
Here ω(k) is the static structure factor of a single chain and h(k) is the Fourier-transformed intermolecular correlation function with which involves all correlations between different chains. The generalization of the Ornstein-Zernike equation [22] to molecular fluids leads to a connection between the functions ω αβ and h αβ : where C γδ (k) is the Fourier-transformed direct correlation function and ρ c the chain density (chains per unit volume). This relation is exact but void as long as we have no second independent equation or boundary condition to determine the direct correlation function. Following the usual arguments, which lead to the PY equation for atomic fluids [23] (hard-core interaction, shortrange direct correlation function), we set h αβ (r) = −1 for |r| < d 0 , where d 0 is the closest approach distance between monomers. Both equations (5) and (6) constitute the RISM theory [13]. The neglect of chain end effects (this approximation is allowed for chain lengths sufficiently large in comparison with the size of an effective monomer) leads to translation symmetry along the chains and therefore to the useful simplification C αβ (r) = C(r), i.e. C(r) becomes independent of the specific position of the monomers along the chains. Taking this simplification into account in (3) and (4), the Ornstein-Zernike equation (5) becomes Note that the chain density ρ c is related to the monomer density ρ by ρ = Lρ c . Equation (7) together with the closure relations h(r) = −1 for |r| < d 0 and C(r) = 0 for |r| > d 0 constitute the RISM theory of a dense system of polymer chains and was treated by different analytical approximations [13] and numerical approaches (based on variational principles) [14,15]. However, a complete analytical approach to the RISM theory of polymer chains will be very helpful for a further discussion of the static structure factor of dense polymer systems. Especially, such a solution can be used as a reference model for a further perturbation theory. Another important reason to search for an analytical representation consists in the investigation of the influence of the structure control parameters, such as the monomer density and chain length. Therefore, an analytical approach to the solution of (7) will be determined and studied in the subsequent sections.

The integro-differential equation representation
It is well known from neutron scattering experiments probing a few deuterated chains in a melt of identical, but hydrogenated, chains [24]- [26] that the end-to-end distances of polymer chains are distributed in a Gaussian manner (screening effect). In general, both the static structure factors of the chains and of the whole polymer system are a result of the inter-and intramolecular interaction, i.e. both depend on density and chain length. However, this 'self-organization', i.e. the simultaneous formation of the structure of the whole system and the configuration of single chains in the system, observed in many numerical simulations on dense polymer systems [29,30], cannot be accounted for by the RISM theory.
The RISM theory represents only an analytical connection between the static structure factors of the molecules and of the total system of the molecules (see equation (7)). It is necessary to introduce the static structure factor of a single chain in the melt as an input of the theory. Therefore, it is reasonable to use the Gaussian distribution (which is confirmed by neutron scattering) for the chain configuration at all densities, although this approximation is exact only for dense systems. Such a chain can be represented by L eff effective segments. The length of the effective segments 0 corresponds to the effective number of segments via L eff 0 = Ld 0 .
Thus, according to (3), ω αβ (k) is the Fourier-transformed Gaussian probability distribution for two monomers at distance |α − β| along the chain. We assume in the following that the length 0 of the effective segments is in a fixed ratio to the repulsion distance d 0 of the hard-core interaction, i.e. 0 = κd 0 . Thus, the static structure factor of a single Gaussian chain is given by ω αβ (k) = exp{− 1 6 k 2 |α − β| 2 0 }. After performing the summation we get with the radius of gyration R 2 g = 1 6 L eff 2 0 = 1 6 Ld 0 0 and the Debye function [27] For sufficiently long chains and a slow variation of the disorder of the monomer composition along the chain we can use the approximation D(x) = (x + θ) −1 , which holds both for large and small values of x (see appendix A). The deviation of this approximation from the original Debye function is at most 10%. Hence (9) and (7), we get or after transformation to position space with the propagator This integral equation is equivalent to the integro-differential equation (see appendix A): with D n = [ξ n − ∇ 2 ]. Note that (11) and (13) (13) in the limit ξ 1 /ξ 2 → 1.

Laplace-transformed integro-differential equation
We expect isotropic solutions h(r) = h(r) and C(r) = C(r) of (11) because of the homogeneity and isotropy of the underlying statistical distribution function of the polymer system in the real space. Therefore, we write Using the new functions and the new closure relations ϕ(r) = 0 for |r| < d 0 and γ(r) = 0 for |r| > d 0 , we get The one-sided Fourier transformation of the reduced one-dimensional integro-differential equation leads to (see appendix B) with Therefore, we have the formal solution Using (20) we can construct a function (z), given by which we will analyse in the next section.

The analytical properties of (z)
We first note that the function (z) is an even function, i.e. (z) = (−z). Using (19) and the closure relation γ(r) = 0 for r > d 0 , we obtain (z = x + iy), i.e. (z) has no singularity anywhere in the finite region of the complex plane. Moreover, (z) behaves as (z) ∼ z −1 in the upper half-plane, whereas (z) grows exponentially as |z| → ∞ in the lower half-plane.
In the same way it follows that (z) has no singularity in the upper half-plane, but a double pole at z = 0 and possibly further poles in the lower half-plane. However, it follows from the special construction of (21) that (z) has no singularity in the whole finite upper half-plane (i.e. without infinity). Because of the z 4 -term, (z) is also regular on the real axis. Moreover, using the symmetry of (z) under the inversion z → −z, there exist no singularities in the finite lower half-plane. Hence, (z) is an even entire function and can be represented by a polynomial of degree 2K (theorem of Rouché) with polynomial coefficients θ n . The degree of the polynomial is determined by the behaviour of (z) as |z| → ∞. Substituting (20) in (21) and performing the limit |z| → ∞, both denominator and numerator are controlled in the upper half-plane by the divergence of (−z), which grows as exp{yd 0 }. Thus, (z) becomes for sufficient large |z| Because of (22), (z) behaves in the upper half-plane in the limit |z| → ∞ as (z) z −1 . Hence, the highest power of (z) is determined by Thus, (z) must be an even polynomial of degree 8 (i.e. K = 4). The coefficients of this polynomial can be obtained from the behaviour of (z) around z = 0. We write the polynomial (z) with another set of coefficients θ n : The substitution of (20) into (21) followed by expansion of this expression in powers of z around z = 0 leads to an implicit representation of the coefficients θ n (n = 0, . . . , 4); see appendix D.
Note that all higher coefficients of the Taylor expansion around z = 0 must vanish identically as a consequence of the theorem of Rouché.

The solution of the integro-differential equation
The combination of (21) and (20) gives now Here and in the following we understand (z) as the even polynomial (26) of degree 8. The function Q(z) can be written as After multiplication of (27) with exp{−izr} (r < d 0 ) and integration over z from −∞ to +∞ along a contour C (with an infinitesimal distance to the real axis) in the upper half-plane, we have the following situation: 1. The function (−z) diverges in the upper half-plane as exp(−izd 0 )/z (see equation (22)). Hence, because of (20), the product (z) exp{−izr} ∼ exp{iz(d 0 − r)} converges in the upper half-plane exponentially to 0 for |z| → ∞, i.e. the contour C can be closed by a semicircle S + of infinite radius in the upper half-plane for all terms of (27), which contain the function (z). 2. The function (z) converges for infinite large |z| in the upper half-plane ( (z) ∼ z −1 ), i.e. the product (z) (z) exp{−izr} ∼ exp{iz(d 0 − r)}/z converges also for |z| → ∞ and 0 r < d 0 . Thus, the integration contour C of the corresponding term in (27) can also be closed by an infinite semicircle in the upper half-plane. 3. The integral 4. The products of exp{−izr} with terms, which contains neither (z) nor (z), are divergent in the upper half-plane and converge in the lower half-plane for r > 0. Therefore, the contour C must be closed by an infinite semicircle S − in the lower half-plane.
Taking into account the orientation of the closed integration contours (i.e. the integrals in the lower half-plane change their sign), (27) becomes Thus, after some algebraic transformations (see appendix C) we obtain the solution with the coefficients λ m , λ m and χ m listed in appendix C. Equation (30) is the solution of (16) and therefore of (13) and (11). The coefficients λ m , λ m and χ m depend on θ n and the density ρ only. Furthermore, the θ n are obtained as solution of the system of nonlinear equations, which are obtained by substituting (30) into (D.6) and inserting these expressions into (D.1)-(D.5), respectively.

Discussion
The structure of (30) and the system of nonlinear equations (D.1)-(D.5) reduces the determination of the static structure factor of a polymer system to a purely algebraic problem, i.e. the solution is possible by using standard methods. Note that the system of nonlinear equations can be simplified to a linear one, if ξ 1 /ξ 2 → 1. This case corresponds to a system of hard spheres (i.e. the RISM equation (7) is reduced to the PY equation [9]) with a well-known solution [28].
On the other hand, we are interested in sufficiently long (or infinitely long) polymer chains when ξ 1 /ξ 2 → 0. This limit leads to another simplification of the system of nonlinear equations, which will be discussed here.  which are linear in the coefficients θ n and the nonlinear equation which is of second order in θ n . Note that the n are linear functions of the coefficients θ n and β is the reduced density β = (π/144 √ 3)ρd 3 0 . These equations can be solved by simple methods (solution of the linear part (31), substitution of these solutions in the reduced nonlinear equation (32), which leads to a simple quadratic equation). Finally, the coefficients θ n can be expressed in terms of the density ρ. Figure 1 shows the behaviour of these coefficients.
The direct correlation function C(r) can be determined by using (30). Unlike the behaviour of a system of hard-core spheres, C(r) diverges for r → 0. This is a characteristic feature of all systems of polymer chains, not only of infinite long chains. The structure of C(r) is shown in figure 2.
The Fourier-transformed direct correlation function is given by Straightforwardly, one can determine the static structure factor which is presented for different densities ρ in figure 3. One can see that the free-chain contributions (especially the divergent behaviour for k → 0) vanish with increasing density. The static structure factor behaves like a system of hard spheres for sufficiently large densities. This expected effect is well known and gives a first verification of the possible validity of the RISM theory of polymer chains.  Another remarkable point is the existence of an upper density ρ u (highest allowed density). For a system of hard spheres, one obtains πρ u d 3 0 /6 = 1, i.e. ρ u corresponds to ideal dense cubic packing. Note that there is no free volume, i.e. the density ρ u is much larger than the density of any possible crystalline lattice. This characteristic effect can also be observed for the RISM theory of polymer chains. Here, one finds πρ u d 3 0 /6 ≈ 2.18485 (under the condition that the length of the effective segments 0 is equal to the repulsion distance d 0 = 1 of the hard-core interaction). The increase of the upper density (in comparison with the hard-sphere result) is connected with the structure of the polymer chain. Because of the used structure factor ω(k) of the single chains, each chain is a random (but not self-avoiding) walk. The self-avoiding condition is valid for the interaction between different chains only. Thus, the overlaps of the self-penetrating chains lead to a latent free volume and consequently to an increase in ρ u .
The thermodynamic equation of state can also derived in the RISM approach. Here, we use the well-known relation [23] (compressibility equation) (ρ c is the density of thermodynamic objects, i.e. the density of polymer chains in the present case). Using ω(0) = ξ 2 /ξ 1 = 1 + L ≈ L (for L → ∞), the static structure factor (34) is given by Hence, for L → ∞  Figure 4 shows the inverse compressibility (∂p/∂ρ) as a function of the density. As expected, there exists a divergent behaviour as ρ → ρ u , i.e. the compressibility vanishes if the density ρ reaches the maximum of the allowed density ρ u of this theory. In figure 5 we show the comparison of our analytical expressions given by the static structure factor with Monte Carlo computations [30] at various densities. In general, this comparison is good, at least over the range of numerical computations.
In conclusion, we find that a reasonable approach of the RISM theory for linear polymer chains is possible, which allows the discussion of the influence of density and chain length on the structure of the system. Such a solution is represented by the algebraic system (D.1)-(D.5). From this point of view, one has an excellent reference model for various perturbation theories, which consider the self-avoiding structure of the polymer chains. (In the present paper, each chain forms a random walk, the self-avoiding condition exists for the interaction between the monomers of different chains only.) Such an investigation will be discussed in a subsequent paper. Our results (for the static structure factor, compressibility) are purely non-numerical expressions for dense polymer systems. These hard-core segment expressions can also be used for a high-temperature perturbation theory for segments possessing soft potential contributions with a hard-core segment polymer reference fluid as used in [31,32].  for a numerical volume fraction 0.58. The curve represents analytic PRISM theory with πd 3 0 ρ/6 = 1.5 and squares represent the numerical MC computations.

Appendix A. Properties of the RISM equation of polymers
The functions h(r), C(r) and J(r) = d 3 r h(r − r )C(r ) have the following asymptotic properties: • The direct correlation function C(r) vanishes for |r| > d 0 .
Hence, the functions h(r), C(r) and J(r) = d 3 r h(r − r )C(r ), respectively, vanish for |r| → ∞. These properties allow the formulation of (11) as an integro-differential equation. After introducing the differential operators we obtain from (12), Therefore, the application of ( D 1 D 2 ) 2 on (11) leads to the integro-differential equation It should be remarked that the order of the propagators G n and the functions h and C is irrelevant in each term of (11). This is a general property of such convolution integrals. Usually, the application of a differential operator to an arbitrary equation leads to a new equation with a larger manifold of solutions. Thus, it is necessary to prove that the set of solutions increases in course of the mapping from the integral equation (11) to the integro-differential equation (A.3).
To do this, we analyse the application of a convolution integral to the differential operator D n , F(r) is acting for h(r), C(r) and J(r), respectively, i.e. the function F(r) vanishes asymptotically for |r| → ∞. Using Stokes's formulae, it follows Note that the propagators G n also vanish for infinitely large distances, i.e. the differential operator D n and the convolution procedure G m * can be exchanged without restrictions, Thus, the application of the convolution procedure G 1 * G 1 * G 2 * G 2 * to (13) subject to the commutation relation (A.4) and the identity (A.2) leads to (11). The clear mapping from (11) to (13) by the application of differential operators and its inverse from (13) to (11) by the application of the convolution implies that both equations have an identical set of solutions. Consequently, the integro-differential equation (13) contains no additional degrees of freedom in comparison with the original integral equation (11). The real cause for this equivalence is the asymptotic behaviour of the functions C, h and J = C * h as |r| → ∞.

Appendix B. Laplace transformation of the integro-differential equation
We use the imaginary Laplace transformation (one-sided Fourier transformation), defined by (B.1) The Laplace transformed expression of the first term of (16) is given by The polynomial P(iz) contains the initial conditions with γ (n) 0 = ∂ n /∂r n (Cr) | r=0 . In general, (16) becomes Here, h * C means the convolution integral whereas the third-order polynomial P G (iz) contains all initial conditions of the three terms of (16). This polynomial can be interpreted as the Laplace-transformed representation of δ-functions and their derivatives. In other words, the initial conditions suggest the existence of a δ-like source term at r = 0. On the other hand, such source terms do not exist for the original differential equation (13). This means, the combinations of all initial conditions within the polynomial P G (iz) vanish by internal cancellations, i.e. P G (iz) ≡ 0. Note that the appearance of initial conditions is only a seeming effect, which is connected with the special choice of the spherical coordinates and the succeeding reduction to the radial coordinate. Furthermore, the convolution integral can be written as The Laplace transformation is most easily evaluated by performing the r integration first, Here, we have used the fact that for r < t the integral kernel vanishes as a consequence of the closure conditions. After integration over r and using (21) we get Considering the definitions, we get (B.10)

Appendix C. Determination of the direct correlation function
From our previous discussion, (z) and (z) have no singularities in the upper half-plane. Therefore, all integrals in (29), which contain these functions, can be simplified by using Cauchy's