Finite Amplitude Stability of Internal Steady Flows of the Giesekus Viscoelastic Rate-Type Fluid

Using a Lyapunov type functional constructed on the basis of thermodynamical arguments, we investigate the finite amplitude stability of internal steady flows of viscoelastic fluids described by the Giesekus model. Using the functional, we derive bounds on the Reynolds and the Weissenberg number that guarantee the unconditional asymptotic stability of the corresponding steady internal flow, wherein the distance between the steady flow field and the perturbed flow field is measured with the help of the Bures–Wasserstein distance between positive definite matrices. The application of the theoretical results is documented in the finite amplitude stability analysis of Taylor–Couette flow.


Introduction
Flows of viscoelastic fluids exhibit the phenomenon dubbed "elastic turbulence" or "inertia-less turbulence". The flows of viscoelastic fluids can become-unlike the flows of the standard viscous fluids-unstable or "turbulent" at very low values of the Reynolds number. This behaviour indicates that the instability or transition to "turbulence" is driven by a nonstandard mechanism. Namely, it is not driven by the nonlinearity due to the inertial term in balance of linear momentum, but it must be attributed to the nonlinearity in the governing equation for the "elastic" part of the Cauchy stress tensor.
The key challenge is to identify the parameter values that prohibit the onset of "elastic" instability or that trigger the "elastic" instability, and to describe the transition scenarios leading from the laminar to the "turbulent" flow. This task requires one to perform some sort of nonlinear stability analysis, since a nonlinear interaction between the finite amplitude perturbations might be decisive.
The phenomenon of "elastic turbulence" has been thoroughly investigated both from the experimental as well as theoretical point of view (see reviews by Petrie and Denn [1], Larson [2], Shaqfeh [3], Morozov and van Saarloos [4] or Li et al. [5]). In particular, the experiments reported by Groisman and Steinberg [6] stimulated enormous research activity regarding the elastic turbulence. On the other hand, theoretical results mainly follow from direct numerical simulations based on various viscoelastic rate-type models (see Dubief et al. [7], Lieu et al. [8], Grilli et al. [9], Page and Zaki [10], Biancofiore et al. [11], Valente et al. [12], Lee and Zaki [13] and Plan et al. [14] for some recent contributions). The need to resort to sophisticated numerical simulations in order to get qualitative insight into the flow dynamics is not surprising.
The reason is that the instabilities in viscoelastic fluids are very likely of subcritical nature (see Meulenbroek et al. [15]). The subcritical nature of the instability implies, as remarked by Morozov and van Saarloos [4], that [Linear stability] (if it exists) is not very relevant for the existence of dynamics of the patterns that typically arise before the instability of the base state occurs.
(Note that a similar issue arises even for the standard Navier-Stokes fluid (see, for example, Baggett and Trefethen [16] for a discussion of several low-dimensional models of subcritical transition).) This means that the linear stability analysis, that is stability analysis with respect to infinitesimal perturbations, is of limited applicability in the investigation of the transition scenarios, albeit it can still provide important insight into the problem. (See, for example, Beris et al. [17], Blonce [18], Grillet et al. [19] and Pourjafar and Sadeghy [20], Pourjafar and Sadeghy [21] for linear stability analysis of flows of viscoelastic fluids described by the Giesekus viscoelastic rate-type model.) Moreover, quoting again Morozov and van Saarloos [4] [Subcritical instability] is governed by all kinds of nonlinear self-enhancing interactions and so there is almost never a simple approximation scheme that allows one to explore the infinite dimensional space of interactions in all details, and determine which direction corresponds to the smallest threshold [for instability]. Thus, in practice, one can explore such situations, in theoretical studies as well as in experiments, only for a given class of perturbations.
On top of that, even if the technique such as weakly nonlinear analysis is apparently successful, then, as Meulenbroek et al. [15] put it, One should also keep in mind that our expansion is only carried out to lowest order in the nonlinearity, so one may wonder about the robustness of these results as long as higher order terms in the expansion are unknown.
In what follows, we want to address the lack of analytical results for the stability problem of flows of viscoelastic fluids subject to finite amplitude perturbations. In particular, using a Lyapunov type technique, we investigate the stability of internal steady flows of viscoelastic fluids described by the Giesekus model, and we derive bounds on the values of the Reynolds number and the Weissenberg number that guarantee the flow stability subject to any (finite) perturbation. The result provides a sufficient condition for stability, hence it can be seen as complementary result to the search for the smallest threshold for instability via approximation methods. The derived bounds are interesting not only on their own. What is perhaps equally interesting is the way the bounds are derived. The derivation heavily relies on the underlying thermodynamical arguments and the notion of energy, which is an approach that seems to be discouraged in the nonlinear stability analysis of viscoelastic fluids (see Doering et al. [22]).
The paper is organised as follows. In Section 2, we describe the Giesekus model, and we briefly comment on its thermodynamical underpinnings. In particular, we identify the energy storage mechanisms and the entropy production mechanisms that are implied by the evolution equations for the Giesekus model. Once the thermodynamical background is summarised, we formulate the governing equations for an internal steady flow (see Section 3), and we proceed with the stability analysis of this non-equilibrium steady state. The stability is analysed using a Lyapunov type functional V neq constructed by the thermodynamically based method proposed by Bulíček et al. [23]. The functional used in the stability analysis of a steady flow v in a domain Ω is constructed in Section 4, and it is given by the formula (1) where v denotes the perturbation of the velocity field, B κ p(t) is related to the perturbation of the stress field, and B κ p(t) is related to the stress field in the steady flow. (See the corresponding sections for the notation.) The fact that Equation (1) can serve as a Lyapunov type functional is closely related to the proper choice of the distance function that characterises the proximity of the perturbation W and the corresponding steady state W. In our case, the distance function is introduced using the Bures-Wasserstein distance dist P(d), BW (A, B) = def Tr A + Tr B − 2 Tr A (see Bhatia et al. [24]), which measures the distance between the symmetric positive semidefinite matrices A and B. Once we generalise Equation (2) to the setting of spatially distributed fields of symmetric positive semidefinite matrices, we exploit the concept of Lyapunov functional, and we derive bounds on the Reynolds number and the Weissenberg number that guarantee the flow stability with respect to any perturbation (see Theorem 1). These bounds are universal for any flow geometry.
The bounds on the Reynolds number and the Weissenberg number are then explicitly evaluated in Section 6 in the case of Taylor-Couette type flow. Further, we also perform direct numerical simulations that allow us to quantitatively document some features of the perturbation dynamics. The results are commented in Section 7.

Governing Equations
The governing equations for the Giesekus fluid (see Giesekus [25]), in the absence of external force, read div v = 0, (3a) where v denotes the velocity, ρ denotes the density and B κ p(t) is an extra tensorial quantity whose physical meaning is given below. Finally, the symbol T denotes the Cauchy stress tensor that is given by the formulae T = mI + T δ , where m denotes the mean normal stress (pressure) and D = def 1 2 (∇v + (∇v) ) denotes the symmetric part of the velocity gradient. Symbols ν, ν 1 , µ and α, α ∈ (0, 1) denote material parameters. Note that, if α = 0, then one recovers the standard Maxwell/Oldroyd-B models. The value α = 0 is however not covered in the presented stability analysis.
The remaining notation is the standard one: d denotes the upper convected derivative, where L = def ∇v, and the symbol A δ = def A − 1 3 (Tr A) I denotes the traceless part of the corresponding tensor. In virtue of the incompressibility constraint in Equation (3a), we have D δ = D. Note that, if one uses a simple substitution S = def µ(B κ p(t) − I), and if one redefines the pressure, p = def −m + 1 3 (Tr S) I, then Equation (3c) and (3d) transform to λ S + S + αλ µ . This is another frequently used form of the governing equations for the Giesekus fluid.

Thermodynamic Basis
The Giesekus model has been originally derived without any reference to thermodynamics. However, we want to design a Lyapunov type functional using concepts from non-equilibrium thermodynamics, hence we need to explore thermodynamical underpinnings of the model. The issue of finding a thermodynamic basis for viscoelastic rate-type models is claimed to be resolved by a plethora of theories for thermodynamics of complex fluids (see, for example, Leonov [26], Mattos [27], Wapperom and Hulsen [28], Dressler et al. [29] or Ellero et al. [30]). (Notably, the treatise by Dressler et al. [29] contains a rich bibliography, and describes the issue from the viewpoint of the GENERIC formalism (see Grmela and Öttinger [31], Öttinger and Grmela [32] and Pavelka et al. [33])). In the present analysis, we exploit the approach proposed by Rajagopal and Srinivasa [34] that is relatively simple and that provides one a purely phenomenologically based concept of visco-elastic response.
The fact that the approach by Rajagopal and Srinivasa [34] closely follows the phenomenological concept of visco-elastic material is best seen in the interpretation of the quantity B κ p(t) that appears in the formula for the Cauchy stress tensor. This quantity can be interpreted as the left Cauchy-Green tensor associated with the elastic part of the fluid response. Using the approach by Rajagopal and Srinivasa [34], the derivation of Maxwell/Oldroyd-B type models was discussed by Málek et al. [35] (see also Hron et al. [36]). More complex viscoelastic rate-type models that document the applicability of the approach in more involved settings are discussed in Málek et al. [37], Málek et al. [38] and in Dostalík et al. [39]. Following Dostalík et al. [39], we know that the Giesekus fluid is a fluid with the specific Helmholtz free energy ψ in the form where θ denotes the absolute temperature, θ ref denotes a constant reference temperature, c iNSE V is a positive material parameter (specific heat capacity at constant volume) and µ is another positive material parameter. The specific Helmholtz free energy describes the energy storage ability of the fluid, and the chosen ansatz is the same as for the standard Maxwell/Oldroyd-B fluid. This implies that the Giesekus fluid and Maxwell/Oldroyd-B fluids differ, from the perspective of the current approach, only in their entropy production mechanisms (see below).
Specifying the Helmholtz free energy as a function of θ and B κ p(t) , one can use the standard thermodynamical identities η = − ∂ψ ∂θ and e = ψ + θη, and obtain explicit formulae for the specific entropy η, and the specific internal energy e Note that adding the kinetic energy to the mechanical part of the internal energy e, that is to the , we can define the specific mechanical energy via Once the Helmholtz free energy, and consequently also the internal energy, is specified, one can derive the evolution equation for the entropy that has the structure ρ dη dt + div j η = ξ, where j η denotes the entropy flux and ξ stands for the entropy production. In the case of Giesekus fluid the entropy production is given by the formula ξ = ζ θ , where (We use the notation A : A = def Tr AA for the scalar product on the space of matrices, and |A| for the corresponding Frobenius norm.) Since B κ p(t) is a symmetric positive definite matrix, it is easy to check that the entropy production in Equation (7) is a nonegative quantity, hence the second law of thermodynamics is satisfied. The fact that B κ p(t) is a symmetric positive definite matrix follows directly from the governing equations in Equation (3) via an argument similar to that of Boyaval et al. [40]. It is also a consequence of the fact that B κ p(t) is in the approach by Rajagopal and Srinivasa [34] constructed as the left Cauchy-Green tensor, which means that B κ p(t) can be decomposed as B κ p(t) = F κ p(t) F κ p(t) . We exploit the positivity of B κ p(t) quite frequently in our analysis.
Finally, we introduce three more important quantities that play crucial role in the construction of Lyapunov type functional via the method proposed by Bulíček et al. [23]. Namely, we introduce the net total energy E tot , the net mechanical energy E mech and the net entropy S of the fluid occupying the domain Ω,

Scaling
The equations in Equation (3) governing the evolution of mechanical variables can be transformed to a dimensionless form by introducing the characteristic length x char , characteristic time t char . (Note that the tensor field B κ p(t) already is a dimensionless quantity.) Using the following relations between the original quantities and their dimensionless versions denoted by stars where the dimensionless Cauchy stress tensor T is given by In Equation (9) . It remains to introduce a scaling factor for the net mechanical energy E mech , which is used for the construction of the Lyapunov type functional in Section 4. Using the Hereafter, we omit the star denoting dimensionless quantities unless otherwise specified. The scaling is chosen in such a way that, if Wi → 0+, then B κ p(t) approaches the identity tensor. Indeed, if Wi → 0+ then, Equation (9c) implies that and the solution of Equation (11) is B κ p(t) = I. Moreover, if B κ p(t) = I, then the second term in Equation (10), that is the elastic contribution to the mechanical energy, vanishes, and the mechanical energy of the fluid reduces to the kinetic energy only. Finally, if B κ p(t) = I, then the additional term in the Cauchy stress tensor in Equation (9d) vanishes. This means that for Wi → 0+ the governing equations in Equation (9) reduce to the standard incompressible Navier-Stokes equations.

Boundary Conditions
The governing equations in Equation (9) must be supplemented with boundary conditions for the velocity v. We are interested in internal flow problems, where one prescribes Dirichlet boundary conditions on a part of the flow domain Ω ⊂ R 3 , and periodic boundary conditions on another part of the domain. Such a domain is usually called the periodic cell. (For example, in the case of flow in between two infinite concentric rotating cylinders, the Dirichlet boundary condition is prescribed on the surfaces of the cylinders, while the periodic boundary condition is prescribed in the direction of the axis of the cylinders.) On the parts of the boundary corresponding to the periodicity directions, say Γ 1 , we therefore prescribe periodic boundary condition for v, while, on the remaining part of the boundary, say Γ 2 , we prescribe the no-penetration and the no-slip boundary condition, where n is the unit outward normal to the boundary of Ω and V is a given velocity in the tangential direction to the boundary. This means that the fluid adheres to the boundary, and, moreover, if V = 0, then, in general, the energy is exchanged between the fluid and its surroundings. Indeed, the balance of the net total energy reads dE tot dt = ∂Ω (Tv) • n ds − ∂Ω j q • n ds, where ∂Ω denotes the boundary of the domain Ω and j q denotes the heat flux. Consequently, if v = 0 on the boundary, then the term ∂Ω (Tv) • n ds does not, in general, vanish or is compensated by the second term on the right-hand side, and the net total energy might even change in time.
Concerning the boundary conditions for the perturbation v with respect to the reference state v (see below), we see that, if v satisfies Equation (12), then the perturbed state v = v + v also satisfies Equation (12) provided that v| Γ 2 = 0.
The periodic boundary condition on Γ 1 is preserved for the perturbation v.
In the following, we frequently use the identity where f : ∂Ω → R 3 is a smooth function such that f fulfills the periodic boundary condition on Γ 1 and f = 0 on Γ 2 . Note that the identity holds even if one part of the boundary, whether Γ 1 or Γ 2 , is not present.

Notation for the Stability Analysis
We are interested in the evolution of the triplet W = def [v, m, B κ p(t) ], which solves the evolution equations in Equation (3). We further use the notation W = [ v, m, B κ p(t) ] for the triplet corresponding to a non-equilibrium steady state solution, and W = [ v, m, B κ p(t) ] for the perturbation with respect to the non-equilibrium steady state. This means that the triplet describing the complete perturbed state is given as a sum of the reference state W and the perturbation W with respect to the reference state (Note that sometimes we work only with the pair W = def [v, B κ p(t) ], since the pressure is insubstantial in our analysis.) The term non-equilibrium steady state is chosen in accordance with the practice in thermodynamics, and it means that the entropy is produced at the steady state W. In particular, the adjective non-equilibrium does not refer to the stability of the steady state.

