An amended approximation of the non-Gaussian probability distribution function

1. Introduction

A statistical theory of a single long-chain polymer molecule is an essential basis for developing a micromechanical model for rubbers. Traditional molecular simulations take into account all chemical aspects at the molecular level, and are therefore very powerful and comprehensive [1]. However, this detailed information is essential only in microscopic dimensions or small-time ranges of up to 1 ms. Due to the computational expenses and limitations to size and calculation time of such simulations, statistical averaging to larger scales is of interest. The majority of the rubber elasticity network models are based on the statistical theory of an idealized chain and in particular on the freely jointed chain (FJC) model. It is the most simplified representation of a real polymer molecule that can randomly fluctuate in space.

In polymer physics, the probability distribution function (PDF) of an FJC with fixed end positions can be described by the random walk theory. The random walk theory describes a random movement of a particle where a direction in each step is independent of all previous steps. If the length of each step is constant, the particle trajectory is one of all possible conformations of the FJC model. The conformational entropy s of the chain is directly related to the PDF and, thus, to the end-to-end distance R by the Boltzmann equation

s = k B ln [ P ( R ) ] ,

(1)

where k B denotes the Boltzmann constant. The free energy W of an ideal polymer chain is derived from its conformational entropy, as

W = − Ts ,

(2)

where T denotes the absolute temperature. Thus, from a point of view of constitutive modeling, finding an approximation for P ( R ) that leads to a minimized error in ln [ P ( R ) ] is of larger interest than minimizing the approximation error of P ( R ) itself.

The first application of the random walk theory to polymer chain statistics was independently proposed by Kuhn [2] and Guth [3]. Kuhn and Grün [4] introduced the inverse Langevin approximation for single-chain statistics for large deformations which appears accurate for relatively long polymer chains. Many constitutive models of rubber elasticity incorporated the Kuhn–Grün approximations of the non-Gaussian theory [5–8]. These are, for example, the three-chain [9], the four-chain [10], the eight-chain [5], the micro-macro unit sphere model [7], and full-network models [8,11,12]. However, deviations of the Kuhn–Grün model from the exact distribution function in the cases of large extensibility and short chains have motivated to search for more accurate alternatives [9,13,14]. Moreover, the inverse Langevin function cannot be represented in an explicit form and requires an approximation, which increases the deviation from the exact distribution function.

As a result, all models proposed so far suffer from the inaccuracies in the Kuhn–Grün approximation. While some approximations outperform the others in the majority of cases, there is no explicitly superior approximation. The ideal choice of the approximation function depends on a variety of factors, such as the chain parameters, or the desired application of the approximation (the probability distribution function itself, the corresponding free energy, or the entropic force developed by the chain). In addition, when making a choice regarding the best approximation for a given application, the applied error measure plays a significant role since the approximation that grants, for example, the best maximal relative error is not necessarily the same that provides the best mean absolute error.

It is hence reasonable to choose an approximation that suits the error measure deemed most critical for the desired application. In this regard, it appears questionable to take one specific approximation as the final solution that fits all conditions. All previously proposed approximations perform well only for a particular number of chain segments while the effect of n remains undiscussed. In this paper, we develop an approach to derive a family of approximations for the free energy of a polymer chain with a broader range of chain segments.

2. Statistics of a single chain

Chain-like molecules such as a single polymer chain in a dilute polymer solution consist of relatively small repeating units (monomers). They are covalently bonded together at some position referred to as cross-linking junctions. In an equilibrium state, the polymer chain can take a definite number of configurations between two end points.

By ignoring the variation of the length of a valence bond from its mean value due to atomic vibrations in the chain molecules [15], one can assume that each bond has a constant length l . The distribution function P ( R ) of the end-to-end distance vector R subject to the following condition (see Figure 1)

R = ∑ j = 1 n r j ,

(3)

can be expressed as [16]

P ( R ) = 1 2 π 2 ‖ R ‖ ∫ 0 ∞ ρ sin ( ‖ R ‖ ρ ) ( sin ( ρ l ) ρ l ) n d ρ .

(4)

OPEN IN VIEWER

Figure 1.

Configuration of a polymer chain with n segments and the first segment placed at the origin of a coordinate system.

If the mean end-to-end distance | | R | | = R · R is much less than the maximum extension of the chain ( | | R | | << nl ) , equation (4) reduces to the Gaussian distribution function

