In the present work, the macroscopic governing equations governing the heat and mass transfer for a general multicomponent system are derived via a systematic nonequilibrium thermodynamics framework. In contrast to previous approaches, the relative (with respect to the mass average velocity) component mass fluxes (relative species momenta) and the heat flux are treated explicitly, in complete analogy with the momentum flux. The framework followed here, in addition to allowing for the description of relaxation phenomena in heat and mass transfer, establishes to the fullest the analogy between all transport processes, momentum, heat, and mass transfer, toward which R. B. Bird contributed so much with his work. The inclusion of heat flux-based momentum as an additional variable allows for the description of relaxation phenomena in heat transfer as well as of mixed (Soret and Dufour) effects, coupling heat and mass transfer. The resulting models are Galilean invariant, thereby resolving a conundrum in the field, and always respect the second law of thermodynamics, for appropriate selection of transport parameters. The general flux-based dynamic equations reduce to the traditional transport equations in the limit when mass species and heat relaxation effects are negligible and are fully consistent with the equations established from the application of kinetic theory in the limit of dilute gases. As an added benefit, for the particular example case of hyperbolic diffusion we illustrate the application of the proposed models as a method to allow the use of powerful numerical solvers normally not available for solving mass transfer models more generally.
Considerable effort has been devoted already in the macroscopic modeling of mass and heat transfer within multicomponent mixtures. Regarding multicomponent diffusion, one can trace its origin to a generalization of the empirically derived Fick's law for a binary system.1–3 However, for a more formal description, one had to wait the development of nonequilibrium thermodynamics (NET).2 This also allowed one to take advantage of the Onsager–Casimir's reciprocal relations4–6 that resulted in significantly reducing the number of independent parameters, n(n-1)/2 instead of (n-1)2,2 where n is the number of diffusing species. In parallel with Fick's work, Maxwell developed an inverse relation correlating the gradients of concentrations to diffusive fluxes7,8 that was later extended to multicomponent systems by Stefan.9,10 This fundamental relation was eventually extended to include pressure and thermal effects by Curtiss and Hirschfelder11 and was finalized shortly afterward in the monumental treatise of “molecular theory of gases and liquids” by Hirschfelder, Curtiss, and Bird.12 More recently, the two approaches (the generalized Fick and the generalized Maxwell–Stefan) were shown to be related in the work of Curtiss and Bird13—in a recent article by Bird and Klingenberg14 these and other historic instances related to multicomponent diffusion have been nicely outlined and summarized. Of interest also is the very recent article by Swaney and Bird,15 where through the use of thermodynamic relations the authors develop equations for energy, mechanical energy, and entropy for multicomponent mixtures.
In parallel to the above works, several efforts have been dedicated to obtaining more general expressions for multicomponent diffusion applicable to arbitrary fluid systems. From the early works, it is worth noting the work of Truesdell16 in developing a more general form of Maxwell–Stefan theory using general principles of the mechanics of arbitrary continua. However, these efforts in modeling multicomponent diffusion within the at the time available continuum mechanics formulation led to indeterminate results when, by theoretical necessity, the modeling of the rate changes of the relative species momenta was also included.17 This unsatisfactory situation was specifically mentioned by Truesdell and Noll in Sec. 130 of their fundamental work, the “non-linear field theories of mechanics” and that statement remained present even in the last edition.17 In retrospect, as part of the outcomes from the present investigation, we attribute that problem in the introduction of species “stresses” when the balance of species momenta was introduced in that theory.17 We feel that species-specific stresses cannot be physically defined from a continuum point of view, as the individual species are not isolated in individual parts of space but are sharing the same space with all others. As a result, any physically imposed constraints following that definition are invalid. This problem issue is common to all such continuum theories that make use of similar assumptions, usually known as standard “mixture” theory. Later work by Lam18 trying to resolve the indeterminacy issue at the continuum level by resorting to collision theory-based explanations led to formulas too complicated to use (on top of the need to use several questionable assumptions). Other works, when attempting to generalize the multicomponent diffusion model, they did that at the constitutive level of fluxes, generalizing the Stefan–Maxwell relations through the use of nonequilibrium thermodynamics at that level—see, for example, Refs. 19–21.
The standard mixture theory and/or the use of species stresses has continued to be used in many recent works that have attempted to extract better expressions for the generalized Stefan–Maxwell relations.22–27 The difference between those works and the more traditional ones mentioned at the end of the previous paragraph that sought to express the diffusion fluxes statically in terms of constitutive relations directly relating them to potentials (such as gradients of the concentration and temperature or the chemical potential) is that in those an additional attempt is made to model their dynamics, following a multi-fluid approach. Such an approach has been followed in several other works in the literature, (primarily the two-fluid case; see, for example, Refs. 28–32) albeit in those works inertia effects and flux dynamics were ultimately neglected. Instead, the emphasis in those works was to capture the effect of nonequilibrium stresses to the flow developing in one of the two-fluid phases, typically corresponding to a polymeric system28–31 but also for suspensions of solid particles.32 Of special importance is previous work33–35 where a non-equilibrium thermodynamics-based, two-fluid theory was presented to describe the structure and rheology of shear banding rodlike micellar solutions. Even if inertia in diffusion fluxes may be often negligible (and a recent work, based on a statistical analysis near equilibrium, claims that this may be what happens for collision-dominated dilute ideal gases36) developing a consistent theory that can account for it may be very worthwhile. First, it is needed to predict finite speed of propagation of diffusive processes within matter, a phenomenon for which there is considerable experimental evidence,37 thus removing the aphysical infinite speed characterizing classical descriptions of heat and mass transfer. Second, it establishes needed consistency to the continuum modeling of transport processes and third it may be the best way to bring in the connection to nonequilibrium thermodynamics—it may also have significant applications under conditions far away from equilibrium as well as it may lead to novel numerical approximations even when flux inertia effects can be safely neglected (an example of which is illustrated further in this work).
In the present work, we elaborate on the modeling of flux dynamics for simultaneous heat and mass transfer in multicomponent mixtures with the objectives (a) correcting the mixture theory by consistently including inertial (kinetic energy) effects to diffusion and heat fluxes that may be important under conditions where an appreciable differential motion between various populations of the chemical species is present, (b) developing the theory for diffusion and heat transfer at the same level (using suitably defined momentum densities) as for momentum transfer. This is worthwhile as it offers a much needed closure to the full equivalence between all transport phenomena to which R. B. Bird has so much contributed, and (c) demonstrating a more straightforward evaluation of the extraction of generalized Stefan–Maxwell relations, when inertia flux effects can be neglected, as their connection to nonequilibrium thermodynamics becomes more straightforward at the momentum densities level. The addition of a momentum density associated with the heat flux allows one to consistently describe not only relaxation effects in heat transfer38–40 but also the mass and heat transfer coupling (the Soret and Dufour effects).
To achieve these objectives, we employ the single generator bracket formulation of non-equilibrium thermodynamics (SGBF-NET).41 This represents one of the simplest formulations of the most modern implementations of nonequilibrium thermodynamics, the GENERIC formulation42 to which it has been shown to be completely equivalent when applied to macroscopic modeling.43 The applicability of SGBF-NET and GENERIC to model transport problems has been shown many times, both in a variety of standard transport cases41,42 but also in many special cases, most recently that corresponding to developing thermodynamically admissible 13-moment equations, as described in another article belonging to the present special Physics of Fluids issue.44 Moreover, it is a thermodynamically consistent formulation that as we show can reproduce the exact functional forms for the heat and mass fluxes obtained from the kinetic theory of dilute gases.12 In this way, for this particular example, not only one can have full expressions (that are non-trivial!) for the dependence of the phenomenological transport coefficients on the system variables, but also one can hope, through the use of minimal adjustable parameters that can in principle be determined experimentally, to be able to extend these developments to more complex cases (e.g., dense gases and/or liquids) where any a priori kinetic theory-based prediction is lacking. However, we also need to note here that one may have a variety of ways through which heat transfer is carried over, depending on the material of interest and the conditions. The form used in the present work is suitable to describe heat transfer as established through collisions and based on the kinetic energy of the participating molecules, exactly as assumed in the kinetic theory of gases.12 Other forms of heat transfer may be possible, necessitating, in general, separate internal variables and fluxes—see, for example, Refs. 38–40 and 45.
In Sec. II, we show how the SGBF-NET can be used to naturally extend the modeling of heat and mass transfer within multicomponent systems to also include the effect of possible heat and different species momenta. In Sec. III, we show how this treatment is compatible in the limit of small differential velocities, to classical descriptions of heat and mass transfer, including the Fourier heat conduction, the Fickian diffusion in a binary mixture, and the Maxwell–Stefan equations, also corrected for coupling with heat transfer (Soret/Dufour effects). In fact, related to the latter case, we also show in Sec. III how the present development leads to equations that are fully compatible to those obtained from kinetic theory for dilute gases12 but also allowing to consistently extent those equations to other systems. This analogy allows us in the dilute gases case to obtain specific values for all the transport coefficients entering the theory. This is demonstrated, for a specific binary gas case, also shown in Sec. III. One further specific example is examined in detail in Sec. IV: That of the binary diffusion of a dilute solute in a liquid solvent. We show there, beyond the compatibility of the flux model to the standard Fickian diffusion in the limit of small differential velocities, the usefulness of the flux approach in setting up very simple but also very powerful numerical approximations, exemplified in the simple problem of transient diffusion in a slab. More specifically, taking advantage of the hyperbolic character of the flux-based model equations, we show how the use of explicit time integration solvers can reproduce results that otherwise would have required significantly more computational effort, demanding implicit numerical time integrations. Finally, in Sec. V, we present our conclusions.
II. SGBF-NET THEORY AND DERIVATION OF MULTICOMPONENT TRANSPORT EQUATIONS
We apply here the single generator bracket formalism of nonequilibrium thermodynamics (SGBF-NET) as developed in Ref. 41. As shown in Ref. 43, this is equivalent for macroscopic systems to the more general GENERIC approach.42 For completion, we start with a very short presentation of the bracket formalism of nonequilibrium Thermodynamics, borrowing from the description offered in a recent publication46 that outlines its influence in elastic fluid modeling—for a more detailed presentation the interested reader is referred to the above-mentioned research monograph.41
The bracket formalism offers a generalization of the Poisson bracket description (which is the most general description of Newtonian mechanics, see, for example, Refs. 47 and 48) from conservative, i.e., Hamiltonian, to dissipative systems. Of critical importance was the development first of Poisson brackets for conservative continuum systems, such as the Euler equations, by Glebsch49 and Arnold.50 The development of the dissipative bracket followed next.51–53 The bracket formalism also extends equilibrium to non-equilibrium thermodynamics, with key assumption that there is a set of variables (superset of the equilibrium ones) in terms of which a local equilibrium can be defined, one of the variables being a local entropy density. Within the bracket formalism, there are four steps that lead to the development of the system of evolution equations governing any given system:41,46
First, one needs to define a set of acceptable state variables, , in terms of which the system is specified. Those are typically field densities, including, at a minimum, those needed to specify the system at equilibrium.
Second, the Hamiltonian needs to be specified, physically identified with an (extended) energy of the system, as a functional (integral over the system volume) of the state variables, , the dependence on the first, second, and higher gradient being optional.
Third, the Poisson bracket (governing the reversible component of the system dynamics), , needs to be specified. It has all the properties of an acceptable Poisson bracket used to specify conservative Hamiltonian dynamics, i.e., it has to be an antisymmetric bilinear functional of the Volterra derivatives for any pair of arbitrary functionals, F, G, that also satisfies the Jacobi identity.41
Fourth, the dissipation bracket (governing the irreversible component of the system dynamics), , needs to be specified. This last component is introduced to account for dissipative phenomena; it is linear in terms of the first functional, but in general nonlinear in terms of the second. It also needs to satisfy (upon linearization, close to equilibrium) certain positive definite symmetric or antisymmetric properties between the various Volterra derivatives of the functionals of F, G depending on the symmetry upon time inversion of the variables with respect to which the Volterra derivatives are taken.41 These properties guarantee satisfaction of the Onsager/Casimir reciprocal relations of linear irreversible thermodynamics4–6,41 as well as the satisfaction of the second law of thermodynamics.
Once those four steps are carried out, the dynamics of any arbitrary functional F are simply specified as
If so wished, the evolution equations for the Eulerian partial time derivatives of anyone of the state variables, , can then easily be deduced by comparing Eq. (1) against the expression obtained through the application of differentiation by parts
and equating the factors in front of each Volterra derivative .41 The key advantage of the bracket approach lies in (a) being able to separate cleanly reversible and irreversible (dissipative) components in the dynamics by simply tracing their provenance to the Poisson and dissipative brackets, respectively, and (b) automatic satisfaction at the macroscopic level of all the irreversible thermodynamic requirements.
In the following, we will apply the bracket formalism to the flux-based modeling of heat and mass transfer. The first step involves the selection of the appropriate state variables. To describe heat and mass transfer between n species, , where , the standard set of variables are: (a) each one of the n species mass densities, , and (b) their corresponding momentum densities
where is the jth species mass-averaged velocity, (c) the entropy density, s. Note that the total mass density, ρ, and its corresponding total momentum density
(where the second equality is used for the definition of mass averaged velocity, ) are then obtained through a straightforward summation of the corresponding species quantities
For reasons that will become obvious in Sec. III, we propose to consider here the following equivalent set of primitive variables: (a) the species mass densities [through which the total mass density ρ, can also be obtained, from the first relationship in Eq. (5)] (b) the total momentum density, , (c), the reduced momentum densities corresponding to each species
(d) the total entropy density, . Note that in this case, the total momentum density () becomes an independent primary variable of the system, reestablishing the status it enjoyed for single component system. However, in exchange, the reduced momentum densities are not all independent. Instead, they must satisfy the following constraint:
In addition, note that, to ensure symmetry and consistency between the treatment of heat and mass transfer, we introduce one more variable, the reduced heat momentum density . This is defined in terms of the heat flux as
where the effective heat mass density, , is a material property with units of mass density entering in the definition of the contribution to the extended free energy density due to the heat flux (see below), and is a corresponding to the heat flux velocity with the physical meaning that it would have resulted to the same overall increase in the entropy if it was added to the fluid velocity. It therefore plays the same role for heat transfer as the relative species velocity does with respect to the diffusive flux of species i. After the addition of this new reduced heat momentum density, the appropriate state variable space can be defined as , with the understanding that not all the reduced species momenta are independent, as they need to obey the homogeneous constraint described by Eq. (7), thus effectively reducing the dimensions of the state space to .
The second step in the bracket formulation involved the description of the total energy (Hamiltonian) of the system. This is provided here as
where the last term is written in terms of the prescribed system variables with representing the equilibrium contribution to the energy density. Note that the heat mass density is not an independent state variable, but rather, an additional material property, . At this point, we cannot be more specific as to what this effective mass density is; from a units point of view, one could alternatively express it in terms of the entropy density, s, the temperature, T, and a characteristic speed, (in an equation, that is, introduced for convenience and shall be used as a proper definition for the characteristic speed, ) as
The additional term correcting the total Helmholtz free energy (at the lowest order and taken here as isotropic) for the presence of nonzero heat fluxes appears in Eq. (9) similar to the term introduced by Coleman et al.54 with the coefficient representing the magnitude of an isotropic tensor in the corresponding general tensorial and anisotropic expression proposed by Coleman et al. [Eq. (1.47) in Ref. 54].
Based on the general expression for the Hamiltonian provided by Eq. (9), one can readily calculate the corresponding Volterra derivatives of the Hamiltonian as
where, as usual, the jth chemical potential and temperature are defined as and j takes the values . Note that the Volterra derivatives with respect to the reduced species momentum densities had to be corrected in order to satisfy the same constraint as the reduced momentum densities, Eq. (7) (which can be expressed as the n-vector generated by the reduced momentum densities to be normal to the (n elements) row vector). Note that there is one-to-one correspondence between the set of variables used here and that involving the species momenta, with their relationships as well as that between the corresponding Volterra derivatives of an arbitrary functional G, summarized for convenience in Appendix A. The important feature of the variables chosen as system variables above is that now the Volterra derivatives with respect to the species densities have their unique contributions arising from the species kinetic energies that are Galilean invariant, i.e., independent of the velocity of the coordinate system in terms of which we represent our system.
Now, if desired, we can also involve additional state variables, for example, as in our previous work,33–35 we can associate with each component k, also a conformation tensor as the mass-weighted average of the second moment of a structural position vector, , corresponding to component k (assumed to be a macromolecule). Note that in such a case, there is an extra contribution to the Hamiltonian through an “elastic” free energy density, that also contributes to the Volterra derivative an extra elastic energy related term . This effect can lead to very interesting phenomena, such as stress induced migration,28 and as such will be studied in a latter work; however, for simplicity it is not considered in the following work presented in this manuscript.
Third, the Poisson bracket, which describes the only reversible process (convection) taking place, can only involve one and only one genuine velocity field, , corresponding to the total momentum, . This is in contrast to the multi-fluid model.35 The reduced momenta, must now be considered as internal variables, and as they are vector quantities need to obey either a covariant or contravariant dynamics to lead to materially objective equations. In this case, the covariant formalism is the appropriate one (as the fluxes are ultimately driven from the gradient of scalar potentials) leading to a Poisson bracket that is cast as follows:41,46
where the Einstein's summation convention over repeated Greek indices (from 1 to 3) is followed here and in the following. The covariant choice is also confirmed later when comparing the results to those of kinetic theory.
Fourth, as the processes of heat transfer and mass diffusion, which are an inherently irreversible process, are driven by the reduced momenta, it is in the dissipation bracket where new terms involving these new variables will appear, replacing previous contributions described directly in terms of the gradients of temperature and the chemical potentials. Therefore, this is where the extension to the bracket formulation of the standard evolution of any arbitrary functional in a multicomponent system, such as for example, described in Sec. 7.3 of Ref. 41, is necessary. The specification of the dissipation is facilitated upon introducing its primitive component. The primitive dissipation bracket, (defined as the dissipation bracket without the entropy correction term), can then be inferred as consisting of two terms, one antisymmetric (therefore leading to no dissipation and entropy correction), , and a symmetric one,
The antisymmetric dissipation bracket is cast so that it faithfully restores the individual species mass density evolution equations to their standard form when the reduced momenta are described in terms of the species mass average velocities [as duly confirmed below—see Eq. (22)]. Similarly, it also introduces the convective effects of the diffusive fluxes to the entropy flux as well as the entropy flux due to the heat flux. The contributions of the species reduced momenta coupling with the species mass balance equations are formally generated from the Poisson bracket that would have been created if we could have multiple species velocities [first term on the right-hand-side of Eq. (12)] with the mass average velocity being replaced by the species velocities and then the Volterra derivatives with respect to the species velocities being expressed in terms of the appropriate elements of the variable space , i.e., the reduced species momenta and the mass average velocity [see Appendix A, Eq. (A3)].
The contributions of the Volterra derivatives with respect to the mass average velocity are the only ones that can be accommodated as part of the Poisson bracket, Eq. (12), based on the convective, reversible dynamics, structure that they represent. We therefore need to allocate the remaining contributions to the dissipation bracket, , representing what constitutes the component of the antisymmetric portion of the dissipation bracket, , that is, related to the coupling of the species reduced momenta with the mass densities. In addition, as the diffusive fluxes also transfer entropy, additional terms need to be added coupling the species reduced momenta to the entropy density. Finally, an equivalent term for the heat flux reduced momentum needs also to appear to model heat transfer. When all these contributions are included, we find
The primitive part of the second, symmetric, contribution of the dissipation bracket, , follows from the general form, which is assumed to be a general bilinear functional with respect to F, H, near equilibrium, as follows:
where with are the shear and bulk viscosity41 and are the i-j generalized binary mass transfer coefficients, , ith generalized Soret/Dufour coefficients, , and a scalar generalized heat conductivity coefficient corresponding to this problem, respectively.
As shown in Appendix B, the expression offered in Eq. (15) is the most general which (a) is bilinear in terms of the Volterra derivatives involved and (b) only considers an isotropic behavior (thus involving only scalar transport coefficients). The latter assumption, reasonable in the case that there are no other, tensorial in nature, structural variables, may not hold when those are present, as for example in the case of polymeric systems—see, for example, Ref. 55. For simplicity, we keep it here with the caveat just noted. Moreover, please note that by construction only the symmetric, off diagonal part of the matrix of the coefficients matters, and these are the only ones considered [for a total of n(n-1)/2 independent components]. Furthermore, note that not all the generalized Soret/Dufour coefficients are independent as, due to the fact that the reduced species momenta are required to satisfy the homogeneous constraint defined by Eq. (7), they can be defined only up to a constant. To leave them fully specified, the Soret/Dufour coefficients have therefore to be required to satisfy the corresponding homogeneous (conservation of linear momentum) constraint
In addition, it is important to note that as a result of the symmetry of the transport coefficients and the constraints described by the above equation, the Onsager/Casimir reciprocal relations are duly satisfied. Moreover, for the non-negative requirement for the corresponding entropy production, all, that is, required to ensure that the quadratic expression for the local rate of entropy production (see below for the derivation) is positive semidefinite. This generates straightforward but rather complex general requirements on the transport coefficients ; suffice to mention though here that in the particular case of a binary mixture those requirements reduce to (a) the transport coefficients to be non-negative and (b) the Soret/Dufour coefficients to satisfy the inequality . We also note that general constraints can be obtained through the use, for example, of a modified Sylvester's criterion: all principal minors of the corresponding square matrix need to have non-negative determinants.56
The entropy correction to the above dissipation bracket is calculated, as usual, based on the constraint of the conservation of the overall system energy (Hamiltonian) H, as follows.41 In general, if the primitive dissipative bracket is offered as
then, the entropy correction term, is defined as
with the (full) dissipation bracket, , provided by the sum of the primitive and entropy correction terms
Corresponding to the above-described Poisson and dissipative brackets, the evolution equations for the state variables are then developed as usual,41 i.e., by requiring the total derivative of any arbitrary functional , also expressed as [see also Eq. (2)]
To be identical to the fundamental equation in the Hamiltonian bracket nonequilibrium thermodynamics, Eq. (1). Using this rationale, the expression for the total time derivative of individual component densities, , would be equal to the coefficient of in the one generator expression for , based on the evaluation of the Volterra derivatives of the Hamiltonian provided by Eq. (11), leading to
which duly reduces the equation to the familiar form of the continuity equation for a single component
Given that the mass densities are independent and additive, performing a summation over all n components in the system of the above equation gives us the overall continuity equation of the system
Through similar reasoning, one can arrive at an expression for the total time derivative of the total momentum density ; that comes out to be
where is a reversible contribution to the stress, as it arises from the reduced species momenta
and is the non-equilibrium pressure of the system defined as
Note that the above expression also accounts for purely nonequilibrium effects arising due to differences in the velocities of the components. However, it is independent of the overall velocity of the system (and therefore Galilean invariant)! Indeed, if one focuses on the contributions that the overall kinetic energy may have on the definition provided by Eq. (26) one can easily find those to be zero
This is an important finding and necessary for the well-posedness of the proposed theory. From Eq. (23) and the definition of total momentum, Eq. (4), one can get the following evolution equation for the mass average velocity:
Similarly, repeating the exercise for the system entropy density, s, we get
where represents the total rate of entropy production and is given as
As mentioned, requiring this expression to give us always non-negative results for the allowed (through constraints) variables, gives as constraints for the various transport coefficients that are involved.
Repeating the exercise for the reduced momenta densities results in the following expression:
Equation (31) can be also written in a slightly more compact form, by combining terms in the definition of the time derivative as well as those involving the Volterra derivatives with respect to the species mass densities, as
with the overscript triangle representing the generalization of a lower convected time derivative for a covariant vector quantity57 weighed by a density. This is a particular case of the general Oldroyd time derivative,46,58 defined as
Note that at equilibrium (i.e., in the absence of any flow and inhomogeneities) the values of the reduced momenta densities , and their lower convected time derivatives, are equal to .
The last equality was obtained by using the identities
Note that, as far as the explicit dependencies on the velocities are concerned, Eq. (35) duly agrees to the results from kinetic theory.59,60 This agreement would not have been possible if another form (i.e., such as that associated with a contravariant vector, or scalar) was used to represent the reversible (i.e., dictated by the Poisson bracket) time derivative rather than the covariant one adopted here. Moreover, in comparison of Eq. (35) with Eq. (130.4) in the last of Truesdell's voluminous work,17 the only difference is that we include a term in lieu of the divergence of a partial stress term which though is left undefined except for the requirement that summing all the partial stresses results in the total stress of the system. In fact, Truesdell and Noll in the last few pages of their work, following Eq. (130.4),17 attempt to make sense of that term based on the work of Adkins and coworkers, with the end result making them feel unsatisfied as it is described in the last paragraph of their work:17 “The foregoing summary makes it clear that the rational theory of diffusion is in its infancy, with many possibilities awaiting search. The uncertainty as to even the overall pattern for the right constitutive equations, reflected by the diversities among the five papers on the subject issued by the same author within a few months, suggest that some basic principle of the subject remains to be discovered.” Thus, in addition to offering a meaningful correction, a key benefit resulting from our investigation is an explicit evaluation of all stress, pressure, and force terms.
This concludes the formal description of the general flux-based heat and mass transfer equations. In Secs. III–IV, we will explore different limits of these equations to (a) validate them, (b) interpret important effects included within them as subcases, and (c) explore the symmetries and mathematical characteristics in developing efficient numerical schemes.
III. LIMITING CASE OF SMALL FLUXES: VALIDATION AND GENERALIZATION OF STEFAN–MAXWELL AND SORET/DUFOUR EQUATIONS AND AN APPLICATION TO A BINARY SYSTEM
First, it is of interest to see what the theory predicts for the limit of small differential velocities (therefore neglecting square velocity terms) and small rates of change (therefore equating to zero the corresponding Oldroyd time derivatives). In that limit, Eq. (34) becomes
where, in obtaining the second equality in the first line, the constraint given by Eq. (16) has been used. Note that, as a result of the constraint provided by Eq. (16), in the above equation one can substitute the absolute species velocities, by the relative ones through which one can easily evaluate the corresponding fluxes, through Eq. (6).
These equations bear high resemblance to the Maxwell–Stefan equations, with corrections for the coupling with heat transfer (Soret/Dufour effects). The best way to check them (as well as to get some idea on possible expressions for the transport parameters entering into these equations) is through a comparison against the results of an independently validated microscopic theory, namely the kinetic theory of gases. Indeed, the original Ref. 12 provides the following two relations for multicomponent heat and mass transfer. First, Eq. (7.6-12) of Ref. 12 provides an expression for the heat flux (after correcting for the enthalpy flux included in (7.6–12) as it represents in the work the “total energy flux”), which after rewriting it slightly to our notation is given as
where is the thermal conductivity, and the are the binary diffusion and thermal diffusion coefficients, respectively.
A physical interpretation shows the specific form of the more general transport coefficients, defined in Eq. (15), as appropriate for ideal gases. The first relation shows that the generalized heat conductivity coefficient is proportional to the entropy density squared, temperature and inversely related to the thermal conductivity. As shown in Eq. (30), it also weighs the square of the heat flux velocity in the expression for the rate of entropy production. The second set of generalized Soret/Dufour coefficients, connecting the reduced momentum densities [first term of Eq. (34) to the heat flux, they involve both the thermal diffusivities and the binary mass diffusivities]. Note that these satisfy the constraint provided by Eq. (16) and the Onsager/Casimir reciprocal relations. Importantly, mapping the general result to the specific case of an ideal gas provides explicit limiting forms for the unspecified generalized transport coefficients.
Second, Eq. (8.1-3) of Ref. 12 provides an expression for the mass transfer which after some rewriting to our notation reads as
where the vector is defined in Eq. (8.1-2) of Ref. 12 that in the absence of external forces reads as
Now, for an ideal multicomponent gas, the chemical potential of component j, , is given as [see, for example, Eq. (6.1.8) of Ref. 61]
Therefore, with the use of Eq. (42), the gradient of the chemical potential, at constant temperature, with temperature (or, equivalently, entropy) being an independent variable of the problem, is given as
leading to the following equality:
If we then use the first line in Eq. (37) that we have shown above to be equivalent to Eq. (38) derived from the kinetic theory [Eq. (7.6-12) from Ref. 12] to eliminate the dependence on the temperature gradient in Eq. (45) it yields
or, using Eq. (39) to eliminate and after multiplying both sides by we find
Finally, this can be written explicitly in terms of velocity differences by using the condition given by Eq. (16) as
Comparison of Eq. (48) with the second line of Eq. (37) shows that the two equations are fully equivalent provided that the following expression is used for the phenomenological coefficients for dilute gases:12
Nicely, the expression is symmetric as required.
Thus, in this limit the mass transfer coefficients, which by Eq. (34) relate the reduced momentum density to the species mass velocities, explicitly depend not only on the mass diffusivities but also on the Soret/Dufour and generalized heat conductivity transport coefficients. This result is to be expected as the relative mass velocities are coupled to the heat flux. This comparison shows (a) the complete equivalence between our model and the one arising from the kinetic theory of dilute gases12 and (b) identifies through the determination of the entropy production (and therefore of the non-negative definite constraints on the generalized mass transfer coefficients) the useful range where this model is thermodynamically admissible.
It remains to see now some typical values for those transport coefficients. Let us first consider the subcase of an ideal binary gas mixture. First, notice there that, because of the constraint expressed by Eq. (16), (as mentioned before), and the only nonzero mass flux coefficient following Eq. (49) is , which is always non-negative:
For dilute binary gas, the significance of the reduced coefficient emerges when one considers the only other criterion that we can develop for the entropy production as described by Eq. (30) to be non-negative besides the coefficients being non-negative
which is of course duly satisfied given the correspondence provided by Eq. (51).
It is of interest to evaluate a typical value of the correction suggested by Eq. (50) on the reduced mass transfer coefficient . To do so, let us use one of the many systems discussed in Ref. 12 for which a complete set of data on the binary diffusion, thermal diffusion, and thermal conductivity are available. A nitrogen–hydrogen, , at 65% of hydrogen binary system seems to fit that constraint: In Table 8.4–11 of Ref. 12 it is reported a thermal conductivity at T = 273.16 as , whereas in Table 8.4-12 a binary diffusivity at the same temperature is reported as . Finally, from Table 8.4-14 for the same system and at about the same temperature a value for the thermal diffusion ratio can be estimated as about . From the definition of the thermal diffusion ratio for a binary gas mixture [Eq. (8.1-13) of Ref. 12]
This results in the following relationship for the coefficient:
Finally, the total molar concentration for an ideal gas at the same temperature is . With all this information, it is straightforward to evaluate the relative correction to the reduced mass transfer coefficient
We therefore see that at least in this case this correction is small and can be neglected, thereby approximating , a practice followed in the literature and also below (whether though this is true in all cases it remains to be seen).
Another case where the equations can be tested for agreement with previous models is in the Fickian binary diffusion of a dilute solute in a liquid solvent in the absence of heat transfer and, as mentioned above, neglecting the correction to the reduced mass transfer coefficient, and using the expression appropriate for the ideal gas case. Under those circumstances, the approximations often made are
leading to the following expression for the right-hand-side and the first term of the left-hand-side of Eq. (47), respectively, as applied to that liquid binary system, for :
where is the mass-based diffusive flux for component , corresponding upon equating the two sides to the following expression for the flux:
A final exploitation of the results obtained here is in delivering an approximation where the fluxes are eliminated as independent variables of the problem but are rather reconstructed as proportional to the gradients of the corresponding potentials (typically followed in heat and mass transfer applications). The Hamiltonian bracket formalism allows us easily to obtain the corresponding equations following an approach that was first introduced systematically in Refs. 63 and 64 in connection to a dilute emulsion and taking as an example that of a damped linear oscillator. It simply suffices to (a) remove those variables from the Poisson expression, Eq. (12), (second and third line in that equation), (b) remove the asymmetric contributions, Eq. (14), from the dissipation bracket, and (c) write the remaining (symmetric) dissipation bracket in terms of the potential gradients. This last part is facilitated by using the constitutive equations provided by Eq. (37) to solve for the fluxes in terms of gradients of the potentials; however, the resulting relations are complicated for the general case. However, things are simpler in the case of the binary gas where the above relationships simplify considerably
leading to the following relations:
This results in the following symmetric dissipation bracket:
The important item to notice from this exercise is that not only do we obtain a the required bilinear form for the dissipation bracket (e.g., Ref. 65) in terms of the gradient of the chemical potential differences and gradients of the temperature; but that all transport parameters involved are now explicitly expressed in terms of the true transport parameters of the problem (i.e., ) with explicit expressions also provided by kinetic theory.12 Even if these correlations have been derived for binary mixtures before65 their correspondence is much clearer, and the approach illustrated here applies to more general, multidimensional systems. Furthermore, this new approach allows for added corrections that may be important in certain systems (such as that involved in the differentiation between the parameters and ). Moreover, this approach and the illustration here are valuable for extending the theory to binary systems applied beyond dilute gases and multicomponent cases, where it is anticipated that the unknown transport coefficients can much better fitted to experimental data using simple empirical models that modify the limiting results obtained from the kinetic theory of ideal gases, rather than starting from a general description that includes all dependences in terms of the state variables open for evaluation. We therefore believe that this connection to the kinetic theory of dilute gases is very valuable for future investigations of multicomponent heat and mass transfer in more complex systems.
IV. BINARY DIFFUSION OF A DILUTE SOLUTE IN A LIQUID SOLVENT
As an illustrative flux-based transport application, it is instructive to consider a simple time-dependent, binary dilute liquid diffusion case, allowing for the time evolution of the full time-dependent solute species momenta. The governing equations in this case are: (1) the species mass densities evolution equation, Eq. (22), applied for species A (solute) assuming zero overall flow, which becomes
and (2) the corresponding species momentum equation, Eq. (60) for a stagnant medium, for which and further, by neglecting higher order contributions, yields the following linearized equation:
with characteristic speed
i.e., a correspondence is established in this limit between solute mass transfer and the wave equation with finite speed of propagation! For a typical low molecular weight liquid (e.g., water, ) and for room temperature,
which is quite appreciable and of the same order of magnitude as the speed of sound in the medium (for water is equal to about 1500 m/s). Note that using this characteristic speed, which is comparable to the average Maxwell–Boltzmann root mean square (RMS) velocity , enables to rewrite Eq. (67) as a hyperbolic equation
which reduces to the usual, linear parabolic diffusion equation when . Therefore, with the new formalism we validate the traditional approach for most applications, except possibly those corresponding to very fast transients.
With this opportunity, we would like to clarify here that this is neither the first time to develop a hyperbolic diffusion equation to transport problems nor the above modeling procedure is the only way to accomplish it. As a matter of fact, there are many ways to obtain a hyperbolic diffusion equation in transport, the simplest one being adding a rate term to the traditional transport relations as first demonstrated for momentum in viscoelasticity by Maxwell66 and later in heat transfer by Cattaneo.67 In mass diffusion, this task has been so far mostly accomplished with the use of an internal variable—see, for example, Refs. 40 and 68 and references therein—although there have been interesting works to develop it from a microscopic point of view using Kramers equation for a Brownian particle.69 However, what we want to stress here is the universality of the approach through which this equation is obtained, applicable to all systems, with the only difference being in the values of the parameters involved. This restores physical expectations from the model, as for example that for a finite speed of signal propagation, without requiring us to impose them a priori. We feel therefore that this deserves a special discussion of some of the practical implications resulting from the use of such an equation, further illustrated below in a specific example of 1D transient mass transfer.
It is also worthwhile to explore Eq. (70), when applied to traditional diffusion problems to observe the time and length scales where the differences between species concentrations predicted by the parabolic and hyperbolic diffusion equations become prominent. It is also illustrative to review the solution methodologies that can be employed to solve the second order linear damped wave equation. Equation (70) can be rearranged to resemble the traditional parabolic diffusion equation or Fick's second law, albeit with a second order time derivative term as shown
Here, is the relaxation time scale for propagation defined as that corresponds to the response time for diffusion when a concentration gradient is imposed. The above equation resembles the damped version of Fourier's law first developed by Cattaneo67 to account for inertia arising from finite speed of propagation in heat transfer.
This problem was further explored by Sieniutycz and Berry70 using variational principles of Onsager type. These variational principles subject the system to a constraint such that the growth of the total entropy generation is minimum. By working out vector couplings occurring in simultaneous heat and mass transfer, the authors were able to obtain a complete set of hyperbolic equations describing heat and mass transfer. In a previous work71 Sieniutycz presented relaxation times for several processes and systems. The relaxation times vary by several orders of magnitude depending on the length scale of the system; for instance, in Brownian diffusion, the relaxation time ranges from to s. As the relaxation time becomes smaller, so does the difference between predictions made by Fickian and hyperbolic diffusion equations.
For the purposes of illustration, consider a dilute solute species diffusing in a solvent in one-dimension through a slab of length , where the solute enters through both boundaries (applied as a step function) at with weight fraction , with the solution initially being at solute weight fraction of . Let us further consider as spatial domain of interest the interval . For this system, the governing equation can be expressed by non-dimensionalizing Eq. (71) using the weight fraction at boundary , the initial weight fraction , the characteristic time , and the corresponding diffusion length as scales to construct appropriate dimensionless variables (denoted here with * superscripts), such that , , and .
The dimensionless system length can then be considered as another measure to characterize the system. The differences between predictions of parabolic and hyperbolic equations are significant when the . Sandler and Dahler59 argue that the differences are only significant at length scales of in liquids or in gases.
Expressing the dimensional variables in Eq. (71) in terms of the above-defined dimensionless equivalents gives us
Exploring the symmetry over the centerline, this problem can be solved in half domain and mirrored about . We can therefore replace the left Dirichlet boundary condition at with a Neumann boundary condition at such that . One additional boundary condition is required to obtain the coefficients, since the equation in second order in time, for which we propose . Solving the above equation using Fourier eigenfunctions such that , the eigenfunctions and eigenvalues, , respectively, based on this combination of homogeneous Neumann and Dirichlet boundary conditions come out to be
where . Using a linear projection of variables along the eigenfunctions, defined as , the transformed equation and boundary conditions become
The solution for this transformed equation is
which can be expressed in the dimensional space as
Figure 1 shows the above solution plotted for an initial solute weight fraction and weight fraction at the boundary . At smaller times, one can observe the progression of the solute front across the spatial dimension. The dotted lines represent the solution of parabolic (Fickian) diffusion problem for the same boundary conditions for comparison. At short times, the parabolic diffusion equation has finite value everywhere, which is nonphysical since solute species cannot travel at infinite speed. This nonphysical behavior of the Fickian diffusion can be clearly contrasted with the improved inertial diffusion model that accounts for the finite speed of propagation in medium.
Another method to obtain an analytical solution for this problem valid at small times is by first considering the transient solution in an infinite domain () of an initial delta function, i.e., subject to the initial conditions,
The solution to this problem has been developed by Gomez72 based on modified Bessel functions of the first kind. The solution was found to be
where are the modified Bessel functions of the first kind and order 0, 1, respectively, defined using a series expansion as
Given the form of initial conditions used (78), the above solution can be considered as a Green's function, and the solution for any other initial value problem
can therefore be obtained as a convolution
For the semi-infinite domain () and the Dirichlet boundary condition applied on the right boundary , the corresponding f function is
Furthermore, because for all values of the solution is equal to 0, Eq. (82) can be simplified to
Evaluating the integral using Eq. (79) gives us the solution for the semi-infinite problem in the following form:
The integrals involved in the expression can be computed either numerically (using Gaussian quadrature) or analytically using truncated series expansion and integration of polynomial expressions—see Appendix C.
Yet another approach to solve the problem of hyperbolic diffusion equation is numerically by using finite differences. In this case, we first need to identify the characteristics of this problem along which the solution propagates, namely and . The numerical solution of hyperbolic problems is restricted by the domain of dependence of the corresponding characteristics. In this case, stability requires (when explicit finite differences are used in time) the time steps () to be smaller or equal than the spatial steps (). For simplicity and for illustration purposes, we propose here a first order finite difference scheme with the time steps () equal to the spatial steps (). This can be most efficiently implemented using staggered grid points for discretization, as shown in Fig. 2. To properly apply the finite differences approximations, we first need to define the following material derivative operators:
where is computed along the (positive) characteristic, and is computed along the (negative) characteristic. Equation (72) can then be written in terms of these operators as
Using the stencil shown in Fig. 2, the above expression can be discretized as follows:
Of course, when a higher accuracy is desired and/or the resolution of additional features present in the equations (additional terms, higher dimensions, and/or nonlinear dependencies) one needs to resort to more sophisticated numerical methods. Again, the equations being hyperbolic allows for a wide selection of really efficient numerical methods—see, for example, Refs. 37 and 73. A similar simplification is carried over to the specification and implementation of the boundary/initial conditions as discussed next.
Second, the propagation from the boundaries/initial conditions of any discontinuities needs to be evaluated numerically or analytically as appropriate for the problem at hand corresponding to the characteristic structure of the governing equation, Eq. (87). For example, in the case of the solution undergoing a step discontinuity from 0 to 1 at the point , this step change propagates along the negative characteristic emanating from this point, satisfying the following equation:
which readily admits the analytical solution . After upgrading the values at the grid points due to the propagation of the boundary discontinuities, the discretized equation can then be iteratively solved for grid point recursively for each layer of time at a time.
Characteristic results are shown, in comparison to the available analytical solutions, in Figs. 3 and 4. In Fig. 3, a comparison of the solutions obtained from finite differences with the analytical solution for semi-infinite domain () obtained in Eq. (85) for various spatial resolutions. A first order convergence of the numerical solution in terms of the common length/time step size is clearly demonstrated. The finite difference and analytical (series) solution to the hyperbolic diffusion problem in a slab has also been compared with numerical and analytical (series) solutions to the parabolic diffusion problem in Fig. 4. As seen there, the numerical solution to the hyperbolic problem shows comparable error to the numerical solution of the parabolic as compared to the common hyperbolic/parabolic solution achieved at times only one order of magnitude larger than the characteristic hyperbolic time .
For transient problems involving very large length scales, the kind one would typically encounter in civil and environmental engineering, it becomes very critical to model the propagation front of diffusing species. Gomez et al.37 have shown how one can leverage the hyperbolic diffusion equation to take advantage of powerful numerical techniques like the discontinuous Galerkin (DG) approach to model convective and diffusive transport in systems involving such disparate length scales. The authors present an example of accidental spillage in a port, demonstrating the capability offered by hyperbolic equations for solving large scale problems efficiently with reduced computational effort.
The capability to express a parabolic problem in a hyperbolic form can therefore be significant even in cases where the hyperbolicity does not directly influence the results, i.e., when , as it allows for the use of efficient hyperbolic numerical solvers. Of course, when and explicit solvers are used the stability requirements on the use of smaller time steps may be excessive. However, in that case, one can use an artificially high relaxation time to alleviate those restrictions, as long as its use does not affect appreciably the results. This choice is not different from the use of artificial diffusion or viscosity in complex fluid flow problems (like viscoelastic and/or turbulent74,75)
The linear hyperbolic diffusion equation, as explained before, is not new—it has appeared several times before in the literature been derived using a variety of approaches, ranging from kinetic theory,59,60 to variational approaches.70,71 What is new though is the systematic development of this equation within the context of nonequilibrium thermodynamics and the capability to extend it to accommodate nonlinear, flow, and/or multicomponent species effects.
In this work, we demonstrated how the use of a systematic nonequilibrium thermodynamics approach, based on four simple steps (selection of state variables, description of the total energy (Hamiltonian), definition of the reversible (Poisson), and irreversible (dissipative) bracket dynamics), can systematically lead without ad hoc assumptions to the complete model equations of complex fluid flow, heat, and mass transfer, including the effects arising from the inertial dynamics of the mass fluxes, within a multicomponent system. It thus establishes to the fullest the analogy and unification between all transport processes, momentum, heat, and mass transfer, toward which R. B. Bird contributed so much with his work,21 making the present work an appropriate tribute to the celebration of his life. Simultaneously, we hopefully showed how his work could be amplified with judicious help from nonequilibrium thermodynamics. Although Bird was clearly aware and a good user of the linear irreversible thermodynamics [see, for example, the last chapter, Chap. 24, of Ref. 21 and an expert in the use of equilibrium thermodynamics to transport processes (as his last paper,15 easily demonstrates) curiously he had not exploited by himself much of the power of the reversible–irreversible structure of either SGBF-NET or GENERIC, albeit he did quote those works in the above-mentioned Chap. 24 of Ref. 21, most likely foreseeing their positive future impact to this field.
Several additional accomplishments are to be noted: First, we demonstrate how to correctly incorporate inertia (flux) effects in the modeling of coupled heat and mass transfer in a method consistent with nonequilibrium thermodynamics while preserving Galilean invariance. This approach incorporates species flux effects through the explicit use of their reduced momenta densities as additional variables. Thereby, this approach resolves a long-standing conundrum in nonequilibrium thermodynamics. Second, this modeling is shown to be successful in extending traditional descriptions of multicomponent diffusion–convection problems leading naturally to general expressions for diffusion and mass transfer without ambiguities or arbitrary assumptions. Those expressions correctly reduce to the extended (including heat transfer) Maxwell–Stefan equations and Fickian diffusion in the corresponding inertial-less limits of ideal gas and binary liquid with dilute solute, respectively. A successful comparison against kinetic theory predictions, in addition to providing for an independent validation of the modeling approach, enables obtaining explicit expressions for all the transport coefficients involved for this limiting behavior, while also providing a paradigm for extending the kinetic theory results to denser and multicomponent systems.
There are many open problems in multicomponent systems where the above-presented analysis may find important application. First, is the issue how to deal effectively with nonidealities and fugacities in extending the dilute gas kinetic theory relations for the transport coefficients. Moreover, as multicomponent systems present critical points where possibly molecular or thermal diffusion coefficients can vanish there is controversy how this can affect the transport properties and what are the limitations in observed behavior given the Onsager reciprocal relations.
Working directly with the mass fluxes (reduced momenta) and their dynamics has distinct advantages: On the one hand, this approach represents the physics more faithfully, allowing for finite speeds of propagation for the disturbances as well as for a wider variety of boundary conditions (described now in terms of the vector fluxes, not on scalar concentrations). On the other hand, it leads to fully hyperbolic equations that are easier to solve numerically through the use of powerful hyperbolic numerical solvers. The traditional parabolic convection–diffusion equation can be obtained in the limit of small relaxation times from its fully hyperbolic, finite signal propagation, equivalent. It is therefore proposed and illustrated here for a specific simple example, how to exploit the hyperbolic formulation resulting from the flux modeling of diffusion–convection problems in order to apply powerful numerical hyperbolic solvers to obtain solutions in their parabolic limits. Additional applications in many areas of chemical engineering are anticipated.
The authors acknowledge funding for this research provided by the National Science Foundation through CBET No. 1804911.
Conflict of Interest
The authors have no conflicts to disclose.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
species/component in the system
conformation tensor corresponding to component
molecular diffusion coefficient between species
Thermal diffusion coefficient of species
Hamiltonian functional density
mass-based diffusive flux of species
thermal diffusion ratio
molar mass/molecular weight
viscous dissipation tensor
universal gas constant
total entropy density
temperature of the system
total momentum density
species momentum density,
species reduced momentum density,
heat momentum density
velocity field of the system or reference frame
species reduced velocity,
weight fraction of species A
molar fraction of species
generalized coefficient of mass transfer between components and ,
reduced mass transfer coefficient,
Soret–Dufour coefficient for component ,
generalized heat conductivity coefficient
energy density functional
species chemical potential
reference chemical potential of species
primitive dissipation functional
total system density
effective heat mass density
volumetric rate of entropy production
characteristic time in hyperbolic diffusion
domain of interest in space
APPENDIX A: CORRESPONDENCE BETWEEN THE OVERSCRIPT CURVED ARC AND REGULAR PRIMITIVE VARIABLES AND THEIR CORRESPONDING VOLTERRA DERIVATIVES OF AN ARBITRARY FUNCTIONAL G
For convenience and easy reference, we summarize here the correspondence between the two alternative and equivalent frameworks. First, it is clear that only the momentum density variables are affected. Second, from the definition of the overscript curved arc variables, Eqs. (6) and (5), their correspondence with the original momentum species densities is as follows:
Similarly, the correspondence between the original and overscript curved arc variables is
From those two sets of relationships, the following relationships can be derived between the corresponding Volterra derivatives of an arbitrary functional G using differentiation by parts:
APPENDIX B: EQUIVALENCE BETWEEN THE MOST GENERAL AND THE PROPOSED BILINEAR DESCRIPTIONS FOR THE (PRIMITIVE) SYMMETRIC DISSIPATION BRACKET IN MULTICOMPONENT DIFFUSION
Considering the variables , the symmetries (parity) upon time reversal, and their physical interpretation [i.e., the definition of the reduced momenta based on the difference of the species velocity from the mass average velocity, Eq. (6)] the most general bilinear description for the (primitive) symmetric dissipation bracket for an isotropic medium is
where with are the shear and bulk viscosity41 and are the i-j generalized binary mass transfer coefficients and a scalar heat conductivity corresponding to this problem. With the caveat that when anisotropicity is present the need to be substituted by symmetric second order tensors, while can have a more complex, anisotropic, structure, the above expression can easily be shown to present the most general bilinear description involving the above-mentioned variables. Of course, higher order dependencies on the Volterra derivatives of the Hamiltonian are also possible to appear in the dissipation bracket, still, close to equilibrium, upon linearization, those expressions also reduce to the above one.
At first sight, Eq. (B1) is different from the expression ultimately used, Eq. (15). We will show in this appendix though that they are equivalent when one considers that the reduced momenta are not independent but that they have to satisfy the constraint specified by Eq. (7). That constraint is easy to show that can also be described by the inner product of the n-vector with the unit n-vector where the Euclidean L2 norm is used for normalization
i.e., the n-vector is restricted to be in a subspace that is perpendicular to with the unit n-vector . Due to this constraint, it is well known that a similar condition also applies41 to the n-vector of the Volterra derivatives with respect to the reduced momenta of any arbitrary functional G (G = H, F, ….),
i.e., the n-vector is restricted to be in a subspace that is perpendicular to the unit n-vector .
Due to this constraint, not all of the n(n + 1)/2 components of (a symmetric, due to Onsager/Casimir reciprocity relations4–6) n × n coefficient matrix are independent. Instead, the only relevant information (that can lead to measurable results) belongs to the subspace of the n × n coefficient symmetric matrix that satisfy the corresponding to Eq. (B3) condition
i.e., without loss of generality we can assume that the n × n coefficient symmetric matrix has a zero eigenvalue along the unit n-vector which becomes one of its eigenvectors.
The most general expression that satisfy this condition can be shown to be provided (through a partial eigenvalue analysis, always possible for a symmetric matrix) by
where is a set of n-1 independent orthonormal n-vectors perpendicular to the unit n-vector , the n-1 zero vector and a general, symmetric matrix. It can easily be confirmed that constructed as indicated by Eq. (B5) duly satisfies the conditions of Eq. (B4).
Although Eq. (B5) is helpful in identifying one way to construct as well as in constructively showing that the actual space where belongs has n(n-1)/2 dimensions [the number of independent components of the general, symmetric matrix ] it does not offer neither the easiest nor the most physical way to define it as its components by necessity depend on the selection of the orthonormal vectors . However, once we have identified, thanks to Eq. (B5) the dimensionality of the space of suffices in order to uniquely define it to decompose it in terms of any set of n(n-1)/2 independent members of that space. As it turns out, those are most naturally constructed based on the decomposition offered in Eq. (15) using the corresponding to each coefficients n × n elementary matrices defined through the specification of their k, l component, , as
It is indeed straightforward to show that each one of those matrices satisfies the conditions corresponding to Eq. (B4) (i.e., it belongs to the right space) and that for they are independent, and therefore, as they are n(n-1)/2 in number they form a basis for that space in terms of which the generic matrix can be described, as
APPENDIX C: THE INTEGRALS APPEARING EQ. (85) CAN BE SOLVED ANALYTICALLY BY FIRST PERFORMING SERIES EXPANSION OF THE MODIFIED BESSEL FUNCTIONS OF THE FIRST KIND, FOLLOWED BY INTEGRATION OF THE POLYNOMIAL TERMS THAT APPEAR IN THE SERIES EXPANSION
For the first integral term , the expansion can be expressed as follows:
Similarly, the second integral term can be expressed as
The above integrals can be used to express Eq. (85) as a series expansion in the following manner: