Behavior of viscoelastic models with thermal fluctuations

Fluctuating viscoelasticity for conformation-tensor-based models is studied at equilibrium, in simple-shear deformation, and in uniaxial extension. The models studied are the upper-convected Maxwell model, the FENE-P model with finite chain-extensibility, and the Giesekus model with anisotropic drag. Using numerical simulations, the models are compared in detail both with each other and with analytical predictions for the Maxwell model. At equilibrium, the models differ only marginally, both in terms of static and dynamic characteristics. When deformed, the average mechanical response of the Maxwell model is unaffected by the strength of thermal fluctuations, while the mechanical response of the FENE-P and Giesekus models show a slight decrease the stronger the fluctuations in simple shear, whereas the decrease in uniaxial extension is marginal. For all models, the standard deviation of the mechanical response increases with increasing strength of fluctuations, and the magnitude of the standard deviation relative to the average for given fluctuation strength generally decreases the stronger the deformation, this effect being stronger for uniaxial extension than for simple-shear deformation.


Introduction
Thermal fluctuations in viscoelastic fluids become important if the length scales of observation and possibly confinement are of the same order of magnitude as the characteristic length scale of the, often meso-scale, constituents of the fluid. This is because the smaller the number of constituents (e.g. polymer chains or segments, liquid crystal rods, colloidal particles) in the volume of observation, the more significant become fluctuations around the average collective behavior. One application where the involved length scales are small, and hence fluctuations can be relevant, is micro-and nanofluidics [1,2]. Another one is microrheology, where the motion of sub-micron sized tracer particles is studied with the use of microscopy to deduce the rheological properties of the suspending fluid [3][4][5][6].
Different routes can be taken for modeling fluids at small scales. Either one can use a particle-based approach, related to e.g. molecular dynamics, dissipative particle dynamics or Brownian dynamics, or one prefers to adopt a field-theoretic approach, related to fluid dynamics. The latter route is what is examined further in this paper.
Several procedures have been proposed for including thermal fluctuations in fluid dynamics. For Newtonian a e-mail: m.huetter@tue.nl b e-mail: m.a.carrozza@tue.nl c e-mail: m.a.hulsen@tue.nl d e-mail: p.d.anderson@tue.nl fluid dynamics, enrichment with thermal fluctuations has been achieved by Landau and Lifshitz [7]. For non-Newtonian fluids, it was proposed to relate the rate-of-deformation tensor to the stress tensor via a memory kernel and to add colored noise to the stress tensor [8,9]. Furthermore, a multi-scale model has been developed in [10,11], which adds elasticity and colored noise to the Newtonian fluid model. Alternatively, enriched smoothed-particle hydrodynamics has been used [12], which can be regarded as a discretized numerical approximation to a continuum model. Recently, Hütter et al. [13] developed a general approach for including thermal fluctuations in conformation-tensor-based viscoelastic models, in accordance with thermodynamic principles. The authors formulated these models in terms of a "square root" of the conformation tensor, namely what they call the contravariant deformation. In Carrozza et al. [14], it has already been shown that by using the contravariant deformation formulation in simulations of viscoelastic fluid flow without fluctuations, the numerical stability is enhanced compared to the conformation-tensor formulation. While the square root of the conformation tensor is not unique, which is discussed in detail elsewhere [15], the formulation in terms of the contravariant deformation appeals because it is closely related to the microstructure and therefore has a more useful physical interpretation than other square roots.
A systematic study of the intrinsic, generally nonlinear, behavior of complex fluids in the presence of thermal fluctuations is a natural next step towards a more generic approach and more flexibility in applications of fluctuating viscoelasticity. This should include variations in the strength of the fluctuations, as well as a comparison of different rheological models, for unraveling the characteristic behavior.
The goal of this paper is to examine in detail the behavior of three viscoelastic models with fluctuations, at equilibrium and in flow, by means of numerical simulation. The three models examined are the upper-convected Maxwell model, the FENE-P model, and the Giesekus model, the formulation of which in terms of the contravariant deformation with fluctuations has been established in Hütter et al. [13]. The reason behind choosing these particular models is the following. When formulated in terms of the conventional conformation tensor and in the absence of fluctuations, the Maxwell model is linear, while the other two models are non-linear, for distinct reasons. The FENE-P model owes its non-linearity to accounting for finite extensibility of the chain in the corresponding free energy. In contrast, the Giesekus model is non-linear because of anisotropic mobility. Since free energy and mobility are the key ingredients for formulating complex-fluid models along non-equilibrium thermodynamic principles, these three models are considered prototypical, and will thus be studied in this paper.
The following notation will be used in this paper. Summations are indicated by Σ (i.e. no Einstein summation convention is used), and the summation indices run over all spatial dimensions, unless indicated otherwise. Furthermore, (e x , e y , e z ) denotes the (right-handed) set of orthonormal basis vectors in Cartesian space.
The paper is organized as follows. The viscoelasticfluid models with thermal fluctuations, formulated in terms of the contravariant deformation, are introduced in sect. 2, and some predictions are formulated. In sect. 3, the viscoelastic models are examined numerically at equilibrium, in simple-shear flow, and in uniaxial extension. The results are discussed and conclusions are drawn in sect. 4.