Governing Equations in a Steady State
The steady state W = [ v, m, B κ p(t) ] whose stability we want to investigate is a solution to the equations in Equation (9) where the partial time derivatives are identically equal to zero. In particular, we assume that the state described by the triplet [ v, m, B κ p(t) ] solves the system subject to the boundary conditions in Equation (12) and the periodic boundary conditions on Γ 1 . Here, the symbol T( W ) denotes the Cauchy stress tensor induced by the triplet [ v, m, B κ p(t) ], that is where D = 1 2 ( L + L ), and L = ∇ v. Note that, if V = 0, that is if no mechanical energy is supplied to the fluid, then the system would admit an equilibrium solution where c is an arbitrary number. (This is the standard ambiguity in the identification of the pressure well known from the case of Navier-Stokes fluid.) Here, we use the adjective equilibrium to emphasise that such a steady state would lead to zero entropy production. Indeed, if B κ p(t) = I and v = 0, then the (mechanical part) of the entropy production in Equation (7) vanishes.
On the other hand, if V = 0, then one must in general expect that the steady fields v and B κ p(t) are spatially inhomogeneous, and consequently the entropy production in Equation (7) is positive.
This means that the system produces the entropy; hence, it is, from the thermodynamical point of view, out of equilibrium. Consequently, as discussed above, we use the adjective non-equilibrium and we refer to the base flow as of non-equilibrium steady state.

Concept of Stability
Concerning the stability of the non-equilibrium steady state, we are interested in its asymptotic stability. If we have a non-equilibrium steady state W that solves Equation (16), then we want to know whether the perturbation W = W + W of the non-equilibrium steady state W tends back to the non-equilibrium steady state W as time goes to infinity. In our case, the evolution of the perturbed state W is governed by the equations in Equation (9) that must be solved subject to the given boundary conditions in Equation (12) and subject to initial conditions The non-equilibrium steady state W is said to be asymptotically stable provided that the solutions that start close enough to the steady state not only remain close enough to the steady state but also eventually converge to the steady state, that is the triplet W converges to W as time goes to infinity, for all sufficiently small initial data v 0 and B κ p(t) 0 (see (Henry [41], Definition 4.1.2) or Yoshizawa [42] for the formal definition and further discussion). Ideally, one would like to obtain stronger results. Namely, one would like to have an unconditional result that states that the non-equilibrium steady state is recovered as time goes to infinity regardless of the choice of initial perturbation. This behaviour is expected if one deals with non-equilibrium steady states that are driven by a small energy inflow that is by a small boundary velocity V , or, in other words, if one deals with non-equilibrium steady states that are not far away from the equilibrium steady state.
The key task in the stability analysis is the choice of a metric/norm on the state space to give a meaning to the statement in Equation (21). Namely, we need to answer the question as how to characterise the distance between W and W, since Equation (21) where dist (·, ·) is a given metric that is not necessarily induced by a norm. Since B κ p(t) is at a given spatial point x ∈ Ω a positive definite matrix, it seems reasonable to design the metric in such a way that it reflects this fact. This means that we have to rely on a metric on the set of positive definite matrices. There are several possible definitions of the metric on these sets (see Appendix A). If we use the Bures-Wasserstein distance, (see Equation (A1) in Appendix A), and if we generalise this concept to the spatially distributed tensor fields, then we can define the distance between W and W as (see Equation (A7) and Appendix A for a discussion of the notation and correctness of this definition). It turns out that this concept of distance nicely fits to the dynamical system we are interested in. The term "stability" is used in many other contexts; hence, we briefly comment on these other notions of stability. In particular, we emphasise what is in the present work not meant by the stability. First, we are not interested in the stability in the sense of continuous dependence on initial data, which is the concept of stability investigated in Dafermos [43] and various subsequent works especially in the theory of hyperbolic systems (see Dafermos [44]). The stability in the sense of continuous dependence on initial data means (see, for example, Schaeffer and Cain [45]) that [. . . ] if the initial data for an initial value problem are altered slightly, then the perturbed solution diverges from the original solution no faster than at a controlled exponential rate.
Apparently, the asymptotic stability we are interested in is a more ambitious concept, since we want the perturbed solution to converge back to the original solution (non-equilibrium steady state). Second, we are not interested in the stability of the steady state subject to infinitesimal perturbations, that is in the linearised stability. We are interested in the evolution of finite amplitude perturbations.
Finally, we emphasise that in our analysis we work with perturbations that are solution to the governing equations in the classical sense. (All derivatives are understood as the classical derivatives, not as generalised derivatives such as distributional derivatives and so forth.) In particular, we do not consider the perturbations that solve the governing equations only in a weak sense, although it is an important issue worth further investigation. The reader interested in the discussion of the state-of-the-art rigorous mathematical theory of equations governing the motion of viscoelastic fluids is kindly referred to the work of Masmoudi [46] or Barrett and Süli [47].

