Fluctuating viscoelasticity based on a finite number of dumbbells

Two alternative routes are taken to derive, on the basis of the dynamics of a finite number of dumbbells, viscoelasticity in terms of a conformation tensor with fluctuations. The first route is a direct approach using stochastic calculus only, and it serves as a benchmark for the second route, which is guided by thermodynamic principles. In the latter, the Helmholtz free energy and a generalized relaxation tensor play a key role. It is shown that the results of the two routes agree only if a finite-size contribution to the Helmholtz free energy of the conformation tensor is taken into account. Using statistical mechanics, this finite-size contribution is derived explicitly in this paper for a large class of models; this contribution is non-zero whenever the number of dumbbells in the volume of observation is finite. It is noted that the generalized relaxation tensor for the conformation tensor does not need any finite-size correction.


Introduction
Fluctuations are particularly important when studying small systems. This also holds for fluids, including complex fluids, e.g., macromolecular and polymeric liquids. Small scales are involved, e.g., in microrheology [1] and micro-and nanofluidic devices [2,3]. For Newtonian fluids, i.e., fluids with a deformation-independent viscosity and a lack of memory, the dynamics on small scales could be described in terms of the fluctuating Newtonian fluid dynamics developed by Landau and Lifshitz [4]. However, this is not sufficient for complex fluids, and thus extensions are needed. For example, the stress tensor has been related to the rate-of-strain tensor by a memory kernel, and correspondingly colored noise has been introduced on the stress tensor [5,6]. Another approach towards modeling fluctuating effects in complex fluids has been taken by Vázquez-Quesada, Ellero, and Español [7] and applied to microrhe-ology [8], in which smoothed-particle hydrodynamics is extended by a conformation tensor that describes the conformation of the small number of polymer chains per volume element. The concept of fluctuating dynamics for the conformation tensor has been extended recently [9], to make it applicable not only to the Maxwell model [10,11], as in [7,8], but to a wider class of models, e.g. the FENE-P model [11,12] and the Giesekus model [11,[13][14][15]. In the approach taken in [9], the Helmholtz free energy in terms of the conformation tensor plays an essential role.
The dynamics of the conformation tensor roots in a finer description, in particular, it can be related to the kinetic theory of dumbbells (e.g., see chapter 13 in [16]). The question addressed in this paper is what lessons can be learned from deriving the dynamics for the conformation tensor with fluctuations from an underlying kinetic description for a finite number of dumbbells. It is pointed out that the dumbbell description already contains the relaxation and fluctuation effects that are relevant also on the conformation-tensor level. This is in contrast to coarse graining from an atomistic description to bead-spring chains or directly to the conformation tensor, e.g., see the work of Underhill and Doyle [17] and of Ilg et al. [18,19], respectively, without fluctuations on the conformation-tensor level.
The paper is organized as follows. In sect. 2, a certain class of kinetic dumbbell models is introduced, based on which a description of fluctuating viscoelasticity in terms of the conformation tensor is derived via a direct route, for a finite number of dumbbells. This route is paralleled in sect. 3, where a thermodynamic approach is taken to arrive at fluctuating viscoelasticity. As part of that, the finite-size correction to the Helmholtz free energy is calculated, and this is found to be essential for finding agreement between the two approaches. In the appendix, three different calculation methods for this free energy are detailed, each of which arrives at the same result. In sect. 4, the relation between the dumbbell models and the multiplicative decomposition of the conformation tensor, which has been discussed recently in the literature [9,20,21], is examined. The paper ends with conclusions and a discussion in sect. 5.
Throughout this paper, the following notation will be used. All summations are spelled out, i.e., no Einstein summation-convention is used for repeated indices. While the symbol ⋅ denotes a contraction of one pair of indices, we use ⊙ for a double contraction: For an order-four tensor A (4) and an order-two tensor B, A (4) The Kronecker delta is given as δ ij , and the Dirac deltafunction as δ(x − y). The dyadic product of two vectors v 1 and v 2 is written as v 1 v 2 .
It can be shown that the probability distribution for the system (1) at equilibrium is given by the Boltzmann distribution with energy Φ; to ensure that, the third term in the square brackets in eq. (1) is essential -if this term was absent, the equilibrium distribution would depend on the mobility tensor, which is unphysical. Note that the mobility M may depend on {Q µ } µ=1,...,N , however, it is assumed in this study that it is the same mobility tensor for all dumbbells µ, and thus it has no subscript µ.
For completeness, it is mentioned that all position vectors are dimensionless, i.e., they are scaled with respect to (w.r.t.) a characteristic constant length-scale, that is omitted throughout this paper for convenience.
In the following, expressions for Φ and M are considered that will allow us to eventually derive a closed evolution equation for the instantaneous conformation tensor Symmetry requires that the potential energy Φ d of a single dumbbell depends only on the (squared) length of the dumbbell vector. If the dumbbells are not interacting, the total potential energy of all dumbbells is given by the sum of the individual dumbbell contributions, i.e., Φ exact = N Φ d (tr (QQ)). However, in this paper, we consider models for which the total potential energy is obtained by interchanging the operations (. . .) and Φ d (tr (. . .)) in Φ exact , i.e., we use the relation where eq. (5) was employed. This mean-field ansatz will be beneficial for deriving closed dynamics forĉ. Obviously, eq. (6) is -for all but one case (see below) -only an approximation to the exact potential for all dumbbells. However, for the purpose of this paper (which is to examine the effect of finite N on the counting of states, and ramifications thereof for the dynamics ofĉ), we employ eq. (6) for defining the class of models we examine. Therefore, we will not go into details about how accurately eq. (6) approximates the exact potential. The energy Φ given by eq. (6) equals the exact energy Φ exact only if the function Φ d (tr (. . .)) is linear, in particular for linear-elastic Hookean springs, Φ d (tr (QQ)) = (H 2)tr (QQ), with H the spring constant of the dumbbell. However, other cases are of interest as well. The dumbbell force generally has the form −Hf Q µ . For example, for a linear-elastic Hookean spring, f = 1, while for a finitely extensible nonlinear elastic spring one can use 1 (see also [12]), with two constants β 1 and β 2 , and which accounts for the finite extensibility in a mean-field sense. Typically, in the FENE-P approximation [12] (see also [29]),ĉ in eq. (7) is used for infinitely many dumbbells. However, since in this paper the focus is on studying systems with a finite number N of dumbbells, we generalize this by using eq. (7) for finite N . The potential that corresponds to the desired expression for the force is given by of which the limit β 1 = 1 and β 2 → ∞ results in the potential for the linear-elastic Hookean spring. Therefore, using the approximation Φ, given by eq. (8) with eq. (7), instead of the exact energy Φ exact is safe for small deformations. For other deformations, it is noted that the potential Φ does include finite extensibility, albeit in a different way than if applied to each dumbbell individually. It is a topic of future research to examine the foundations of the meanfield approximation Φ for finite N thoroughly. In this paper, this approximate expression for the energy forms part of defining the kinetic models that are examined in the following, and it is suitable for deriving closed dynamics for c, in the same spirit as the FENE-P approximation has been introduced earlier for infinite N [12]. As far as the mobility tensor M is concerned, its dependence on the dumbbell vectors {Q µ } µ=1,...,N is restricted to a dependence on the instantaneous conformation tensor, A particular realization of that is with friction coefficient ζ, and the parameter α is used to adjust the amount ofĉ-dependence. In case of imposed deformation, a non-zero value for α thus results in anisotropy of the friction tensor, where the anisotropy is introduced in a mean-field sense. The form (10) corresponds to the widely used Giesekus model for anisotropic drag [13][14][15]. For the Giesekus model, one typically uses eq. (10) for infinitely many chains, i.e., N → ∞, to render the model solvable. However, the mobility of a dumbbell is affected primarily by the other dumbbells in its vicinity, and therefore the finite-N generalization (10) is reasonable. Allowing for the dumbbell potential energy Φ d (tr (. . .)) in eq. (6) to be nonlinear and/or the dumbbell mobility tensor M to depend on the conformation tensor effectively introduces mean-field type couplings of the individual dumbbells. In practice this implies that either (a) the N dumbbells must be in the vicinity of each other so they can interact, or (b) they diffuse rapidly enough in space to effect such interactions. In the case of the mobility tensor, the implied physics is reasonably clear: the assumption is that the average orientation of surrounding dumbbells affects the mobility of any given test dumbbell. Having this rationale in mind suggests that this mean-field mobility Table 1. Overview of the parameters in the potential Φ given by eq. (7) and eq. (8), and the mobility tensor M , eq. (10), for the three models Dumbbell Parameters: makes more sense for a finite number N of dumbbells, than it does for infinitely many. In the case of the potential energy, the microscopic physics of the implied coupling between dumbbells is less clear. Still, we find it worthy of note that a potential energy function exists from which the FENE-P model can be derived exactly.
In order to highlight the overall structure of the modeling in the remainder of this paper, the general forms eq. (6) for the potential Φ and eq. (9) for the mobility tensor M will be used. These general results can then be reduced to the Hookean dumbbell model, the FENE-P model, and the Giesekus model, respectively, by appropriate choices for the forms and parameters for the potential Φ and the mobility tensor M , see table 1.

Transition from dumbbells to the conformation tensor
Given the definition of the instantaneous conformation tensorĉ, eq. (5), and using the Itô interpretation of stochastic calculus [22,30], one has in general where ⟨. . .⟩ Itô dt implies that in dQ µ only terms involving the Wiener increments are kept and subsequently reduced according to the rule (see Table 3.1 in [22]) Applied to the above class of models, the SDEs (1) for the dumbbell vectors {Q µ } µ=1,...,N can be transformed into an SDE for the conformation tensor, (13) in terms of the dumbbell potential energy Φ d and the dumbbell mobility tensor M , where it has been assumed that M ⋅ĉ =ĉ ⋅ M . It is pointed out that the symmetry ofĉ must be taken into account explicitly when calculating the partial derivatives of M (see [9] for details). The quantity dĉ f denotes the thermal fluctuations, (14) with B given by eq. (2) for a general mobility tensor M , eq. (9). It can be shown that the fluctuations have the properties The SDE (13) for the conformation tensor with fluctuations obeying the statistical properties given by eq. (15) and eq. (16) is the benchmark to which the thermodynamic treatment further below will be compared, for the general class of models described by energy Φ, eq. (6), and mobility tensor M , eq. (9). If the potential Φ and the mobility tensor M are of forms more general than eq. (6) and eq. (9), respectively, the dynamics for the conformation tensor would not close automatically, in which case one would have to employ procedures of coarse graining [24,[31][32][33]. The procedure presented in this sec. 2.2 has also been followed in [7] for deriving theĉ-dynamics for the Hookean and FENE-P models.
It is noted that, in the limit of many dumbbells, i.e., N → ∞, not only do the fluctuations dĉ f become insignificant. In addition, also the "thermal drift", i.e., the secondlast contribution on the right-hand side (r.h.s.) of eq. (13) vanishes as well, and thus the conventional deterministic dynamics for the conformation tensor is recovered. The Hookean dumbbell model -which results in the Maxwell model for the conformation tensor -, the FENE-P model, and the Giesekus model are sub-cases of the SDE (13) with noise (14) when choosing the parameters β 1 , β 2 , and α appropriately, see table 1.

Fluctuating viscoelasticity derived using thermodynamics
The main idea in taking a thermodynamic approach to modeling dynamical systems is that one can concentrate on key ingredients for the static and dynamic properties, and the thermodynamic approach makes sure that these ingredients are processed towards the final model in a consistent way. For the dumbbell dynamics (1), the key ingredients are the potential Φ and the mobility tensor M . In contrast, for the dynamics of c, the key ingredients are the Helmholtz free energy density ψ and the generalized relaxation tensor Λ (4) (see [9] for details). In particular, the dynamics for the fluctuating conformation-tensor c is given by the SDE where div c Λ (4) denotes the divergence of Λ (4) in c-space, the order-four tensor B (4) satisfies and dW is a tensor with increments of independent Wiener processes (see [9] for further details). The structure of the SDE (17) with eq. (18) is completely analogous to the one for the dumbbell models, eq. (1) with eq. (2). In the following, we derive expressions for ψ and Λ (4) , based on those for Φ and M , respectively.

Free energy density ψ(c) for finite N
The Helmholtz free energy Ψ = Ψ (c) for the symmetric conformation-tensor c for a finite number N of dumbbells is given by Ψ = −k B T ln Z, with the canonical partitionfunction where D is the number of spatial dimensions. The Kdimensional Dirac δ-function makes sure that only those states in {Q µ }-space are accounted for that are compatible with the conformation tensor c. Sinceĉ is symmetric by definition, see eq. (5), only K = D(D + 1) 2 independent conditions are needed (instead of D 2 ); no more conditions are required for properly restricting the integration in {Q µ }-space.
It is pointed out that δ (K) is actually a δ-function in c-space in the sense that ∫ δ (K) (ĉ − c)d K c = 1. Using this latter relation, the integral of Z(c) over all (symmetric) conformation tensors c reduces to which is the conventional canonical partition function in the absence of the constraintĉ = c.
As we restrict our attention to energy functions Φ which depend on {Q µ } µ=1,...,N only by way ofĉ, see eq. (6), the canonical partition-function Z is related to the microcanonical partition-function Γ by way of with Different procedures for calculating the dependence of Γ on c explicitly are discussed in Appendix B, one based on deriving a differential equation for Γ , another one with a more geometrical interpretation, and a third one using a scaling argument. Following any of these procedures, the result for finite N is with a c-independent prefactor Γ 0 . Based on eq. (21) with eq. (23), the Helmholtz free energy density ψ = Ψ V per volume V becomes with and number density n = N V , and where we have omitted a c-independent additive constant, which is irrelevant for the formulation of the dynamics of c according to eq. (17).
Since Φ is proportional to N , it is evident that the first two contributions on the r.h.s. of eq. (24) are independent of the size of the system, for given number density n. In contrast, the third contribution, ∆ψ, does depend on the size of the system, in particular it becomes more relevant the smaller the system. To the best of our knowledge, this finite-size correction to the Helmholtz free energy (density) has not been derived earlier.
Using H = k B T and G = nk B T , the Helmholtz free energy density (24) with eq. (25) for the three models discussed in table 1 agrees with standard literature (e.g., see [10,11,16]), and with what has been used in the fluctuatingviscoelasticity approach in [9], with the important difference of the finite-size correction ∆ψ. Using, in contrast to our procedure, a continuous (representative of N → ∞) distribution for the dumbbell vector Q (e.g., see [16]), the thereby-derived Helmholtz free energy density corresponds to eq. (24) where the finite-size correction ∆ψ is absent.

Relaxation tensor Λ (4)
In the dumbbell dynamics (1), structural relaxation is expressed as −M ⋅ (∂Φ ∂Q µ ) dt. In the c µ -dynamics (17), structural relaxation is expressed as −Λ (4) ⊙ (∂ψ ∂c) dt, with an order-four relaxation tensor Λ (4) [9]. In the translation from M to Λ (4) , a reduction of variables and the volume of the system V are involved, the latter being necessary since Φ is an energy while ψ is an energy density. The relation between M and Λ (4) is given by (see also [31], and sect. 6.4 in [24]) where the contractions run over the components of Q µ and M , and ⟨. . .⟩ is the average over {Q µ }-space for given c.
Since M depends on the positions {Q µ } µ=1,...,N only by way ofĉ, taking the average is thus equivalent to replacinĝ c by c everywhere in eq. (27). It is to be noted that there is no finite-size correction in Λ (4) . Using again G = nk B T and with ζ = 4k B T λ, the relaxation tensor Λ (4) given by eq. (27) for the three models in table 1 turns out to agree with the standard expressions in the literature (e.g., see [10,11,16]) , and with what has been used in [9] in the context of fluctuating viscoelasticity.

Application to fluctuating viscoelasticity
According to the general procedure in [9], represented in eq. (17), and using the Helmholtz free energy density (24) with eq. (25) and the relaxation tensor (27) withĉ → c, one observes that the results in [9] need to be amended by including the finite-N contribution In particular, one obtains for the complete fluctuating dynamics where the symmetry of c has been taken into account when calculating div c Λ (4) (see [9] for details). See Appendix A for explicit exemplary applications of this equation. The symbol dc f denotes the fluctuating contribution given by where b satisfies the condition In general, B in eq. (30) has dimensions D ×P with P ≥ D (as described in sect. 2.1), b has dimensions D × P ′ with P ′ ≥ D, and therefore dW has dimensions P ′ × P ; for practical purposes, one may choose P = P ′ = D. The tensor dW consists of increments of statistically independent Wiener processes, with the properties When using component notation, eq. (33) implies that any two of the components of dW are independent from each other.
A direct comparison shows that the SDE (13) derived directly from the dumbbell model and the SDE (29) derived via the thermodynamic route, respectively, agree. It is noted that the expressions for the fluctuations, dĉ f in eq. (14) and dc f in eq. (30), respectively, have a different form. However, it can be shown that they have the same statistical properties. First, both representations are linear superpositions of increments of Wiener processes, and second, for both representations the average is given by eq. (15) and the covariance by (16). The difference in the expressions for the fluctuations is not a short-coming of the approach; it rather reflects the non-uniqueness of the decompositions (18) and (31). The non-uniqueness of the decomposition (31) can actually be utilized for relating the expressions (14) and (30) for the noise in even more explicit terms. Specifically, choosing b to be D × N with the column vectors of b equal to Q µ √ N (µ = 1, . . . , N ) (see sect. 4 for a further elaboration), and setting the row vectors of dW equal to dW µ (µ = 1, . . . , N ), one finds that the expressions (14) and (30) are identical.
In deriving the relaxation term in eq. (29), i.e. the first term on the r.h.s. that is proportional to M , from the general form eq. (17), one notices the following: The prefactor 1 N in Λ (4) given by eq. (27) is cancelled by the prefactor N in the derivative ∂(ψ − ∆ψ) ∂c, for ψ given by eq. (24) with Φ according to eq. (6). If one chose to not cancel these factors, one would observe that (1 N )M (with factor 1 N ) is the relevant mobility on the conformationtensor level not only for the thermal drift and the fluctuations (by way of the covariances (16)), but also for the relaxation.
The importance of the finite-size correction ∆ψ in the free energy density, eq. (25), for the evolution equation (29) is pointed out. In particular, dc ∆ψ exactly cancels those contributions from div c Λ (4) that are related to the derivative of the explicit factorsĉ in Λ (4) , eq. (27). If ∆ψ was neglected, agreement between the SDEs (13) and (29) could not be achieved.
In Appendix A, the dynamics for the conformation tensor with fluctuations, eq. (29) with eq. (30), is presented explicitly for three models, namely for the Hookean dumbbell / Maxwell model, the FENE-P model, and the Giesekus model.

Comments on the multiplicative decomposition of c 4.1 Eliminating degrees of freedom
Above, the relation has been established between the dynamics formulated in terms of dumbbell vectors, on the one hand, and in terms of the conformation tensor, on the other hand. In this section, the relation of the dumbbellvector description to the multiplicative decomposition of the conformation tensor (e.g., see [9,21]), is examined, whereb has dimensions D × P ′ with P ′ an arbitrary dimension. This decomposition can be written in the formĉ withb P ′ ,ν the ν-th column vector ofb P ′ (see also [20]). In general, one obviously must require P ′ ≥ D for this decomposition to be complete for arbitrary conformation tensorĉ.
In view of the expression (5) for the conformation tensor, a natural choice is P ′ = N withb N,µ = Q µ √ N , relating the dynamics ofb N directly to that of the dumbbell vectors Q µ . In the following, we focus on Hookean dumbbells, i.e., the Maxwell model, for illustrative purposes. The evolution equations (1) for Hookean dumbbells, using the potential (8) and mobility tensor (10) with the parameters given in table 1, become where the identifications λ = ζ (4H) and H = k B T have been made. Eq. (36) translates directly into the dynamics ofb N , where dW N has dimensions D × N with its µ-th column vector given by dW µ . Let us now compare this result with the dynamics for the "square root" b D of the conformation tensor c, i.e. c = b D ⋅ b T D where b D has dimensions D × D, as derived in [9] and amended in [21], The close relation between b D and the elastic, i.e., recoverable, part of the deformation gradient in solid mechanics has been discussed in [9]. For N → ∞, the dynamics of the column vectors of b 3 agree with the treatment proposed in [34]. Two major differences between eq. (37) and eq. (38) are apparent. First, the relaxation in eq. (37) drives the column vectors ofb N to zero, while , considering e.g. D = 3, according to eq. (38) 1 det b 3 , and b 3,3 to b 3,1 × b 3,2 det b 3 , respectively, i.e., the column vectors of b 3 become orthonormal in the course of relaxation. The second major difference between eq. (37) and eq. (38) relates to the absence of the thermal drift (the third term on the r.h.s. of eq. (38)) in eq. (37). Both differences, in relaxation and thermal drift, are a hallmark of eliminating degrees of freedom when going to a reduced description of the dynamics. More specifically, both of these contributions are tightly related to the entropy, i.e., to the counting of states in configuration space, see sect. 3.1. It is pointed out that the difference in relaxation does not depend on the value of N , while the thermal drift clearly is a finite-N effect, i.e., it is related to the fluctuations. However, despite these differences between eq. (37) and eq. (38), one should keep in mind that they both result in the same dynamics for the conformation tensor.

Rotational dynamics
In [21], it has been discussed that, due to the non-uniqueness of the decomposition c = b D ⋅ b T D , there is on-going rotational dynamics in eq. (38) even at equilibrium, with a relaxation time that depends approximately linearly on N . In the following, an attempt is made to rationalize this N -dependence of the rotational relaxation time in terms of the dynamics for the dumbbell vectors Q µ , eq. (36). Consider the linear combination Based on eq. (36), its dynamics is given by where the Wiener process increments defined by dW = ∑ N µ=1 γ µ dW µ satisfy with γ 2 ≡ ∑ N µ=1 γ 2 µ . Let us consider equilibrium, κ = 0. In order to study the orientation dynamics, we write R = Rn, with R and n the length and the orientation vector of R, respectively. By using Itô calculus [22,25], it can be shown that the SDEs for R and n are given by The contributions proportional to γ 2 originate from the second-order term in the Itô calculus 2 . In particular, one observes in the orientation dynamics (44) that the relaxation time for the orientation vector n is given by As an example, consider the case D = 3. To get an idea about the rotational dynamics of the column vectors b 3,µ of b 3 , which satisfy c = ∑ end, think of a Voronoi tessellation on the sphere, generated by the six "poles" that themselves are generated as intersections of three orthogonal axes with the sphere. Let us now consider two of these Voronoi sectors, V i and Vī, which are on opposite sides of the sphere. It is reasonable to assume that the number of vectors Q µ which have their orientation in these two Voronoi sectors together can be written as ϕN , where ϕ is independent of N . After mapping all vectors Q µ with orientation in Vī into V i by way of Q µ → −Q µ , which is just making use of the symmetry of dumbbell description, all ϕN vectors with orientation in V i are averaged to obtain R, i.e., γ µ = 1 (ϕN ) for these vectors and γ µ = 0 for all others. In this case, one finds γ 2 = 1 (ϕN ), and thus the relaxation time for rotation (45) is increased w.r.t. λ, namely as λ n ∝ N λ, in agreement with the observation in [21]. For the construction of a suitable vector R, the Voronoi sectors have been used, instead of, e.g., considering a random selection of vectors Q µ , for two reasons: First, the construction with Voronoi sectors allows to construct three such vectors (representative of the column vectors of b 3 ) that are clearly linearly independent from each other. And second, the length of R is proportional to the length of the dumbbell vectors Q µ with a prefactor that is of order O(1 2), i.e. independent of N .

Discussion and conclusions
The focus of this paper has been on deriving viscoelasticity in terms of a conformation tensor with fluctuations, based on the kinetic theory of dumbbells. This has been achieved by identifying the conformation tensor with the arithmetic average over a finite number N of dumbbells, eq. (5), and using two alternative routes for deriving the dynamics: a direct approach using stochastic calculus, and a thermodynamic approach, in which the Helmholtz free energy plays a key role. It has been shown that these two approaches agree only if a finite-size contribution to the Helmholtz free energy of the conformation tensor is taken into account. The main messages of this paper are therefore the following: -If the number N of dumbbells is finite, the commonly employed expressions for the thermodynamic potentials need to be corrected: Using statistical mechanics (see Appendix B), one finds that the conformational entropy must be corrected by replacing N by N −D −1 (with D the number of spatial dimensions), which in turn modifes the Helmholtz free energy, see eq. (24) with correction term eq. (25). -The thereby obtained finite-size correction in the free energy is crucial for guaranteeing compatibility between the dynamics of the conformation tensor and of the underlying dumbbells.
While these general conclusions have been established in general terms for a large class of models, they have also been exemplified for the Hookean (Maxwell) model, the FENE-P model, and the Giesekus model (see Appendix A); the dynamics for the conformation tensor with fluctuations for these three models is summarized in table 2. When discussing a model with fluctuations, the deterministic counterpart serves as a benchmark. For all models discussed in this paper, one recovers the known deterministic models in the thermodynamic limit N → ∞, i.e., if the number of dumbbells in the volume of interest V diverges, keeping the number density n = N V constant. Beyond that thermodynamic limit, however, there is also an interest in the behavior of the average conformation tensor ⟨ĉ⟩ for finite N , i.e., in the presence of fluctuations. For most models studied in this paper, the nonlinearities do not allow to obtain a closed form equation for ⟨ĉ⟩ based on the stochastic differential equation (SDE), eq. (13) and eq. (29), for the fluctuating conformation tensorĉ. The notable exception to this rule is the Hookean dumbbell (Maxwell) model, with potential energy Φ d (tr (QQ)) = (H 2)tr (QQ) and mobility tensor M = (2 ζ)1. Taking the average of the SDE (A.1) over different realizations of the fluctuations, one observes that the average conformation tensor ⟨ĉ⟩ obeys the same differential equation as its deterministic (N → ∞) counterpartan observation that generally does not hold for nonlinear models. In particular, one finds for the Hookean dumbbell model at equilibrium ⟨ĉ⟩ eq = 1, where H = k B T has been used. As a word of caution, it is pointed out that the conformation tensorĉ that minimizes the Helmholtz free energy density, eq. (24) with eq. (25), is given bŷ c min = (1 − (D + 1) N )1, i.e., it does depend on the finite size (N ) of the system. The fact thatĉ min ≠ ⟨ĉ⟩ eq is not a contradiction; it merely points out that the distribution of thermal fluctuations around the minimum is not symmetric.
The relevance of the finite-size correction of the Helmholtz free energy, ∆ψ given by eq. (25), has been discussed primarily in the context of formulating dynamics with fluctuations for the conformation tensor c. However, it is also of immediate consequence for the formulation of fluctuating viscoelasticity in terms of the "square root" b 3 , where c = b 3 ⋅ b T 3 . In [9], a thermodynamic approach has been taken towards deriving the dynamics of b 3 , based on the dynamics of c. Therefore, if the thermodynamic potential for the c-dynamics contains a finite-size contribution, the same holds true also for the thermodynamic potential for the b 3 -dynamics, see sect. 4.3 in [9] for details. Beyond these implications for the thermodynamics of a b 3 -formulation, the kinetic models for N dumbbells have also been employed in this paper to give an explanation for the existence of a rotational relaxation time proportional to N in the fluctuating dynamics of b 3 , which has been observed earlier [21].
MH acknowledges stimulating discussions with Hans Chris-tianÖttinger, particularly in relation to the interpretation of eq. (38). PDO is grateful to Georgetown University and the Ives Foundation for support. A preliminary version of the derivation in Appendix B.2, for the N → ∞ limit, was developed by DJR during previous discussions with Joseph Peterson and Gary Leal.

Author contribution statement
Peter Olmsted initiated this collaboration of the authors, and Daniel Read spotted the inconsistency between the results in [9] and a direct kinetic-theory approach. Markus Hütter has made the main contributions to the main part of the manuscript, as well as to Appendices A, B.1 and B.3.1, while Daniel Read and Peter Olmsted have developed Appendices B.2 and B.3.2, respectively. All authors have discussed in detail about all parts of the paper, and helped to bring the paper in its final form. All the authors have read and approved the final manuscript.

Conflict of interest
The authors declare that they have no conflict of interest.

A Dynamics of conformation tensor for three exemplary models
For completeness, three concrete realizations of the general dynamics for the conformation tensor, eq. (29) with eq. (30), are provided in this section.

A.1 Hookean, i.e., Maxwell, model
Using the potential Φ in eq. (8) with eq. (7) in the limit β 1 = 1 and β 2 → ∞, and with the mobility tensor M given by eq. (10) for α = 0 for Hookean dumbbells, the general equation (29) turns into the Maxwell model with fluctuations (see also [9,21]), where the fluctuations are determined by eq. (30) with eq. (2). In particular, one may choose both b and B to have dimensions D × D and to be symmetric, i.e., b = √ c and B = √ M = 2 ζ 1. This result for the stochastic dynamics of the conformation tensor is identical to what has been derived in [7].

A.2 FENE-P model
Using the potential Φ in eq. (8) with eq. (7), and with the mobility tensor M given by eq. (10) for α = 0, the general equation (29) turns into the FENE-P model with fluctuations (see also [9,21]), The fluctuations are determined by eq. (30) and eq. (2) where, as for the Maxwell model, one may choose both b and B to have dimensions D × D and to be symmetric, i.e., b = √ c and B = √ M = 2 ζ 1. Table 2. Dynamics for the conformation tensor c with fluctuations; see Appendix A for the derivation. Symbols are explained in the text. Note that the square root of a tensor is symmetric. The increments dW of Wiener processes satisfy ⟨dWij(t)⟩ = 0 and ⟨dWij(t) dW kl (t ′ )⟩ = δ ik δ jl δ(t − t ′ ) dt dt ′ (see eq. (32) and eq. (33)).

Model Conformation dynamics Noise
Maxwell

A.3 Giesekus model
Using the potential Φ in eq. (8) with eq. (7) in the limit β 1 = 1 and β 2 → ∞, and with the mobility tensor M given by eq. (10), the general equation (29) turns into the Giesekus model with fluctuations (see also [9,21]), where the fluctuations are determined by eq. (30) with eq. (2). One may choose both b and B to have dimensions D × D and to be symmetric, i.e., b = √ c and B = √ M , where M is given by eq. (10) for α ≠ 0. The therebyobtained expression for the fluctuations differs from that used in [9], while sharing the same statistical properties. The difference in the expressions for the fluctuations reflects the non-uniqueness of the decomposition (18).
Of the three models discussed explicitly in this Appendix A, the Giesekus model is the only one for which there is a thermal drift in the dynamics ofĉ, i.e., the second-last contribution on the r.h.s. of eq. (A.3). It is present only if the mobility of the dumbbells is anisotropic, α ≠ 0. Beyond the N -dependence of the noise dĉ f , the thermal drift is the second explicit consequence of the finite size of the system.

B Calculation of partition function Γ
In this Appendix, three procedures are presented for calculating the microcanonical partition function Γ , eq. (22), for a finite number N of dumbbells.

B.1 Procedure 1: Differential equation
We start by noting that the microcanonical partition function Γ , eq. (22), contains only K = D(D + 1) 2 (rather than D 2 ) Dirac δ-functions, representative of the constraints related to the independent components c ij with 1 ≤ i ≤ j ≤ D. The following notation is introduced: . With eq. (B.2), the microcanonical partition function can be written in the form One can show that where ∆ ′ ij denotes the derivative of ∆ ij w.r.t. its argument. Furthermore, one can derive the following relations, where Q µ,k is the k th component of the µ th dumbbell, of which there are three non-zero cases: In calculating eqs. (B.7)-(B.9), we have made use of the identity xδ ′ (x) = −δ(x), with x =ĉ ij − c ij , which implies To proceed, we use eq. (B.5) for i = j and eq. (B.7) to derive an expression for 2c ii (∂Γ ∂c ii ), where we have used ∆ . The remaining integral is then re-written by performing an integration by parts w.r.t. Q µ,i ; the corresponding boundary-terms can be neglected if the values of the components of c are finite. This leads to where, again, we have used ∆ . Using the product rule for calculating the derivative ∂∆ ∂Q µ,i and rearranging terms results in where the quantities G ji and G ⋆ ij are given by where we have used eq. (B.8) and eq. (B.9) (with i and j interchanged), and then eq. (B.4) and eq. (B.5). Combining eq. (B.12) with eq. (B.14) and eq. (B.13), one obtains which is a differential equation for Γ . To solve this equation, consider the ansatz For calculating the partial derivatives of this ansatz w.r.t. c ij , the following needs to be kept in mind. On the one hand, since c is symmetric, there are only K independent variables, rather than D 2 , which is in line with the strategy adopted, e.g., in eq. (B.2) and eq. (B.3). This implies that there are only K partial derivatives to be calculated in eq. (B.15), namely w.r.t. c ij with 1 ≤ i ≤ j ≤ D. On the other hand, the tensor c contains elements c ij and c ji = c ij . Hence, for i ≠ j, one finds for any function f (c): where "no-sym" emphasizes that the corresponding derivative is taken without enforcing the symmetry of c, i.e., considering all D 2 components of c as independent variables. For the derivative of ansatz (B.16) one thus obtains where the factor (2 − δ kl ) originates from the fact that there are only K independent components in c, c kl with k ≤ l; particularly, for k < l, the variable c kl appears at two off-diagonal positions in the full, symmetric matrix c. With this, it can be shown that the ansatz (B.16) is indeed a solution of the differential equation (B.15) if The c-independent prefactor Γ 0 can not be determined with this approach. However, since Γ 0 results only in an additive contribution to the Helmholtz free energy Ψ , it turns out to be irrelevant for the dynamics of the conformation tensor c.

B.2 Procedure 2: Geometry
We begin with a geometrical interpretation of the microcanonical partition-function Γ (c) defined via eqs. (21) and (22). Substituting from eq. (21) back into the full partition function eq. (20) we find: Comparing the two lines of this equation we recognise that we may interpret the differential quantity Γ (c) d K c as being the volume, within the DN dimensional {Q µ }-space, that is within an increment d K c of c (i.e. where the K independent components of c are each varied within an interval dc α of their base value). Evaluation of the dependence of this volume on c yields Γ (c) up to a constant prefactor. Since this derivation requires some visual imagination, we first demonstrate how this may be reasoned in the specific case D = 3 before indicating how the calculation may be generalised to arbitrary D.
We also note that, for D = 3, there are K = 6 independent components of the symmetric c tensor: c xx , c yy , c zz , c xy , c xz and c yz .
We proceed by evaluating the dependence on c of two separate volumes. We first evaluate V const , which is the (3N −6)-dimensional subvolume of {Q µ }-space within whicĥ c (given in terms of {Q µ } by eq. 5) is held exactly equal to c. Then, for each point within V const , we find the volume dV 6 which is the 6-dimensional subvolume swept out by varying c xx by dc xx , c yy by dc yy , etc. Such excursions must all be perpendicular to the subvolume V const (in the {Q µ }-space), because contours of fixed c xx are perpendicular to the gradient direction of c xx (etc.). So, we can evaluate Γ (c) d 6 c as the product of V const and dV 6 : Evaluation of these volumes is most straightforward in the co-ordinate frame in which c is diagonalised (which can always be done since c is symmetric). In the diagonal frame, c takes values λ x , λ y and λ z along the diagonal, but is zero in the off diagonal components. In terms of the N -dimensional vectors X, Y and Z these constraints can be expressed as: where Λ i = N λ i , i.e. the vectors X, Y and Z are restricted to the surfaces of N -dimensional hyperspheres of radii √ Λ x , Λ y , and √ Λ z respectively, whilst simultaneously being held to be mutually perpendicular.
The volume V const is the (3N − 6)-dimensional volume swept out by rotating the X, Y and Z vectors subject to the above constraints, the rotations being within the Ndimensional space of these vectors. A brief analogy may help at this stage: the surface of a sphere of radius r is found by summing up small tiles formed by varying the polar co-ordinate angles θ and φ by small increments dθ and dφ. The distance along the sphere surface moved during such increments is rdθ and rdφ (multiplied by a geometric factor sin θ which is irrelevant to the scaling with r). So, the scaling of surface area with r can be found from the product of these lengths, r 2 dθdφ. Likewise, the volume V const is the sum over a tiling of small incremental volumes, made by rotating the vectors X, Y and Z by small angles dθ k in each available direction, whilst keeping their length fixed and retaining their mutually perpendicular orientation. Since X has N dimensions, there are in total (N −1) directions in which X could be rotated whilst keeping the length of X fixed. One such rotation will rotate the X vector towards the Y vector. In this case, the Y vector must also rotate by the same angle, so as to maintain the perpendicular condition X ⋅ Y = 0. Hence, for rotation angle dθ 1 , X rotates so that its end sweeps out a length dl X = √ Λ x dθ 1 perpendicular to X (in the direction of Y). Likewise Y also rotates so that its end sweeps out a length dl Y = Λ y dθ 1 perpendicular to Y (in the direction of −X). These changes in X and Y result in a total length moved in {Q µ }-space which is Similarly a second rotation direction of the X vector is available, towards the Z vector (in which case Z must also rotate). Likewise, a third rotation carries Y towards Z without rotating X. By similar arguments, the total length moved in {Q µ }-space for these is: These three rotation directions having been dealt with, there remain (N − 3) further rotation directions available for the X vector, all of which allow X to remain perpendicular to both Y and Z. For each of these rotations of X, the length moved in {Q µ }-space is Similarly, there are (N −3) rotations available to each of Y and Z, each of which leaves the other vectors unchanged, giving lengths: A single "tile" in the volume V const is obtained by sweeping through each of the above (3N − 6) incremental rotations. The total volume V const is obtained by adding together all such tiles as the vectors are rotated. We require only the dependence of V const on c (ignoring prefactors). Hence, by taking the product of the lengths swept out by a set of incremental rotations, we obtain: We now turn to the volume dV 6 , the 6-dimensional subvolume swept out by incrementing c xx by dc xx , c yy by dc yy , etc. away from a single point within V const . We first calculate the length in {Q µ }-space traversed by incrementing c xx by dc xx in a direction perpendicular to the contour of constant c xx . We note that c xx = X ⋅ X N, so the direction perpendicular to the contour of constant c xx is ∇c xx = (2X N, 0, 0) .
Moving a distance dl xx in this direction changes c xx by Inverting this, the distance moved in {Q µ }-space when incrementing c xx by dc xx is and similarly for dc yy and dc zz . We next note that so the direction perpendicular to the contour of constant c xy is ∇c xy = (Y N, X N, 0) Moving a distance dl xy in this direction changes c xy by dc xy = ∇c xy dl xy Inverting this, the distance moved in {Q µ }-space when incrementing c xy by dc xy is dl xy = N Λ x + Λ y dc xy and similarly for dc xz and dc yz . Hence, multiplying these six lengths together, Evaluating Γ (c) d 6 c = V const dV 6 , and hence Γ (c), we find that all factors of form Λ x + Λ y cancel, and we are left with But detc = λ x λ y λ z , and so: with dependence on detc as required.

B.2.2 Generalisation to D dimensions
Most of the above formalism carries directly over to D dimensions. In calculating V const , each vector of type X has a total number of (N − 1) rotations available (whilst preserving length), but (D − 1) of these rotate the vector towards one of the others. So, the total number of rotations available for each vector, which do not also require rotating one of the other vectors, is (N − D). The number of rotations requiring two vectors to rotate (i.e. of the form X to Y, Y to Z etc., but avoiding double counting) is M = D(D − 1) 2. So, V const is of form: where there are M terms of form Λ x + Λ y . In calculating dV K , there are D diagonal constraints and M = D(D − 1) 2 off-diagonal constraints in c (a total of K = D(D + 1) 2). Each diagonal constraint gives a factor of form √ Λ x in the denominator of dV K , whilst each off diagonal constraint gives a factor of form Λ x + Λ y . Hence: In the product Γ (c) d K c = V const dV K , all M factors of form Λ x + Λ y are present once in the numerator, and once in the denominator, and so cancel. Hence:

B.3 Procedure 3: Scaling
In this procedure, scaling arguments are employed for the calculation of the microcanonical partition function Γ , eq. (22), in real space and in reciprocal space, respectively.

B.3.1 Scaling in real space
Similar to the procedure described in sec. B.2, it is again chosen to describe the {Q µ }-space in terms of the Ndimensional vectors X i , 1 ≤ i ≤ D, where the µ-th component of X i equals the i-th component of Q µ .
In what follows, we consider the coordinate system in which the conformation tensor c is diagonal, with eigenvalues λ i (1 ≤ i ≤ D). The microcanonical partition function Γ , eq.
where Γ 1 comes from the substitution of variables in the volume element, and Γ 2 comes from the substitution of variables in δ (K) . It is to be noted that the integral Γ 3 does not depend on c. Rather, we find Γ 3 = Γ (1). Therefore, one obtains for the c-dependence of the microcanonical partition function where there are L = D(D−1) 2 zero off-diagonal elements. Hence we introduce D Dirac δ-functions for the diagonal constraints, and L δ-functions for the off-diagonal constraints. We enforce these with Fourier transforms, leading to u a (λ a −ê a ⋅ĉ ⋅ê a ) − i where the second sum is over all L distinct pairs of eigenvectors (b < c). Inserting the definition forĉ (eq. (5)) and interchanging ∑ µ and ∑ a in the exponential, Γ (c) can be written as where the matrices U and V are constructed from the constraint fields u a and v α as follows (represented in the basis {ê a }) for D = 3 and we have rescaled v by a factor of 2. Next, we scale the u and v integrals according to We finally rescale Q µ =Q µ ⋅ Λ, where Λ is the diagonal matrix Λ aa = √ λ a , Λ ab = 0(a ≠ b). Hence, we find Γ (c) = 2 L (2π) K (det c) where the integral I is independent of the conformation tensor c, andŪ +V is equal to (B.40) without the factors 1 √ λ a λ b .