General aspects
Consider a complex fluid, the microstructure of which shall be characterized by a symmetric and positive definite conformation tensor c, representative of, e.g., the conformation of polymer coils. Given a model of viscoelasticity formulated in terms of this conformation tensor, the contravariant deformation b introduced by Hütter et al. [13] is related to c by way of where the superscript "T" denotes the regular operator (matrix) transpose. It has been shown that it is advantageous to use b as a fundamental variable when incorporating thermal fluctuations in a viscoelastic model, see [13] for details. A physical interpretation of this finding can be given as follows. As stated in sect. 1, fluctuations originate from the number of constituents of the fluid in the volume of observation being relatively small. In the case of polymer chains, a detailed description would make use of vectors for describing the structure of the chains, e.g. in terms of N end-to-end vectors of chains or of chain segments, with N being relatively small. There are thermal fluctuations on these vectors, because the entities making up the polymers (e.g. monomers, coarse-grained beads) are agitated by the surrounding molecules. Therefore, fluctuations are naturally represented on the level of vectors.
Since the conformation tensor c can be interpreted as the second moment of the distribution of the end-to-end or segment vector [16,17], the decomposition (1) suggests that b is more closely related to the vector interpretation than c is. And therefore, fluctuations are more readily incorporated on the level of b, rather than on the level of c. Furthermore, the kinematics of b, i.e. its behavior in flow, mimics that of a contravariant vector. More precisely, the column vectors of b are chosen to obey contravariant behavior [13]. In the following, the upper-convected Maxwell model, the FENE-P model, and the Giesekus model are presented in compact form; for more details, the reader is referred to the original development [13]. To include thermal fluctuations properly, the general equation for the non-equilibrium reversible-irreversible coupling (GENERIC) [17][18][19] has been employed [13], whereby it is guaranteed that the fluctuation-dissipation theorem is respected.
All three models have a characteristic relaxation time λ and a shear modulus G. Throughout this paper, timerelated quantities are scaled with the relaxation time λ, while stress-related quantities are scaled with the shear modulus G, in order to make these quantities dimensionless. For example, t will be used to denote the dimensionless time, L stands for the dimensionless (transpose of the) velocity gradient, and σ is the dimensionless stress tensor. The dynamics of b is presented in the form of stochastic differential equations (SDE). For all SDEs reported in this paper, it is understood that the Itô interpretation of stochastic calculus [20,21] is used.