P ( R ) = ( 3 2 π n 〈 l 2 〉 ) 3 / 2 exp ( − 3 ‖ R ‖ 2 2 n 〈 l 2 〉 ) .

(5)

Most of the problems in polymer physics are treated based on the Gaussian chain model (5), which is physically unreasonable for real chains [15] in particular in the case of large deformations.

The probability distribution function is subject to the normalization condition generally written by [17]

∫ − ∞ ∞ ∫ − ∞ ∞ ∫ − ∞ ∞ P ( x , y , z ) dx dy dz = ∫ 0 ∞ 4 π r 2 P ( R ) d r = 1 .

(6)

Equation (4) is given in the integral form and is thus not convenient for practical applications and numerical computation, particularly for large values of n . Rayleigh [18] obtained a solution for small values of n , e.g., n = 6 . Chandrasekhar [19] evaluated the integral equation (4) for any finite value of n separately. Nevertheless, the calculations become very tedious for larger values of n . There have been several attempts to get an exact solution for equation (4). Treloar [20], Wang and Guth [9], Nagai [21], and Hsiung et al. [22] represented it in form of a series as

P ( R ) = 1 2 n + 1 ( n − 2 ) ! π l 2 ‖ R ‖ ∑ p = 0 k ≤ ( n − | | R | | / l ) / 2 ( − 1 ) k ( n k ) ( n − 2 k − ‖ R ‖ l ) n − 2 ,

(7)

which is also inconvenient for practical calculation, especially for large values of n . For this reason, it is usually approximated by a closed-form expression. In the following section, the most well-known approximations of the exact probability distribution function are discussed and compared with respect to their accuracy.

3. Appraisal and comparison of different approximations

In order to avoid tedious and computationally expensive calculation of the integral form (4) or its series expansion (7), several approximations based on different methods and approaches have been proposed. It is clear that a robust and comprehensive approximation must be capable of fitting to the exact solution for all ranges of the variables | | R | | , l , and n or at least to their physical ranges for a specific application. In order to verify the effects of these parameters, we introduce an end-to-end extension ratio t ≡ R / nl such that 0 ≤ t ≤ 1 , and where the abbreviation R = | | R | | is used. In the following, different methods approximating the exact solution are discussed and compared.

3.1. Asymptotic analysis

Based on the asymptotic analysis [23], Kuhn and Grün [4], James and Guth [24], and Flory and Rehner [10] obtained an approximation for equation (4) for the case n >> 1 and t << 1 . Their results are the same but follow distinctive approaches resulting from the non-Gaussian distribution with limited chain extension. The Kuhn–Grün approximation is given by

P KG ( t ) = 1 ( 2 π n l 2 ) 3 / 2 ( sinh L − 1 ( t ) L − 1 ( t ) exp ( t L − 1 ( t ) ) ) n ,

(8)

where L − 1 ( t ) represents the inverse of the Langevin function defined as L ( x ) = coth ( x ) − 1 / x . An accurate approximation of this function can be based on a high-order series expansion (see, e.g., Itskov et al. [25]), on a non-rational or a rational function. Recently, a number of approximations have been proposed for it (26–28). In the limit of n → ∞ , equation (8) results in the Gaussian distribution (5).

An amended Kuhn–Grün approximation proposed by Jernigan and Flory [14] aims to describe polymeric network chains in the context of rubber elasticity at high extensions t ≈ 1 as

P aKG ( t ) = C aKG P KG ( t ) ( L − 1 ( t ) t ) ,

(9)

where C aKG denotes the normalization factor satisfying equation (6).

3.2. Steepest descent method

The steepest descent or saddle-point method is an extension of Laplace’s method for an integral approximation. By the steepest decent, Wang and Guth [9] obtained the solution of equation (4) for n >> 1 as

P WG ( t ) = P KG ( t ) ( L − 1 ( t ) t k 2 − 1 2 ) ( 1 + q ( t ) + ⋯ ) ,

(10)

with q ( t ) = ( k 4 − 3 nk 2 2 ) / 8 k 2 2 , where k 2 and k 4 are given by

k 2 = 1 − t 2 − 2 t L − 1 ( t ) k 4 = 3 nk 2 2 − ( 5 − 14 t 2 + 9 t 4 − 22 t L − 1 ( t ) + 30 t 3 L − 1 ( t ) − 2 ( L − 1 ( t ) ) 2 + 30 t 2 ( L − 1 ( t ) ) 2 + 6 t ( L − 1 ( t ) ) 3 ) .