Concept of Lyapunov Functional
Let us briefly recall the concept of Lyapunov functional (see Henry [41]). We consider a system of governing equations in the form where X is a steady state, that is F( X) = 0, and where · st denotes a norm on the underlying state space. We say that the functional V ( X X) is a strict Lyapunov functional of the steady state X provided that: 1. There exists a neighbourhood of X such that the functional is bounded from below by a function f of the distance between the steady state X and the perturbation X, that is where f is a continuous strictly increasing function such that f (0) = 0 and f (r) > 0 whenever r > 0. 2. The time derivative of V ( X X) is negative and bounded from above by a function g of the distance between the steady state X and the perturbation X, that is where g is a continuous strictly increasing function such that g(0) = 0 and g(r) > 0 whenever r > 0.
If the given system of governing equations admits a strict Lyapunov functional near the state X, then we know that the steady state X is asymptotically stable (see (Henry [41], Theorem 4.1.4)). This means that the solution X = X + X that starts in the neighbourhood of X satisfies While the concept of Lyapunov type functional is very simple, it is difficult to apply in a particular setting. The main difficulty is to find a Lyapunov type functional. (Note that in the infinite dimensional setting one can not easily exploit LaSalle's invariance principle since it requires precompactness of the trajectories, which is a qualitative property that goes beyond our assumption regarding the existence of the classical solution. Consequently, having a Lyapunov type functional with Equation (26b) replaced by the mere negativity everywhere except at the equilibrium, d dt V ( X X) < 0 for X = 0, is not a viable option.) Fortunately, since we are interested in equations describing a physical system, we can try to search for the functional using physical concepts.
If we were dealing with the stability of a homogeneous equilibrium steady state in a thermodynamically isolated system, then a Lyapunov type functional could be constructed using the net entropy S and the net total energy E tot functional. It can be shown that the appropriate Lyapunov type functional is in this setting reads where θ is the spatially homogeneous temperature in the equilibrium steady state (see Bulíček et al. [23] for details). The observation that Equation (28) can be used as a suitable Lyapunov type functional for stability analysis of homogeneous equilibrium steady states in a thermodynamically isolated systems or mechanically isolated systems immersed in a thermal bath is well known (see, for example, Šilhavý [48] or Grmela and Öttinger [31], Öttinger and Grmela [32]), and in the continuum thermodynamics setting it dates back to the works of Coleman [49], Gurtin [50], and Gurtin [51]. (The core idea can be found in earlier works, see especially Duhem [52].) Unfortunately, the same functional cannot be used in stability analysis of non-equilibrium spatially inhomogeneous steady states in thermodynamically open systems. This fact is clear from Equation (28) itself. If one works with a spatially inhomogeneous steady states, then θ in Equation (28) is a function, and Equation (28) does not define a functional at all.

Construction of Lyapunov Type Functional for Stability Analysis of a Spatially Inhomogeneous Steady State
Recently, Bulíček et al. [23] proposed a method for construction of Lyapunov type functionals for stability analysis of non-equilibrium spatially inhomogeneous steady states in thermodynamically open systems. In the ongoing analysis, we use the same ideas as in Bulíček et al. [23], but we restrict ourselves to the mechanical quantities only. This is a matter of convenience, since we are interested in mechanical quantities only, and the temperature evolution has no feedback on the mechanical part of the system of governing equations. Consequently, we do not need to work with the Lyapunov type functional for the full thermomechanical problem, and we can construct a simpler Lyapunov type functional solely for the mechanical quantities.
Using the net mechanical energy functional E mech introduced in Equation (10), one can see that the net mechanical energy in a thermodynamically closed system must decay in time. Consequently, the functional can serve as a Lyapunov type functional for stability analysis of equilibrium spatially homogeneous state in Equation (19). Following the methodology outlined in Bulíček et al. [23], we use the Lyapunov type functional for the equilibrium steady state in Equation (29), and we define the candidate for the Lyapunov type functional for the non-equilibrium steady state as where (This is essentially the affine correction trick introduced in a different context by Ericksen [53] (see also Dafermos [43]).) The (dimensionless) explicit formulae for the individual terms in Equation (30) read Using Equation (31) in Equation (30), we get, after some algebraic manipulations, the explicit formula for the proposed Lyapunov type functional It remains to show that the functional in Equation (32) has the properties in Equation (26) introduced in Section 4.1. First, we show that the condition in Equation (26a) holds for a neighbourhood of W, which means that we have to specify a norm on the state space.
The suitable norm is the norm introduced in Appendix A in Definition A3. This norm is a "standard" Lebesgue type norm. However, it turns out that it is convenient to use this norm for a characterisation of the evolution of a "shifted" state. The idea is the following. If we are given a constant-in-time tensor field B κ p(t) , which corresponds to the steady solution of Equation (9), and a state This shifted state seems to be ideal for the investigation of the perturbations to the steady state B κ p(t) . Indeed, the steady state for Equation (9) Now, instead of investigating the behaviour of the perturbation B κ p(t) with respect to B κ p(t) , we can investigate the behaviour of the shifted perturbation Z κ p(t) with respect to the identity I. Lemma 1 (Relation between the proposed Lyapunov functional and a norm). Let W and W = W + W denote two states of the system governed by equations in Equation (9), and let Z and Z denote the corresponding shifted states. Furthermore, let · st denote the norm introduced in Definition A3. Then, there exists a positive constant D(Ξ) such that Proof. We note that Equation (32) for the Lyapunov type functional can be rewritten as where we use the cyclic property of the trace and the properties of the determinant. Making use of Lemma A3 and the inequality in the integrand of the last term in Equation (36), we see that The less straightforward part of the analysis of properties of proposed Lyapunov type functional V neq is the evaluation of its time derivative dV neq dt . The formula for the time derivative is derived via a lengthy algebraic manipulation described in Appendix B, and it is given below in Lemma 2. Note that, although we are working with a thermodynamically open system, the formula for the time derivative does not contain boundary terms. Lemma 2 (Explicit formula for the time derivative of the Lyapunov type functional). Let W and W = W + W denote two states of the system governed by equations in Equation (9), where the state W is the steady state, that is it solves Equation (16). The formula for the time derivative of the functional V neq ( W W ) introduced in Equation (32) reads We note that the terms are negative provided that v = 0 and B κ p(t) = O. If we were able to show that these damping terms are strong enough to balance all the remaining terms on the right-hand side of Equation (38), we would get the desired result (Equation (26b)) concerning the negativity of the time derivative. This should be possible at least for sufficiently small Reynolds number Re and Weissenberg number Wi. The hypothesis follows from the observation that as Re and Wi tend to zero, then the magnitude of the damping terms increases, and it should outgrow the other terms in Equation (38) that do not depend on Re and Wi. This observation is consistent with the expectation that low Reynolds number and low Weissenberg number flows are stable. Now, the objective is to show that the hypothesis is true, and that the proposed functional indeed satisfies the condition in Equation (26b). In the quantification of the "sufficient smallness" of the Weissenberg number Wi and the Reynolds number Re, we aim at a simple but very rough estimate based on the elementary use of Friedrichs-Poincaré, Cauchy-Schwarz, Young and Korn (in)equalities (see Nečas et al. [54] or Evans [55] or any other standard reference work on function spaces). A precise characterisation of the Reynolds number and the Weissenberg number that guarantee the negativity of the time derivative, and hence the stability, could be obtained by a variational technique known from the standard energy method (see Joseph [56] or Straughan [57]). This is however beyond the scope of the current contribution.
Lemma 3 (Estimate on the time derivative). Let W and W = W + W denote two states of the system governed by equations in Equation (9), where the state W is the steady state, that is it solves Equation (16). Then, there exist constants C 1 ( W, Re, Wi, Ξ, Ω) and C 2 ( W, Re, Wi, Ξ) such that the time derivative of the functional V neq ( W W ) introduced in Equation (32) can be estimated from above as where we denote and where λ min (·) denotes the minimal eigenvalue of the corresponding matrix and C P denotes the domain dependent constant from Friedrichs-Poincaré inequality.
Proof. See Appendix C.
Lemma 4 (Estimate on the time derivative in terms of the norm on the shifted state space). Let W and W = W + W denote two states of the system governed by equations in Equation (9), where the state W is the steady state, that is it solves Equation (16), and let Z and Z denote the corresponding shifted states (see Equation (33)). Let us further assume that the constants C 1 and C 2 in Lemma 3 are negative. Then, there exists a positive constant C(C 1 , where V neq ( W W ) denotes the functional introduced in Equation (32) and · st is the norm introduced in Definition A3.
Proof. In virtue of Lemma 3, we already know the estimate in Equation (40). Making use of Friedrichs-Poincaré inequality v 2 L 2 (Ω) ≤ C P ∇v 2 L 2 (Ω) , where C P is the domain dependent constant. We see that, if C 1 < 0 and C 2 < 0, then Equation (40) implies Next, we use Lemma A4, which implies that where Z κ p(t) denotes the shifted state introduced in Equation (33). Consequently, we see that which means that Equation (43) can be further manipulated to the form The inequality in Equation (46) gives in virtue of the definition of the norm · st (see Definition A3 and the proposition in Equation (42)). (Recall that the transformation in Equation (33) implies that Z κ p(t) = I.)

Main Result
Using the estimate from Lemma 4 and the relation between the metric and the functional V neq (see Lemma 1), it is straightforward to prove the following theorem.
Theorem 1 (Sufficient conditions for unconditional asymptotic stability). Let the pair W = [ v, B κ p(t) ] solve the governing equations for the steady state in Equation (16) with the boundary conditions in Equation (17). If the Reynolds number Re, the Weissenberg number Wi and the dimensionless shear modulus Ξ are such that the constants C 1 and C 2 introduced in Equation (41) satisfy then the spatially inhomogeneous non-equilibrium steady state W is unconditionally asymptotically stable, holds for any initial perturbation W, where the metric dist (·, ·) is the metric introduced in Equation (A7).
Proof. We first investigate the stability in the shifted state space (see Equation (33)), that is we investigate perturbation Z with respect to Z = [ v, I]. We introduce the functional V neq ( W W ) (see Equation (32) and the equivalent expression in Equation (36)). The functional satisfies the condition in Equation (26a) (see Lemma 1). Furthermore, if C 1 and C 2 are negative, then Lemma 4 implies that the functional V neq ( W W ) decreases along the trajectories in a desired manner, that is it satisfies Equation (26b). Consequently, the functional V neq ( W W ) is a genuine Lyapunov type functional for any neighbourhood of the steady state Z, hence Z is unconditionally asymptotically The convergence in the norm on the shifted space seems to be an obscure characterisation of the approach to the equilibrium. However, if we exploit the definition of the shifted state (see Equation (33), and the equality in Equation (A17) proved in Lemma A3), we see that where the last inequality follows from the estimate in Equation (A10) in Lemma A2. Here, E is a positive constant that depends on W, and dist (·, ·) is the metric introduced in Definition A2 (Equation (A7)). This metric is a natural one if we restrict ourselves to the set of positive definite tensor fields. The inequality in Equation (49) then implies Equation (48).
We note that if we want to investigate the spatially homogeneous steady state W = [ v, B κ p(t) ] = def [0, I], that is if we set the boundary condition V = 0, then and the condition in Equation (47)  Having the Lyapunov type functional given by Equation (32), it is interesting to see how the functional works in the case of close to the equilibrium setting, that is for B κ p(t) ≈ I, and for small perturbations, that is for small B κ p(t) . We see that if B κ p(t) is small and if B κ p(t) is close to the identity, then and the proposed Lyapunov type functional V neq can be approximated as V neq ≈ V naive where The functional V naive might be the first candidate for the Lyapunov type functional if the stability is investigated using the popular "energy method" (see, for example, Straughan [57]). The functional is clearly nonnegative, and it vanishes if and only if the perturbation vanishes. Moreover, the candidate V naive for the Lyapunov type functional is much simpler than V neq . Indeed, the proximity of the perturbation to the non-equilibrium state [ v, B κ p(t) ] is measured using the standard Lebesgue space norms, and V naive does not depend on the value of B κ p(t) . Therefore, it seems that V naive is a good candidate for the Lyapunov type functional for the analysis of arbitrary spatially inhomogeneous non-equilibrium steady state [ v, B κ p(t) ].
The inappropriateness of V naive for the stability analysis is in fact apparent even in a very trivial setting. Let us consider the spatially homogeneous equilibrium steady state B κ p(t) = I, v = 0 in a mechanically isolated container, that is we set V = 0 in the boundary condition in Equation (12). If we use the (exact) evolution equations for the perturbation velocity in Equation (A28b), and if we evaluate the time derivative of V naive , then we get (see also Equation (A34)). The last term on the right-hand side of Equation (53) can be evaluated using the (exact) evolution equation for B κ p(t) (see Equation (A29)). Substituting Equation (A29) into Equation (53) and using the fact that B κ p(t) = I and v = 0, yields Using the standard manipulation we see that the second term on the right-hand side of Equation (54) dv. (56) Let us now consider an initial perturbation that is chosen is such a way that Ω Ξ Tr( D B κ p(t)

2
) dv > 0, which can certainly be done. This positive value will dominate the right-hand side of Equation (56) provided that the Reynolds number and the Weissenberg number are large enough. Consequently, V naive will (initially) increase, and it would be useless as the Lyapunov type functional unless we a priori limit ourselves to small perturbations.
On the other hand, if we use the functional V neq in the case B κ p(t) = I and v = 0, then we immediately see that the constants C 1 and C 2 in Equation (40) are negative, and that the equilibrium steady state is asymptotically stable with respect to any perturbation and any value of the Reynolds and the Weissenberg number! (Note that Guillopé and Saut [58] obtained only a conditional stability result in a Sobolev space norm for the equilibrium rest state B κ p(t) = I and v = 0 (see their Corollary 3.5 and assumptions of Theorem 3.3).) Based on the analysis presented above, we can therefore claim that we have indeed benefited from a well constructed Lyapunov type functional V neq and the choice of metric. Unlike the naive Lyapunov type functional V naive , the proposed Lyapunov type functional V neq seems to properly reflect the nonlinearity of the governing equations and the related energy storage mechanisms and the entropy production mechanisms.

Taylor-Couette Flow
Let us now consider a viscoelastic fluid described by the Giesekus model introduced in Section 2 with α = 1 2 , and let us investigate the stability of steady flow in the standard Taylor-Couette flow geometry (see Figure 1). The objective is to show as how the theory introduced above works in a specific setting. The choice α = 1 2 is motivated by the simplicity of the expressions for the corresponding steady state.
The fluid is placed in between two infinite concentric cylinders of radii R 1 and R 2 , with R 1 < R 2 . The cylinders are rotating with the angular velocities Ω 1 (inner cylinder) and Ω 2 (outer cylinder) along the common axis. The geometry naturally leads to the use of cylindrical coordinates (r, ϕ, z); the normed basis vectors are denoted as gr, gφ and gẑ (see Figure 1). Since the domain is unbounded in the z-direction, we henceforth consider a periodic cell where h > 0 is arbitrary, and we use the notation Γ 1 = def {(r, ϕ, z) ∈ R 3 | R 1 < r < R 2 , 0 ≤ ϕ < 2π, |z| = h} for the top and bottom base, and Γ 2 = def {(r, ϕ, z) ∈ R 3 | r ∈ {R 1 , R 2 }, 0 ≤ ϕ < 2π, |z| < h} for the cylindrical walls of the domain. The flow is driven by the rotation of the cylinders along the common axis.

Base Flow-Non-Equilibrium Steady State
The first task in the stability analysis is to find the steady solution to the governing equations. This solution is the spatially inhomogeneous non-equilibrium steady state W as introduced in Section 3.2. The characteristic length and characteristic time have been chosen as x char = def R 1 , t char = def 1 Ω 1 . We use the periodic boundary condition on Γ 1 and the no-slip boundary condition for velocity field v on Γ 2 , that is v| r=R 1 = R 1 Ω 1 gφ, v| r=R 2 = R 2 Ω 2 gφ. These boundary conditions are consistent with the requirements on boundary conditions specified in Section 2.4. In their dimensionless form, the boundary conditions read where we have introduced two dimensionless parameters η = def R 1 R 2 and ζ = def Hereafter, we work exclusively with the dimensionless variables and thus, for the sake of simplicity, we omit the star denoting them.
Since the problem has the rotational symmetry, we search for the steady non-equilibrium state in a special form. Namely, the solution to Equation (16) subject to boundary conditions in Equation (58) is sought in the form Note that the chosen ansatz for the velocity field automatically satisfies the incompressibility condition. The assumptions lead to the following expressions for the velocity gradient, the symmetric part of the velocity gradient, the convective term, the divergence of B κ p(t) , and the upper convected derivative of B κ p(t) : where we introduce the angular velocity ω(r), vφ(r) = def ω(r)r. Using Equation (60), we see that the governing equations for the velocity field in Equation (16b) reduce to Assuming that dω dr = 0 in (R 1 , R 2 ), Equation (61b) can be solved for Brr, Brφ, Bφφ and Bẑẑ. However, for general α ∈ (0, 1), the formulae for the aforementioned quantities are too complex to be written down here. Let us note however that for α = 1 2 the formulae simplify significantly; the solution to Equation (61b) which satisfies the condition of B κ p(t) being positive definite in this case reads where we denote c = def 2Wi r dω dr . Substituting Equation (62) into the second equation in Equation (61a) then yields an ordinary differential equation for the angular velocity ω supplemented by the boundary conditions ω| r=1 = 1, ω| r= 1 η = ζ, which follow from Equation (58) and the fact that vφ(r) = ω(r)r. Equation (63) together with the boundary conditions constitute a boundary value problem which needs to be solved numerically.

Explicit Criterion for the Stability of Spatially Inhomogeneous Non-Equilibrium Steady State
Here, we explicitly compute constants C 1 , C 2 defined by Equations (41a) and (41b) for the Taylor-Couette problem and for the specific values of the dimensionless numbers Ξ, Re and Wi. Let us recall that, for the sake of simplicity, we set α = 1 2 , and we consider the steady tensor field B κ p(t) given by Equation (62). We fix the values for the geometric parameter η and angular velocities ratio ζ as η = 1 2 and ζ = 2. The angular velocity ω is obtained by solving Equation (63) which is a boundary-value problem for a second order nonlinear differential equation. The problem was solved numerically using the function solve_bvp from SciPy library version 1.0.0, which implements a fourth-order collocation algorithm with the control of residuals as described in Kierzenka and Shampine [59]. With the angular velocity ω in hand, we immediately get the steady velocity field v = ω(r)rgφ, and the steady left Cauchy-Green tensor field B κ p(t) through Equation (62). The plots of the velocity field and the components of B κ p(t) are shown in Figure 2.    Having computed the steady velocity field v and the corresponding steady field B κ p(t) , we can evaluate the constants C 1 and C 2 in the estimate in Equation (40). The gradient of v as well as the gradient of B κ p(t) are again computed numerically from the obtained numerical solution. Finally, the Poincaré constant for the cylindrical annulus is determined via an explicit solution of the corresponding eigenvalue problem −∆u = λu for the Laplace operator with Dirichlet boundary condition, which leads, for the geometrical parameter η = 1 2 , to the value C P ≈ 0.1025. The resulting stability regions in the Re-Wi plane are shown in Figure 3 for a fixed value of the dimensionless shear modulus Ξ. As one might expect the spatially inhomogeneous steady state is indeed unconditionally asymptotically stable if the Weissenberg number Wi and the Reynolds number Re are small enough. Re Re  Unconditional asymptotic stability is granted provided that C 1 < 0 and C 2 < 0, numerical values of constants C 1 and C 2 are evaluated using Equation (41). Giesekus parameter α = 1 2 and problem parameters η = 1 2 and ζ = 2.

Numerical Experiments-Evolution of Various Initial Perturbations
Finally, we document the theoretically predicted behaviour by numerical experiments. The numerical experiments allow us to quantitatively track the evolution of key quantities such as the net kinetic energy, and also to quantitatively monitor the energy exchange between the fluid and its surroundings.
The governing equations were numerically solved using standard techniques. The weak forms of the governing equations were discretised in the space using the finite element method, while the time derivatives were approximated with the backward Euler method. The two-dimensional domain Ω was discretised by regular quadrilaterals. The mesh divided the annular region Ω into 80 pieces in the radial direction, and in 720 pieces in the azimuthal direction. The corresponding total number of degrees of freedom in all numerical experiments was over 1.3 × 10 6 . The velocity field v and the B κ p(t) field were approximated by biquadratic Q2 elements, and the pressure field m was approximated by the piecewise linear discontinuous elements P1 d (see Korelc and Wriggers [60] for details). The finite element pair that was used for the velocity/pressure fields satisfied the Babuška-Brezzi condition, the finite element for B κ p(t) field was chosen to be the same as for the velocity in order to provide rich enough finite element space for the solution. The same finite elements were chosen for the two-dimensional simulation of other viscoelastic rate-type fluids (Oldroyd, Burgers and their various nonlinear versions) by Hron et al. [61] and Málek et al. [62]. In the three-dimensional case, low order elements can be used to decrease the overall cost of the calculation (see Tůma et al. [63]).
The numerical scheme was implemented in the AceGen/AceFEM system (see Korelc [64] and Korelc [65]). The main advantage of the system is that it provides automatic differentiation used for the computation of the exact tangent matrix needed by the Newton solver that treats all nonlinearities. The final set of linear equations was solved by the direct solver Intel MKL Pardiso. The stopping criterion for the Newton solver was set to 10 −9 .
Using the numerical scheme, we are ready to study the behaviour of various perturbations to the non-equilibrium steady state. In all scenarios described below, we use the dimensionless parameters and we fix the geometric parameter η and angular velocities ratio ζ as in Section 6.2, that is η = 1 2 and ζ = 2. The chosen values of η, ζ and Ξ correspond to the stability diagram shown in Figure 3a.
The values of Reynolds number and Weissenberg number are outside the region where we have proven the decay of the proposed Lyapunov type functional. Nevertheless, as we show below, the Lyapunov type functional is, in the cases being investigated below, still a decreasing function. First, we start from the homogeneous steady state solution [v, B κ p(t) , m] = [0, I, 0], and we let the system to spontaneously evolve up to the time instant t = 1000. (More precisely, the initial condition is v = 0 inside the domain Ω, and Equation (17) holds on the boundary of Ω. After the first computational time step, which is chosen as ∆t = 0.05, we get on the discrete level a divergence-free velocity field with the appropriate boundary condition. This discrete velocity field provides us a consistent initial condition for further computations. Therefore, we formally start the evolution not at t = 0, but at t = 0.05.) At this time instant, the system is almost relaxed and is close to the steady solution. The solution at t = 1000 is used as a starting point for solving the steady governing equations (without the time derivatives) and the spatially inhomogeneous non-equilibrium steady state is obtained just in two Newton iterations. (The finite element solution coincides with the semi-analytical steady solution obtained in Section 6.1. This among others provides us a tool for the code verification.) Consequently, the finite element solution is in what follows used as the spatially inhomogeneous non-equilibrium steady state W.
Having obtained the numerical representation of the spatially inhomogeneous non-equilibrium steady state, we proceed with two scenarios concerning the specification of the initial perturbation.

Scenario A-Localised Perturbation of the Left Cauchy-Green Field
In the first scenario, we keep the initial velocity field perturbation equal to zero, while the initial perturbation in B κ p(t) is localised in space (see the first snapshot in Figure 4). Since the system is fully coupled, the perturbation in the B κ p(t) field triggers for t > 0 a nontrivial evolution of the velocity perturbationṽ (see Figure 5). This can be observed also in the plots showing the evolution of the net elastic stored energy and the net kinetic energy (see Figure 6).
Bκ p(t) Figure 4. Scenario A, snapshots of | B κ p(t) | at different time instants. Finally, we also investigate the time evolution of the proposed Lyapunov type functional V neq and the naive Lyapunov type functional V naive , and the net mechanical energy flux going through the boundary of Ω (see Figure 6). Although we work with the Reynolds number/Weissenberg number pair outside the guaranteed stability region, we see that the value of Lyapunov type functional V neq still decreases in time, and that the perturbation vanishes for t → +∞. This indicates that the estimates on the time derivative of the proposed Lyapunov type functional are, at least for a class of perturbations, too strict and they might be improved. One should also note that the "net kinetic energy" of the perturbation, that is the functional Ω 1 2 | v| 2 dv, does not decrease for all t > 0 (see Figure 6b). In fact, it experiences a transitional growth, and such a transient growth can be observed even for the Reynolds number/Weissenberg number values within the stability region. This is a natural observation. The elastic energy stored in the material can be released in the form of the kinetic energy. It is only the combination of the elastic energy and the kinetic energy that appears in the Lyapunov type functional that leads to a quantity that decays at any time.
Further, the net mechanical energy flux through the boundary fluctuates around the value that corresponds to the non-equilibrium steady state, and then it reaches the value that corresponds to the spatially inhomogeneous non-equilibrium steady state (see Figure 6d). This can again happen even if the Reynolds number/Weissenberg number take values within the stability region.

Scenario B-Global Perturbation of the Velocity Field
In the second scenario, we start with a nonzero velocity perturbation v, and the B κ p(t) field is kept unchanged, The initial velocity v is prescribed as v| t=0 = Ωrgφ, where the angular velocity is the arithmetic mean of the two angular velocities Ω = Ω 1 +Ω 2 2 . (Formally, we apply the same procedure as discussed in the previous section. The initial condition is v = Ωrgφ inside the domain Ω, and Equation (17) holds on the boundary of Ω. The actual computation starts after the first (formal) time step, when the discrete velocity field is divergence-free and it fulfills the boundary condition.) Again, as in the previous case, the non-zero perturbation in one unknown field ( v) triggers for t > 0 a nontrivial evolution of the other unknown field ( B κ p(t) ) (see Figures 7 and 8).
Bκ p(t) Figure 7. Scenario B, snapshots of | B κ p(t) | at different time instants. Moreover, this numerical experiment is instructive for yet another reason. In Figure 9c, we plot the time evolution of the values of the functionals V neq and V naive . Clearly, the functional V naive (see Equation (52)), which is a naive candidate for the Lyapunov type functional, experiences a transitional growth. Interestingly, the proposed complex Lyapunov type functional V neq is still a decreasing function, although the Reynolds number/Weissenberg number values are outside the region, where we have actually proven the decay of the functional. This further indicates that the functional V naive is indeed not a good candidate for a Lyapunov type functional (see also Section 5 for further discussion).

Conclusions
We have investigated the stability of spatially inhomogeneous non-equilibrium steady states (flows) of viscoelastic fluids described by the Giesekus model. We have derived bounds on the values of the Reynolds number and the Weissenberg number that guarantee the flow stability subject to any finite perturbation. The stability has been investigated using the Lyapunov type functional given by the formula A few observations concerning the proposed Lyapunov type functional are at hand. First, the proposed Lyapunov type functional has a relatively complicated form. In particular, it is not quadratic in the perturbation B κ p(t) , and it depends on the spatially inhomogeneous non-equilibrium state B κ p(t) . This makes it remarkably different from a naive Lyapunov type functional of the form which might be a first try if the stability problem were analysed using the popular "energy method". However, as we have shown, the complicated structure of the proposed functional V neq leads to a relatively simple and well structured formula for its time derivative, which in turn allows one to formulate conditions that guarantee the negativity of the time derivative. Furthermore, the complicated structure of the proposed functional V neq also leads to a simple relation between the functional and the metric on the set of spatially distributed symmetric positive definite matrices. Second, the Lyapunov type functional has been used in the investigation of stability of solution to the complete system of nonlinear governing equations. In particular, the evolution equations for the perturbation have been investigated without any simplification. This makes the present approach different from the "energy budget" analysis (see, for example, Joo and Shaqfeh [66], Byars et al. [67], Ganpule and Khomami [68], Smith et al. [69], Karapetsas and Tsamopoulos [70], Pettas et al. [71], and especially the work by Grillet et al. [19] who investigated the Giesekus model). The "energy budget" analysis, although valuable in the discussion of the nature of the instability mechanisms, is based on the linearised momentum equation for the perturbation and linearised constitutive equation for the "polymeric stress". Consequently, the standard "energy budget" analysis is, unlike the present approach, a tool that cannot be used in the finite amplitude stability analysis of the complete system of nonlinear governing equations. One might also note that, despite the complexity of the proposed Lyapunov type functional, the formula for its time derivative is in fact quite simple compared to the formulae in the "energy budget" analysis. This happens even though the "energy budget" formulae paradoxically stem from various simplifications of the original system of governing equations.
Third, the Lyapunov type functional has been designed using thermodynamical arguments. In fact, the proposed Lyapunov type functional has been constructed from the net mechanical energy functional E mech (see Equation (10), via Equation (30)). This makes the construction quite general, and one might speculate that a similar approach is very likely applicable to other popular viscoelastic rate-type models such as the PTT model (see Phan Thien and Tanner [72]) or the FENE-P model (see Bird et al. [73]), as well as complex viscoelastic rate-type models with, for example, stress diffusion terms (see Málek et al. [37] or Dostalík et al. [39]). Further, the construction of the Lyapunov type functional has been based on the method proposed by Bulíček et al. [23], and this method is speculated to work even for complex coupled thermo-mechanical systems. This naturally calls for the investigation of the applicability of the method in more complex settings such as flows of viscoelastic rate-type fluids with temperature dependent material parameters.
Fourth, thermodynamical type considerations such as the identification of the energy storage mechanisms and entropy producing mechanisms are known to play an important role in the rigorous mathematical theory of nonlinear partial differential equations governing the motion of viscoelastic fluids (see, for example, Hu and Lelièvre [74], Boyaval et al. [40], Barrett and Boyaval [75], Barrett and Boyaval [76], Barrett and Süli [47] or Bulíček et al. [77]). On the other hand, rigorous mathematical analysis of long-time behaviour of viscoelastic fluids is usually done without a direct appeal to thermodynamics, and the available results are quite limited especially if one considers thermodynamically open systems (see, for example, Guillopé and Saut [78], Nohel and Pego [79], Jourdain et al. [80] or Renardy [81]). (Usually, only stability of unidirectional steady flows in simple geometries is considered.) Consequently, the approach proposed in the current contribution might be of interest from the rigorous mathematical perspective as well. This means that one should deal with the weak solution to the governing equations, and that one should investigate the applicability of the presented arguments for a solution/perturbation that has only the smoothness that can be actually proven.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Distance between Positive Definite Matrices and Its Generalisation to Spatially Distributed Tensor Fields
The function dist P(d), BW (A, B) defined as and referred to as the Bures-Wasserstein distance, is known to be a metric on the manifold of symmetric (Hermitian) positive semidefinite matrices P(d) ⊂ R d×d of arbitrary dimension d (see Bhatia et al. [24]).
(See also the work of Bhatia [82] for a thorough discussion of properties of symmetric positive semidefinite matrices.) Another possibility regarding the definition of a metric on the manifold of symmetric positive definite matrices is (see Bhatia [82]), where |·| denotes the Frobenius norm, |M| = def (Tr (MM )) 1 2 . (See the work of Graham [83] for comments on the use of Equation (A2) in the context of fluid mechanics. For further discussion regarding the metric on the set P(d), see also Hiai and Petz [84] and Hiai and Petz [85].) Using the metric in Equation (A1) or Equation (A2), one can introduce a metric on the set of spatially distributed symmetric positive semidefinite matrices. where we use the triangle inequality for the Bures-Wasserstein/δ 2 metric dist P(d), · (A, B). Using the Hölder inequality in both integrals on the right-hand side, we arrive at which, upon dividing both sides by dist P Ω (d), · (A, B), yields the desired result. Now, we are in position to define a metric that will allow us to characterise the distance between the states Y = [v, B κ p(t) ] and Y = [v, B κ p(t) ] of the dynamical system given by equations in Equation (3). Note that in Definition A2 we do not assume that some of the states Y or Y is a steady state.
Definition A2. Let Y and Y be the states of the dynamical system governed by equations in Equation (9). The function defines a metric on the set of states. Here, dist P Ω (d), BW (·, ·) denotes the metric on the set of spatially distributed symmetric positive semidefinite matrices introduced in Definition A1, and · L 2 (Ω) denotes the standard norm on the Lebesgue space L 2 (Ω).
Proof. It is straightforward to verify that Equation (A7) indeed defines a metric. The three fundamental properties of the metric-nonnegativity, symmetry and identity of indiscernibles-are clearly satisfied in virtue of the properties of the L 2 (Ω) norm and the metric dist P Ω (d), BW (·, ·). The triangle inequality follows from the discrete Minkowski inequality p , (see, for example, Evans [55]), and the triangle inequalities for the metric induced by the L 2 (Ω) norm and the metric dist P Ω (d), BW (·, ·).
Besides the metric in Equation (A7), we also make use of the conventional norm/metric introduced via Equation (A8). (The nomenclature "shifted state space" is explained in Section 4.2.) The fact that Equation (A8) defines a norm is straightforward to show using the same argument as in the proof of correctness of Definition A2.
Definition A3 (Norm/metric on the shifted state space). Let Z = def v, Z κ p(t) and Z = def v, Z κ p(t) be the states of the dynamical system governed by the equations in Equation (3) with the transformation introduced in Equation (33). The function defines a norm on the corresponding shifted state space, and the function is the metric induced by the norm in Equation (A8).
Lemma A2 (Pointwise inequality for Bures-Wasserstein metric I). Let A, B ∈ P(d), and let A be an invertible matrix. If |·| denotes the standard Frobenius norm, then Proof. Bhatia et al. [24] (in their Theorem 1) showed that the Bures-Wasserstein metric can be equivalently defined as where U is a matrix that belongs to the group of unitary matrices U(d) of dimension d × d. Making use of the submultiplicativity of the Frobenius norm, we see that Following the argument given in Bhatia et al. [24], one can show that the maximum value of the term Tr [U X + X U] is attained for the unitary matrix U that corresponds to the unitary matrix V in the polar decomposition of X, X = VP. The matrix V in the polar decomposition of X can be found explicitly as V = X (X X) − 1 2 , which yields If we set U = def V in the last term on the right-hand side of Equation (A12), then we get max U∈U(d) Moreover, if |·| denotes the Frobenius norm, then Proof. The proof follows from the spectral decomposition. Matrix A − 1 2 BA − 1 2 is a symmetric positive definite matrix, hence it is unitarily similar to a diagonal matrix with the eigenvalues of A − 1 2 BA − 1 at the diagonal. Consequently, if {λ i } d i=1 are eigenvalues of A − 1 2 BA − 1 2 , then the right-hand side of Equation (A16) reads while on the left-hand side of Equation (A16) we use the definition of the metric (see Equation (A1)), and we get Instead of Equation (A16), we therefore need to prove the following inequality for the eigenvalues which reduces to However, if λ i > 0, which is our case since we are dealing with positive definite matrices, then it is straightforward to check that Indeed, the function on the left-hand side of Equation (A22) is a function that vanishes for λ i = 1, while its derivative reads − 1 , which implies that the function is increasing for 0 < λ i < 1, and it is decreasing for λ i > 1. The point λ i = 1 is therefore the point where 2(1 − √ λ i ) + ln λ i reaches its global maximum, and Equation (A22) holds. Equation (A21) and consequently also Equation (A20) then follows via summation of Equation (A22) over the individual eigenvalues. The equality in Equation (A17) is a straightforward consequence of the definition of the Frobenius norm, and the definition of the Bures-Wasserstein metric (see Equation (A15)).
Lemma A4 (Inequality for the Frobenius norm). Let A, B ∈ P(d), let A and B be invertible matrices, and let |·| denote the Frobenius norm. Then, where we have used the submultiplicativity of the Frobenius norm, |XY| ≤ |X| |Y|. Next, we use the spectral decomposition of the symmetric positive definite matrix A − 1 2 BA − 1 2 , and we find that

Appendix B. Formula for the Time Derivative
Let us now derive Equation (38) for the time derivative of the proposed Lyapunov type functional. Straightforward differentiation of Equation (32) under the integral sign yields dVneq dt In the following, we treat the individual terms of Equation (A26) separately. To proceed further, we need to formulate the evolution equations for the perturbation W, which gives us the formulae for the partial time derivatives of v and B κ p(t) .

Appendix B.1. Evolution Equations for Perturbation
The perturbed field W = W + W must satisfy the governing equations in Equation (9), which means that and and where we have used the notation L = ∇ v, L = ∇ v and similarly for the symmetric part of the velocity gradient D and D. Now, we are in position to exploit the fact that the non-equilibrium steady state W solves Equation (16). Using Equation (16) in Equation (A27) yields where and which can be in virtue of Equation (16c) further simplified to In the subsequent analysis, it is however more convenient to work with Equation (A28d) instead of Equation (A29). Equation (A29) is exploited only in Section 5. The system in Equation (A28) of evolution equations for the perturbation W must be solved subject to the boundary condition in Equation (13) and periodic boundary condition on Γ 1 .
Having identified the formulae for the time derivatives of v and B κ p(t) , we can go back to Equation (A26), and we can start to evaluate the individual terms on the right-hand side of Equation (A26) for the time derivative of the Lyapunov type functional.

Appendix B.3. Second Term of Equation (A26)
Using the evolution in Equation (A28d) for B κ p(t) yields We can immediately see that the second and the third terms vanish due to the incompressibility condition in Equation (A27a) and the invariance of the trace under cyclic permutations. The first term on the right-hand side of Equation (A35) can be shown to vanish as well via the standard manipulation (The last equality again follows from the Stokes theorem, the identity in Equation (14) and the incompressibility condition in Equation (A27a). Moreover, we have also used the fact that u • (∇ ln det A) = Tr A −1 (u • ∇) A .) Finally, we see that Multiplying from the left, taking the trace and integrating over the domain Ω yields Consequently, using the invariance of the trace under cyclic permutations and rearranging the terms we obtain the identity Having the identity, let us now go back to Equation (A26), and let us manipulate the third term on the right-hand side of Equation (A26). Employing the evolution equation in Equation (A28d) for B κ p(t) yields The first term on the right-hand side of Equation (A40) reduces to where we use a similar manipulation as in Equation (A36). Moreover, the expression can be further transformed as follows Thus far, we have found that which-upon exploiting the identity in Equation (A39) in the first term-reduces to where we again use the incompressibility condition.
dv, (A45) and using the resolvent identity A −1 we can rearrange the next to the last term in Equation (A45), and we see that the time derivative of the proposed Lyapunov type functional is indeed Equation (38).
where u is a (smooth) vector field that vanishes on Γ 2 and satisfies the periodic boundary condition on Γ 1 , and D u denotes the symmetric part of the corresponding gradient ∇u. The Friedrichs-Poincaré inequality reads u 2 L 2 (Ω) ≤ C P ∇u 2 L 2 (Ω) , where C P is the domain dependent constant, u is a (smooth) vector field that vanishes on Γ 2 and satisfies the periodic boundary condition on Γ 1 , and w 2 L 2 (Ω) = def Ω |w| 2 dv denotes the standard Lebesgue space norm, and |w| denotes the standard Euclidean norm. Now, we are in position to find an estimate on the right-hand side of Equation (38). In estimating the time derivative in Equation (38), we can completely ignore the next-to-the last nonpositive term on the right-hand side of Equation (38). We can thus write This means that we lose a term that has the negative sign, and that the estimate on the stability range will be more demanding than it would have to be. We also need to restrict ourselves to α = 0.
Let us now bound the individual terms involved in Equation (A47). The first term can be in virtue of Korn equality rewritten as The third term on the right-hand side of Equation (A47) can be estimated using the spectral estimate 2 for the symmetric matrix D, where λ min ( D) denotes the smallest eigenvalue of D and λ max ( D) denotes the largest eigenvalue of D at the given spatial point x.
The spectral estimate yields Further, using the Poincaré inequality, we get The last term in Equation (A47) can be estimated as where λ min ( B κ p(t) ) denotes the smallest eigenvalue of the given symmetric positive definite matrix B κ p(t) at the given spatial point x. (Note that λ min ( B κ p(t) ) is in virtue of the positivity of B κ p(t) a positive number.) Estimates on the rest of the terms are obtained easily by the application of Cauchy-Schwarz and Young inequalities and the submultiplicative property of the matrix norm. First, we group the second and the fifth term in Equation (A47), and we get The fourth term in Equation (A47) is in virtue of the Poincaré inequality estimated as where we introduce the notation ∇ B κ p(t) 2 = def ∑ 3 k,l,m=1 give us Equation (40).