Thermodynamics for finite N
The essential building blocks in the thermodynamic approach taken in [13] are the free energy per unit volume ψ and the relaxation tensor Λ (4) . In [22], these two quantities have been derived by statistical mechanics, departing from the underlying vector description, for a finite number N of chains or chain segments, respectively. It turned out that there is a finite-N correction Δψ to the thermodynamic limit ψ ∞ of the free energy density [22], with where G = nk B T for entropy elasticity has been used with number density n, and D stands for the number of spatial dimensions. While the free energy used in the applications in [13] thus needs to be refined, according to eq. (2) and eq. (3), the relaxation tensor requires no correction for finite N . As a result of including Δψ in the thermodynamic derivation [13] of the three models presented below, the corresponding drift terms are modified. It is noted that the corresponding models presented below are fully compatible with the c-dynamics derived directly on the basis of the corresponding vector-descriptions, see [22] for details.
This finite-N correction is a manifestation of the smallness of the volume element considered. In the following, models will be presented that are based on such finite-N extension of macroscopic "bulk" expressions for the free energy. It should be mentioned, however, that by doing so the presence of walls -which alter the polymer conformations (e.g., see [23])-is not taken into account. Therefore, the models as presented below are applicable to small volume elements in the context of micro-/nanofluidics and microrheology where small volume elements are required, as long as the volume element is not immediately adjacent to the confining surface; there, strictly speaking, another free energy function would have to be used.

Upper-convected Maxwell model
The b-representation of the fluctuating dynamics of the upper-convected Maxwell model, in the following simply referred to as the Maxwell model, and the corresponding stress tensor are given by (see [13], and using eq. (2) and eq. (3)) with ν = 1 − (D + 1)Θ. The four contributions on the right-hand side (r.h.s.) of eq. (4) represent, in this order, deformation, relaxation, thermal drift (closely related to the fluctuations), and fluctuations. The dimensionless parameter Θ quantifies the importance of the thermal fluctuations with respect to the characteristic elastic energy, and is defined by Θ = k B T /(GV), with k B the Boltzmann constant, T the absolute temperature, and V the size of the volume of observation. The smaller the volume of observation, the more important are the thermal fluctuations. For entropy elasticity, one finds Θ = 1/N . In eq. (4), the symbol dW denotes the dimensionless increments of a multicomponent Wiener process, representing white noise. This second-order tensor is specified by the following averages and covariances: where t and t denote two moments in time, 1 (4) is the fourth-order unit tensor with components [1 (4) ] ijkl = δ ik δ jl , and where the covariance of two quantities A and B is denoted by A; B ≡ AB − A B . Equations (6), (7) show that dW is uncorrelated in time and that its components are independent of each other.

FENE-P model
The dynamics and the stress tensor are given by (see [13], and using eq. (2) and eq. (3)) These equations differ from the Maxwell model in sect. 2.3 only by the dimensionless factor f = β/(β + 3 − tr(c)), which describes the finite extensibility of the chain, with β a constant. This model reduces to the Maxwell model for β → ∞.

Giesekus model
The dynamics and the stress tensor are given by (see [13], and using eq. (2) and eq. (3)) with dW 1 and dW 2 increments of two statistically independent multicomponent Wiener processes, which individually obey the statistics described by eqs. (6), (7). The dimensionless parameter α determines the magnitude of the anisotropic drag, and D is the dimensionality of the model. The Giesekus model reduces to the Maxwell model for α → 0.