(11)

In a similar way, Yamakawa [16] also derived the following expression

P SD ( t ) = P KG ( t ) ( L − 1 ( t ) ) 2 t ( 1 − ( L − 1 ( t ) csch L − 1 ( t ) ) 2 ) 1 / 2 ,

(12)

where csch represents the cosecans hyperbolicus. It should be noted that equation (12) is valid over the whole range of 0 ≤ t ≤ 1 .

3.3. Series expansions

Most of the proposed approximations focus mainly on the simplified and specific solutions of the exact problem like series for the inverse of Langevin function. Most of the Taylor series expansions are developed around t = 0 and are only valid for limited extensions of the polymer chain. In case of n >> 1 (very long chains), the Taylor series for ( sin x / x ) n is written as

( sin x x ) n = exp ( − n ∑ k = 1 ∞ B 2 k − 1 ( 2 x ) 2 k 2 k ( 2 k ) ! ) = exp ( − nx 6 ) ( 1 − n x 4 180 − n x 6 2835 + 7 n 2 − 12 n 453600 + ⋯ ) ,

(13)

where B k denote the Bernoulli numbers. Using equation (13), equation (4) becomes

P series ( R ) = 3 4 π n l 3 ( 6 π n ) 1 / 2 exp ( − 3 R 2 2 n ) { 1 − 3 ( 5 − 10 R 2 n + 3 R 4 n 2 ) 20 n + 1015 − 9660 R 2 n + 13734 R 4 n 2 − 5364 R 6 n 3 + 567 R 8 n 4 5600 n 2 − 9 ( 315 − 1260 R 2 n + 1134 R 4 n 2 − 324 R 6 n 3 + 27 R 8 n 4 + . . . ) 1400 n 3 } .

(14)

3.4. Other approximations

Other approximations proposed in the literature are either the modified and enhanced version of the aforementioned approximations or formulated based on a combination of the previous approaches. For a model corresponding to a random walk on a diamond lattice, James and Guth [24] obtained an expression (unpublished but mentioned by Wang and Guth [9]) for the case of perfectly flexible chains as

P JG ( t ) = C JG [ 2 exp ( − 2 fA ) ( sinh 2 ( f ) + 5 8 + cosh ( f ) sinh ( f ) A ) ] 3 n 2 ,

(15)

where A = 3 t and f must be extracted from

A = sinh ( f ) ( sinh 2 ( f ) + 1 4 ) 1 / 2 .

(16)

Malacarne et al. [29] proposed an approximation based on the non-Gaussian Tsallis distribution which is an improved version of the Gaussian model. Their proposal is efficient and simple to implement for a large number of monomers. Accordingly, the normalized distribution takes the form as

P Mal . ( R ) = ( 1 π n 2 l 2 ) 3 / 2 Γ ( 3 n 2 ) Γ ( 3 n 2 − 3 2 ) ( 1 − R 2 n 2 l 2 ) 3 n / 2 − 5 / 2 ,

(17)

where Γ represents the gamma function.

Ilg et al. [13] proposed an estimation based on a one-dimensional version of a finitely extensible dumbbell model. The proposal utilizes trigonometric functions to approximate the inverse Langevin function in the formula by Kuhn and Grün [4] as

P Ilg ( R ) = C Ilg ( nl π R sin ( π R nl ) ) 2 ,

(18)

where C Ilg is the normalization factor. Recently, Khiêm and Itskov [30] revised the approximation by Ilg, equation (18), as

P KI ( R ) = C KI ( nl π R sin ( π R nl ) ) 9 π 2 ,

(19)

where C KI denotes the normalization factor.

Very recently, Morovati and Dargazany [31] proposed multiplicative corrections of the Kuhn–Grün approximation (8) as

P MD 2 ( t ) = P KG ( t ) [ exp ( − 1 . 75 t ) ( 1 − t ) 2 ] ,

(20)

P MD 3 ( t ) = P KG ( t ) [ exp ( − 2 t + 0 . 29 t 3 ) ( 1 − t ) 2 ] .

(21)

Note that all proposed approximation functions must satisfy the normalization condition (6).

4. Comparison and discussion of existing approximations