Analytical calculations and predictions
In this section, analytical calculations and predictions are made, in order to be able to rationalize in more detail the results of the numerical simulations presented further below. The two non-linear models, FENE-P and Giesekus, are not amendable to a detailed analytical analysis, particularly when fluctuations are included. In contrast, the Maxwell model is more straightforward to analyze, which is presented in the following. To that end, it is useful to present, in addition to the b-formulation (4), also the c-formulation of the Maxwell model (see [13], and using eq. (2) and eq. (3)), where D is the dimensionality of the model. The dynamics of b, (4), has the following properties. In the absence of flow (L = 0) and fluctuations (Θ = 0), one finds from the dynamics (4) that the contravariant deformation must obey b = b −1,T at equilibrium, while from the dynamics of c (12) one finds c = 1. The fact that other states than b = 1 are admissible is a consequence of the decomposition (1) not being unique. To any b that satisfies eq. (1) for given c, an orthogonal transformation can be multiplied to the right of b without affecting c (see also [15] for a more detailed discussion). In the presence of small fluctuations, 0 < Θ 1, we can expect c to be close to 1 at all times, while b, however, may keep rotating randomly in equilibrium by virtue of the fluctuations. Transferring the considerations about isotropically distributed random rotations described in appendix A to the evaluation of b · b T eq = c eq , one expects for the average and variance of the components of b, where c eq = 1 has been used (see eq. (17) below). A further prediction that one can make concerns the correlation time for fluctuations in b. Particularly, it can be shown (see appendix B) that there are two (dimensionless) correlation times, namely, which are related to stretch (s) and rotation (r), respectively, with ϑ < D. In contrast to appendix B, ϑ in eq. (16) effectively contains an average over many states around which the fluctuations are examined. Equation (15) shows that stretch decorrelates with the characteristic relaxation time of the polymer, while the decorrelation of rotation is delayed by a factor 2/((D − ϑ)Θ), see eq. (16). The properties of the dynamics of c, (12), can be assessed as follows. Taking the average of eq. (12), and realizing that the fluctuation contribution on the r.h.s. vanishes because of the Itô interpretation (meaning that b and dW are uncorrelated), one obtains for the evolution of the average Therefore, the solution for the average in the presence of fluctuations is equal to the solution to the deterministic model, i.e. to eq. (12) with Θ = 0. Particularly, at equilibrium c eq = 1. If one considers start-up of simple shear, with Weissenberg number Wi = λγ for constant shear-ratė γ for t > 0, the conformation tensor is given by In contrast, if one examines start-up of uniaxial extension, with Weissenberg number Wi = λε for constant extensionrateε for t > 0, the conformation tensor is given by with w 1 = 1 − 2Wi and w 2 = 1 + Wi, for −1 < Wi < 1/2. Based on eq. (12), also an evolution equation for the covariance tensor c; c can be derived, using Itô calculus [20,21], which leads to in component notation. For the stationary states of simple shear (18) and uniaxial extension (20) specified above, and hence also for equilibrium (Wi = 0), these equations can be solved in closed form. All non-zero components of the covariance tensor are listed in table 1, for D = 3. In particular, it is noted that the covariance between any two different components of the conformation tensor vanishes at equilibrium (Wi = 0) and in uniaxial extension, but not in simple shear. In view of the expression for the stress (5), the statements made above about the conformation tensor transfer readily to statements about the stress tensor. Particularly, we note in passing that the variance of the off-diagonal components of the stress tensor at equilibrium are given by σ i =j ; σ i =j eq = Θ, which closely resembles fluctuation formulas for the stiffness matrix of solids [24], when written in dimensional form.
Finally, the following can be anticipated about the decorrelation of components of the conformation tensor. The discussion of b has brought forward two time scales, related to stretch and rotation, respectively. Since the conformation tensor, in contrast, has a well-defined equilibrium state and fluctuations around it will thus not suffer from continued rotations, it is anticipated that only the stretch relaxation-time, λ s , will be observed for the conformation tensor. Table 1. Components of the covariance tensor c; c for the Maxwell model in the stationary states of simple shear and uniaxial extension, respectively. Other components can be constructed by using the symmetries c ji = cij and c kl ; cij = c ij ; c kl , and are zero otherwise. The variances at equilibrium can be obtained from either of the two modes of deformation by setting Wi = 0.

Variance
Simple shear Uniaxial extension

General aspects
In this section, the behavior of the models presented in sect. 2 is examined numerically for D = 3. For the FENE-P model, the parameter β describing the finite extensibility of the chain is chosen to be β = 50 [25,26]. For the Giesekus model, the parameter α representing the anisotropy in mobility is set to α = 0.1 [27]. Forward-Euler time-discretization with step size Δt is applied to the equations in sects. 2.3-2.5 to obtain the solution for the contravariant deformation tensor b at equilibrium as well as in simple shear and uniaxial extension, by imposing a constant velocity gradient L, see eq. (18) and eq. (20). To start the time stepping, the initial condition b(t = 0) = 1 is used in all calculations, which is the square root of the average conformation tensor at equilibrium for the Maxwell model, see eq. (17). Random numbers to generate the Wiener-process increments dW are drawn from a standard normal distribution with a variance equal to Δt (see eqs. (6), (7)) using the Ziggurat method [28]. For this method, the random numbers have a period of 2 32 − 1.

Equilibrium
We start with analyzing the time-correlation of the fluctuations in b and c. The time-correlation function of a time-dependent quantity A is given by where τ is the dimensionless time-difference, or lag time. At equilibrium, obviously, the explicit dependence on the instance of time t disappears, and only the time difference τ is relevant. This correlation function will be studied in the following for components of b and c.
In view of the discussion in the previous section, one expects that two relaxation processes are present in the dynamics of b, with dimensionless relaxation times equal to unity (see eq. (15)) and proportional to 1/Θ (see eq. (16)), representative of stretch and rotations, respectively. This has been taken into account in setting up the simulations. Particularly, two types of simulations have been performed. First, the dynamics of b has been simulated for t sim = 800/Θ, which is preceded by an equilibration of duration t equil = 15/Θ that is not included in the calculation of the correlation function. Second, the dynamics of b has also been examined when after every time step a re-initialization (RI) of the spurious rotation is performed by setting b equal to the symmetric squareroot of c, which should get rid of the rotational relaxation altogether. For this modified b-dynamics as well as for the dynamics of c, t sim = 800, since the strength of the noise Θ should not affect the correlations. For these cases, no equilibration has been used for the Maxwell model, since the initial condition represents its average solution. In contrast, for the FENE-P and Giesekus models, t equil = 5 has been used. The time step was set to Δt = 10 −2 , which is small enough in view of the anticipated relaxation times. The corresponding number N s of samples, i.e. statistically independent trajectories, of the ensemble has been chosen carefully as to not interfere with the period of the random number generator; the detailed values are discussed in appendix C, in table 3.
As an example, the correlation function for b xx for the FENE-P model without RI is shown in fig. 1. In the semilogarithmic representation in the inset, the presence of two relaxation processes is clearly visible. In order to quan- tify the corresponding relaxation times, representative of stretch (s) and rotation (r), the following function is fit to the numerical results, Notably, for the time-correlation functions of components of b with RI and for components of c, there is no rotational relaxation (not shown in figure), and therefore only one exponential is used to represent the data (i.e., γ = 1). The fit is performed only in a certain range 0 ≤ τ ≤ τ fit , which is chosen in such a way that both relaxation processes give significant contributions. Particularly, for the b-dynamics without RI, τ fit = 1/(2Θ) at Θ = 0.03, and τ fit = 1/Θ at Θ = 0.1. For the b-dynamics with RI and for the cdynamics, τ fit = 1, independent of Θ. The relaxation times obtained by performing a leastsquares fit of eq. (24) to the time-correlation function C A obtained by numerical simulation are listed in table 2. All models show a relaxation processes that has a relaxation time of order unity, related to stretch (see eq. (15) Several observations can be made in fig. 2. First, for all three models, both diagonal and off-diagonal components of b fluctuate around zero, the standard deviation being slightly larger than 0.5. Particularly, for the Maxwell model, the numerically determined standard deviations differ from the analytical predictions based on eq. (14) by less than 1%. And second, as expected, the off-diagonal components of the conformation tensor c for all three models fluctuate around zero, while the diagonal components fluctuate around unity. There are only marginal differences between the models. For the Maxwell model, the averages and standard deviations determined numerically deviate from the analytical predictions by less than 1%.
Based on the results in table 2 and fig. 2, one observes that the differences between the models are quite small, which suggests the following. Apparently, the region in phase space probed by the thermal fluctuations is relatively small, more specifically, too small to feel a significant effect of the finite extensibility and anisotropic drag, which are higher-order effects. The observation that c is not affected by Θ can be rationalized as follows. For the Maxwell model, this fact has been proven in sect. 2.6. For the FENE-P model, which differs from the Maxwell model only in terms of the dimensionless factor f in the restoring force (see sect. 2.4), one observes that for β = 50 this factor f is approximated well by f 1 + (tr(c) − 3)/β, which is close to unity at equilibrium for | tr(c)| β, i.e., the FENE-P model is represented well by the Maxwell model under these circumstances. The Giesekus model differs from the Maxwell model only in terms of the mobility tensor, but not in terms of the Helmholtz free energy. Therefore, the Giesekus model has the same probability distribution at equilibrium as the Maxwell model (see also [13,22]), and therefore the average of the conformation tensor must be the same for these two models.

Simple shear and uniaxial extension
In this section, the effect of fluctuations on the rheological behavior is studied, for simple shear (18) and uniaxial extension (20). For most of simulations under deformation, the time step used is Δt = 10 −2 . However, in some cases Δt = 10 −3 has been used, namely i) for the FENE-P model in simple shear at Wi = 1 at Θ = 0.1, and at Wi = 3 and Wi = 10, as well as in uniaxial extension at Wi = 3, and ii) for the Giesekus model in uniaxial extension at Wi = 3 and Θ = 0.1. When using these time steps, a tenfold refinement of the time step changed the results for the averages and standard deviations by less than 1%. The number of samples are N s = 1.9×10 5 for the Maxwell and FENE-P models, while N s = 9.5 × 10 4 has been used for the Giesekus model, the latter model requiring twice as many random numbers as the former two models.
An example of start-up flow is shown in fig. 3. In order to compare the three models and rates of deformation and to highlight the effect of fluctuations more easily, the mechanical response is normalized by its value in the absence of fluctuations; those latter values can be found in appendix D for completeness.
Based on fig. 4 and fig. 5, the following observations can be made about the relative influence of the fluctuations. For the Maxwell model, increasing the fluctuationstrength Θ does not affect the average response, as expected based on eq. (19), independent of the Weissenberg number Wi. In contrast, the average response of the FENE-P and Giesekus models does change (decreases) upon increasing Θ, particularly so for the higher values of Wi, i.e. when the non-linearity in these models becomes significant. Both of these models show a comparable decrease of the average shear stress and average first normal-stress difference upon increasing Θ (approx. 4.5% at Θ = 10 −1 and Wi = 10), with the exception that the decrease of the shear stress for the FENE-P model is weaker (approx. 1% at Θ = 10 −1 and Wi = 10). For all models, the standard deviation increases with increasing strength of fluctuations Θ at all Wi, obviously, and in some cases the standard deviation becomes even larger than the average itself, see the square symbols in fig. 4 and  fig. 5. And, the relative magnitude of the standard deviation for given Θ generally decreases the higher the Weissenberg number Wi. It is noteworthy that the standard deviation relative to the average decreases more strongly when increasing Wi for the Giesekus model than it does for the other models. Finally, it is mentioned that, for the Maxwell model, the difference between the numerical results presented in fig. 4 and fig. 5 and the analytical predictions for the average (19) and standard deviation, based on table 1, is on the order of 1% or less.
The first normal-stress difference in uniaxial extension is shown in fig. 6. It is noted that the Maxwell model is examined only for relatively small Weissenberg numbers Wi, since no stationary solution exists for Wi ≥ 1/2, see eq. (21). In contrast to the case of simple-shear deformation, increasing the fluctuation-strength Θ in uniaxial extension does not leave the average response of only the Maxwell model unaffected, as expected based on eq. (21), but has also only marginal effect on the averages of the other two models. The only exception is that for both the FENE-P model and the Giesekus model there is a weak decrease of the average for Θ = 0.1 at Wi = 0.3 by approx. 2%. Similar to simple shear, the standard deviation in uniaxial extension increases with increasing strength of fluctuations Θ at all Wi, obviously, and also here the standard deviation becomes even larger than the average for some cases. The relative magnitude of the standard deviation for given Θ decreases the higher the Weissenberg number Wi, this effect being stronger for uniaxial extension than for simple-shear deformation. For the Maxwell model, the difference between the numerical results presented in fig. 6 and the analytical predictions for the average (21) and the standard deviation, based on table 1, is on the order of 1% or less.
It has been shown analytically in sect. 2.6 that the average response of the Maxwell model is not affected by the strength of the thermal fluctuations Θ. Formulating clearcut analytical predictions for the FENE-P and Giesekus models is cumbersome for two reasons. First, taking the average of the corresponding SDE for c (see [13] for details) and considering stationary state results in an equation that is not closed in terms of c . And second, one does not have knowledge about the probability distribution of c either. Both of these issues occur because there is non-linearity and multiplicative noise in the corresponding SDE.
We now return to the cases for which the standard deviation is larger than the magnitude of the average, see the square symbols in fig. 4, fig. 5, and fig. 6. In order to identify the border-line for this to occur, we seek pairs of values (Wi, Θ) that satisfy σ xy 2 = σ xy ; σ xy , (only for simple shear), One can study these conditions for a given model by numerical means, in general, i.e., by running a large set of simulations for various combinations of Wi and Θ, and subsequently extract the combinations (Wi, Θ) that satisfy the above conditions. For the Maxwell model, a more direct route can be taken, since analytical expressions for the averages (see eq. (19) and eq. (21)) and standard deviations (by way of table 1) are available. Using these relations for D = 3, the combinations (Wi, Θ) satisfying the conditions (25) and (26) can be determined readily; the solutions are shown in fig. 7 for simple-shear deformation, and in fig. 8 for uniaxial extension.
In both figures, only values Θ ≤ 1 are shown, but on physical grounds one can argue that reasonable choices should even be significantly smaller than unity. Particularly for entropy elasticity, the strength of the fluctuations, Θ, is equal to the inverse of the number of chains (or chain segments) in the volume of observation [13]. Therefore, it seems that for Θ O(10 −1 ) the modeling approach taken in this paper is not appropriate, and rather a particlebased approach should be taken instead. The curves shown in fig. 7 and fig. 8 are monotonously increasing for Wi > 0, and monotonously decreasing for Wi < 0.
The standard deviation is smaller than the magnitude of the average in the regions in the lower left and lower right corners in fig. 7 and fig. 8, whereas the standard deviation is larger than the magnitude of the average in the upper center. Naturally, the separation lines between these regions depend on the criterion looked at, and therefore there are separation lines associated with the shear stress and with the first normal-stress difference, respectively, in fig. 7. It is noted that, for simple shear ( fig. 7), the first normal-stress difference sets more stringent conditions than the shear stress for having the standard deviation smaller than the magnitude of the average. The graphs also contain the cases studied in fig. 4, fig. 5, and  fig. 6, for the Maxwell model. Particularly, points which are above the respective separation lines in fig. 7 and fig. 8 have indeed been identified already in fig. 4, fig. 5, and fig. 6 as the ones with standard deviation larger than the average.