The above discussed approximations are compared in Figure 2 for the case of n = 10 . Both the normalized PDFs and the corresponding normalized entropic energies are shown along with their respective relative and absolute errors. Generally, relative errors tend to provide a better assessment of the approximative performance. However, this is not the case for the entire range of R / nl as the exact PDF and all its approximations converge to zero with R / nl approaching one. This leads to extremely small values and large relative errors. In cases like these, the relative error is not meaningful. Hence, it is reasonable to judge the approximative capabilities by utilizing the absolute and relative errors side-by-side.

OPEN IN VIEWER

Figure 2.

Comparison of the different normalized PDFs ( n = 10 ) and the corresponding normalized energies along with their respective relative and absolute errors.

Accordingly, the steepest descent, equation (12), and the three-term Morovati and Dargazany approximation, equation (21), significantly outperform the others for the example case of n = 10 . However, this conclusion cannot be generalized for arbitrary n without further investigation. All previously proposed approximations have been proven to perform well in single cases chosen by the respective authors leaving other important domains of n undiscussed. Thus, a deeper investigation of the approximations over a reasonable range for n is necessary to ensure proper comparability.

In the following, all approximations discussed in the previous sections are evaluated for every integer n within n ∈ [ 3 , 120 ] . For each case, the absolute and relative errors were computed. The mean and the maximum of these errors are plotted in Figure 3. Extreme extensibility R / nl → 1 is less important for practical applications as the chain force tends to infinity. Therefore, in this study, the errors are only evaluated in a range limited to R / nl < 0 . 95 . This limitation allows the comparison of different models in a more practically relevant and realistic way.

OPEN IN VIEWER

Figure 3.

Comparison of the mean and average errors of different normalized PDFs (top row) and normalized energies (bottom row) over a range of n . Note that some curves with large errors have been removed while cropping the subfigure to the area of interest.

According to the results plotted in Figure 3, the approximative qualities of the approaches vary significantly for different n . While there are some approximations that outperform the others in the majority of cases, there is no clear superior approximation in the whole domain and with respect to every considered error measure. In order to find the most suitable approximation, one hence has to base the choice on the specific value of n of the polymer chain studied, whether the PDF itself or the corresponding entropic energy should be approximated, and the error measure deemed most important for the considered application. Since the exact solution (7) explicitly depends on the choice of n , a suitable approximation can be expected to do as well. However, all previously discussed approximations lack this feature.

5. An improved approximation based on a Padé approximant

Many previously discussed approximations share a common structure containing of up to three parts. The core of these approximations is the Kuhn–Grün approximation P KG ( t ) corrected by a factor K cor ( t , n ) . Finally, the normalization condition (6) is enforced using a scaling factor C , which results in a generalized approximation formula as

P ( t , n ) = C P KG ( t ) K cor ( t , n ) .

(22)

In the present work, a correction function K cor ( t , n ) based on a Padé approximation is proposed. For a given function f ( t ) , the Padé [32–34] approximant of order [L/M] is a rational function of the form

[ L M ] f ( t ) = P ( t ) Q ( t ) = ∑ i = 0 L p i t i 1 + ∑ j = 1 M q j t j = p 0 + p 1 t + p 2 t 2 + ⋯ + p L t L 1 + q 1 t + q 2 t 2 + ⋯ + q M t M ,

(23)

where L and M denote the integers determining the order of the polynomials in the nominator and denominator, respectively. The choice of q 0 = 1 does not limit the generality but provides uniqueness of the Padé approximant for given L and M .

In the present work, we set L = 3 and M = 2 , and consider the maximal extensibility of a polymer chain at t = 1 (where P → 0 ). Accordingly, a PDF based on the Padé approximation can be given by

P ( t , n ) = C P KG ( t ) K cor ( t , n ) , K cor ( t , n ) = a ( n ) t 3 + b ( n ) t 2 + c ( n ) t + 1 ( 1 − t ) 2 ,

(24)

with three parameters a , b , and c each depending on n . C can be calculated in a straight-forward fashion using the normalization condition (6). For convenience, C can also be factored into K cor ( t , n ) , which leads to q 0 = C and linearly scales a , b , and c . Practically, this can be achieved by either optimizing for a , b , and c and then calculating equation (6) or directly applying the normalization as boundary condition while optimizing for a , b , c , and C . Both procedures will provide the same solution.

Usually, the parameters a , b , and c are obtained from partial sums of the Taylor series expansion around some point [35]. This is not feasible in the present context, where the function f ( t ) to approximate is not given in an explicit form and hence cannot be utilized. For this reason, the Levenberg–Marquardt method [36,37] implemented in MATLAB [38] was applied to obtain the parameters a , b , and c for every particular value of n .