Discussion and conclusions
The effect of thermal fluctuations on viscoelasticity has been examined, by analytical and numerical analyses. The models examined are based on a conformation tensor in principle, but earlier work has shown that it is beneficial to formulate the fluctuating extensions of such models in terms a "square root" of the conformation tensor, namely the so-called contravariant deformation [13]. The models examined are the most standard and simple model of all (Maxwell), a model that accounts for non-linearity in terms of the thermodynamics (FENE-P, with finite chainextensibility), and a model accounting for anisotropic mobility (Giesekus). From a non-equilibrium thermodynamics viewpoint, they are prototypical for a wide class of conformation-tensor-based models.
While analytical predictions can be formulated for the Maxwell model due to its simplicity, the other two models could only be studied numerically. The simulation results can be summarized as follows. All three models display two relaxation processes related to stretch and rotation, respectively. As deformation is imposed, the average mechanical response of the Maxwell model is unaffected by the strength of thermal fluctuations Θ. In contrast, the average response of the FENE-P and Giesekus models decreases slightly in simple shear upon increasing Θ, particularly so for the higher values of the Weissenberg number Wi; this decrease is marginal in uniaxial extension. For all models, the standard deviation increases with increasing strength of fluctuations Θ for given Wi. And, the magnitude of the standard deviation relative to the average for given Θ generally decreases the higher the value of Wi, this effect being stronger for uniaxial extension than for simple-shear deformation.
Fluctuations in the rheological response are ubiquitous in molecular-dynamics simulations of polymer solutions and melts at finite temperature, see e.g. [29][30][31][32]. For example, it has been observed that, in the steady state, the magnitude of fluctuations in the shear viscosity and in the first normal-stress coefficient decreases upon increasing the shear rate in the shear-thinning regime [29,30,32]. Representing the results of this paper in unscaled form, the same qualitative behavior is found for both the FENE-P model and the Giesekus model. A more detailed comparison with molecular-dynamics simulations in a future study will be valuable in order to highlight the capabilities and shortcomings of the coarse-grained approach in this paper.
Our results show that there are conditions under which, in stationary flow, the standard deviation of the mechanical response is larger than the magnitude of the average. This occurs particularly for low Weissenberg number Wi and relatively high strength of the fluctuations Θ. These results have been obtained from averaging over a large number of trajectories. However, assuming ergodicity, which seems justified for the models at hand, the same results could have been obtained by time-averaging over a single sufficiently long trajectory in stationary state. For the example of the shear stress (but completely analogously for the first normal-stress difference), this implies that when imposing simple-shear deformation with Wi > 0, the shear stress can become negative temporarily, which can be interpreted as the viscosity being temporarily negative, while the average shear stress and average viscosity are positive. Shear banding is often associated with a non-monotonic relation between the rate of deformation and the stress [33]. The relation between these findings and ours needs further studies; in our case, it is just the fluctuations around the average that are unusual, and also we do not observe a negative slope in the stress vs. strain-rate relation, but rather the ratio between stress and strain-rate becomes negative temporarily. Furthermore, in the literature, molecular simulations of polymeric glasses revealed the existence of negative moduli, however, in that case this occurred in small regions that remained rather stable in the course of time [34]. The consequences of our findings about the large standard deviation need to be investigated in further studies. at finite temperature, also levels of the free energy higher than the minimum are explored. From the shape (asymmetry) of the free-energy function in eq. (B.3), one can infer that the average b is stretched with respect to its size at the minimum of the free energy. This implies that the time scale for the decorrelation of the stretch is equal to the characteristic time scale of the polymer (the latter being λ in a dimensional representation). In contrast, the decorrelation of the rotation is delayed by a factor 2/((D − ϑ)Θ). It is to be expected that, for ϑ to represent a typical state at a given strength of the thermal fluctuations Θ, the value of ϑ decreases slightly upon increasing Θ, since in this case higher free-energy values are probed and the asymmetry of the free-energy function has a stronger effect. Admittedly, the analysis above is a simplification of the true dynamics in several respects. On the one hand, we can anticipate that b will show arbitrarily large on-going rotations in the course of the equilibrium dynamics, and therefore considering fluctuations about b = b 0 1 only and looking at fluctuations that are small are simplifications. On the other hand, the above result would need averaging over the various states populated at equilibrium. Despite these shortcomings, the above result gives an indication of there being two relaxation processes, related to stretch and rotation, respectively. 10 -1 10 0 10 1 Fig. 9. Stationary-state shear stress in simple-shear flow at Θ = 0, σ xy,Θ=0, for all three models. The Weissenberg numbers Wi are specified underneath each measurement.

Appendix C. Number of samples N s for equilibrium calculations
The random number generator is called whenever the random increments in the dynamics are determined, i.e., in every time step, over the entire duration of the simulation, and for all samples. As mentioned in the main text, the number of samples N s has been chosen such that, for the given time step and simulation time, the period of the random number generator is not surpassed. Table 3 lists the values of N s used for the simulations on which the timecorrelation functions are calculated. Several observations can be made about table 3. First, the simulations of the Giesekus model have about half the number of samples compared to the other two models, because the Giesekus model requires twice as many random numbers per time step. Second, for a given model, the b-dynamics without RI has about a factor Θ fewer samples compared to the other dynamics, because its respective simulation time for each sample is longer by a factor 1/Θ. And third, for the same reason, only for the b-dynamics without RI does the number of samples depend on Θ, whereas the number of samples for the other dynamics does not depend on Θ.
Keeping in mind that the equilibrium averages and variances for b have been calculated without RI in the dynamics, the same observations made about table 3 also apply to the number of samples reported in table 4.  the values at Θ = 0 for all three models and all deformations studied in the main text are shown correspondingly in figs. 9-11 in this appendix.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.