The specific target function for the optimization procedure has a special significance discussed in the following. While a least-squares-based approach targeting the PDF itself may be an obvious choice, it is neither the only nor necessarily the best one. Depending on the desired application, minimizing the relative error or the mean/maximum of one of the error norms can be advantageous. For example, it appears reasonable to focus on the optimization of ln ( P ) rather than on P itself due to the higher significance of the normalized energy compared to the pure PDF in the context of constitutive modeling of polymers.

For the sake of this exemplary study, we minimized the least-squares of the absolute error of ln ( P ) . The so resulting values for a , b , and c are shown in Figure 4 as functions of n . This allows for perfect tailoring of the approximation to a specific application. In the relevant range of n , the parameters a , b , and c show a behavior, which can be easily fitted to simple functions, allowing quick and easy calculation (without any new optimization). The corresponding functions, as depicted in Figure 4, can be represented as

a ( n ) = − 0 . 42412 + 0 . 44 n ,

(25)

b ( n ) = 1 . 2675 + ln ( 25 . 8579 n ) 36 . 8362 ,

(26)

c ( n ) = − 1 . 8427 − n 4018 . 4835 .

(27)

OPEN IN VIEWER

Figure 4.

Parameters a , b , and c obtained for up to n = 120 . The optimization target was to minimize the least-squares of the absolute error of ln ( P ) . In addition, an approximation for each curve is given.

6. Results

In the following, the Padé approximation-based approach featuring the functional dependence on n , as proposed in the previous section, is compared to the steepest descent approach and both the two-term and the three-term Morovati and Dargazany approximation which revealed the best performance in comparison to other approximations in a wide range of n .

In Figure 5, the maximal absolute error of ln ( P ) and the maximal relative error of P are plotted versus n . One observes that the proposed approximation demonstrates the best performance in the whole range of n . In addition, the proposed approach leads to an error less dependent on n . While the advantage is especially prominent for the absolute error of ln ( P ) , which was set as the optimization target, it is still present in the PDF itself. Except for very short chains ( n < 7 ) , the relative error of the PDF is never higher than 1.3% and the absolute error of ln ( P ) never exceeds 0 . 01 .

OPEN IN VIEWER

Figure 5.

Comparison of the Padé-based approximation as given in equation (24) to the best approximations from the literature. The left picture shows the absolute error of ln ( P ) (which was set as optimization target), while the right shows the resulting relative error of P .

In a wide range of n , the proposed Padé approximation-based approach is capable to provide stronger predictions in comparison to the others even when used with constant a , b , and c . Figure 6 shows the resulting error, when utilized with averaged constant values for a , b , and c . While the dependency on n enhances the approximation for n < 20 and 100 < n , the effect is minimal for chains of medium size.

OPEN IN VIEWER

Figure 6.

Comparison of the proposed Padé-based approximation as given in equation (24) with averaged parameters ( a = − 0 . 4136 , b = 1 . 4640 , and c = − 1 . 8587 ) to the best approximations from the literature.

7. Conclusion

In the current work, we showed that most of the previously proposed approximations of the non-Gaussian probability distribution function (4) of the polymer chain length perform well only for a particular number of chain segments n . The effect of n on the approximation quality was investigated. Although some approximations outperform others, in the majority of cases, there is no clear superior approximation in the whole domain and for every considered error measure. The proper choice depends on the particular chain under investigation (determined by n ) and the desired application. As a result, there exists no versatile unique approximation to meet all the requirements.

In this paper, we derived a family of approximations of the non-Gaussian probability distribution function on the basis of the Padé method. The developed method is able to incorporate the polymer chain length and can be tailored to any specific application. For example, formulation (24)–(27) provides the minimal absolute error of ln ( P ) . In a wide range of n , the proposed approach appears to be the most accurate in comparison to other approximations known from the literature. A simplified approach, which neglects the dependency on the number of chain segments n , was shown to still outperform most of the alternative approximations.

The proposed approximation function is very robust for the modeling of polymer blends made up of short chains and also can minimize errors in simulating multimodal polymer networks which contain a wide range of polymer chain sizes [39]. Moreover, the non-Gaussian theory of polymer networks has been widely utilized to model the swelling behavior of cross-linked gels with finite extensibility of networks of high charge densities. The accuracy of such models can be considerably enhanced by the use of the proposed approach [40,41].