\section{Introduction}
\label{sec:introduction}

Consider two concentric cylinders, the outer rotating with respect to the inner one. The gap between the cylinders is filled with a non-Newtonian fluid. In such an experiment either the force required to rotate the outer cylinder (shear stress) or speed of the rotation (shear rate) can be controlled, the other being measured, see Fig.~\ref{fig.exp}.\footnote{We assume that the measured shear stress and shear rate corresponds to the actual properties of the fluid within the gap.} An interesting experimental observation is that when varying shear stress, shear rate behaves continuously, while when varying shear rate, shear stress exhibits a jump \cite{boltenhagen.p.hu.y.ea:observation}.

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.6\textwidth]{boltenhagen-hu-fit-arrow.pdf}
  \caption{Steady-state stress/shear--rate behavior for a 7.5/7.5 mM TTAA/NaSal solution from the controlled shear stress ($\Delta$) and controlled shear--rate ($\bullet$) experiments by Boltenhagen et al.~\cite{boltenhagen.p.hu.y.ea:observation}. The fit (-\,-) was obtained from the constant applied stress data using one-dimensional form of the constitutive relation~\eqref{eq:le-roux-rajagopal}. \textbf{When increasing shear stress, shear rate varies continuously} (although not monotonically). \textbf{When increasing shear rate, shear stress exhibits a jump.} The dashed arrow indicates the experimentally observed jump in shear stress while the dotted arrow the predicted (see Sec.~\ref{sec.stability}) jump when lowering shear rates (hysteresis). Note that the vertical axis expresses the strain rate, in contrast to the usual habit, and the S-curve in this figure is thus referred to as vertical.}
  \label{fig.exp}
\end{figure}

We ask the following two questions:
\begin{enumerate}
\item How to describe the relation between shear stress and shear rate by means of a constitutive relation compatible with non-equilibrium thermodynamics?
\item What is the reason that when controlling the shear stress, the response (shear rate) is continuous, while when controlling the shear rate, the response (shear stress) is discontinuous?
\end{enumerate}

The first question is answered in Sec.~\ref{sec.CR}, where an appropriate model is fitted to the experimental data and then reformulated in terms of gradient dynamics. The second question is partially answered in Sec.~\ref{sec.conjugate}, where the multivalued conjugate gradient structure is identified. Subsequently, in Sec.~\ref{sec.CR-thermodynamics}, the dynamics is lifted to an extended state space. This enables us to discuss various regimes of the dynamics and its stability, metastability and instability with respect to perturbations of the constitutive relation---the CR-stability, and to complete answering the second question in Sec.~\ref{sec.stability}, which is the essential result of this paper. Moreover, we propose that the phenomenon of critical heat flux in Sec. \ref{sec.CHF} could be approached by the method of non-convex dissipation potentials.

Novelty of this paper lies in the following points. A dissipation potential is identified that yields constitutive relation~\eqref{eq:le-roux-rajagopal}. The dissipation potential is not convex everywhere (it is convex in the vicinity of equilibrium), and CR-stability of the dynamics generated by this dissipation potential (with respect to perturbations of the dissipation potential) is discussed in the context of mesoscopic multiscale thermodynamics. The conjugate dissipation potential, which is obtained by a form of Legendre transformation, is multivalued, which, to our best knowledge, is mentioned in the thermodynamic literature for the first time whatsoever. Moreover, a thermodynamic reason is identified which favors seeking constitutive relation in the form $\ensuremath{\tensorq{D}}({\traceless{\cstress}})$, i.e. shear rate as a function of shear stress. Finally, the experimental behavior presented in Fig.~\ref{fig.exp} is explained. A model based on a non-convex dissipation potential is also proposed for the critical heat flux in Sec \ref{sec.CHF}.

\section{Constitutive relation}
\label{sec.CR}

The experimental dependence of shear rate on shear stress plotted in Fig.~\ref{fig.exp} can be described by a model of non-Newtonian fluid introduced by Le~Roux and Rajagopal~\cite{roux.c.rajagopal.kr:shear}
\begin{equation}
  \label{eq:le-roux-rajagopal}
  \ensuremath{\tensorq{D}} = \left[ \alpha \left( 1 + \beta \absnorm{{\traceless{\cstress}}}^2 \right)^n + \gamma \right] {\traceless{\cstress}},
\end{equation}
where $\alpha$ and $\beta$ are positive numbers, $\gamma$ is a non-negative number, $n$ is a real number, $\ensuremath{\tensorq{D}} =_\mathrm{def} \frac{1}{2} \left( \ensuremath{\nabla \vecv} + \transpose{\ensuremath{\nabla \vecv}} \right)$ is the symmetric part of the velocity gradient and ${\traceless{\cstress}} =_\mathrm{def} \cstress - \left( \frac{1}{3}\Tr \cstress \right) \identity$ is the deviatoric part of the Cauchy stress ${\tensorq{T}}$. Setting $\gamma = 0$, constitutive relation~\eqref{eq:le-roux-rajagopal} reduces to the stress power-law model studied by Málek et al.~\cite{malek.j.prusa.ea:generalizations}.  The fit is plotted in Fig.~\ref{fig.exp}.

Note that only the deviatoric part of the stress tensor is taken into account and that the fluid is considered to be effectively incompressible. It is possible that abandoning such a constraint could lead to other interesting behavior, but the purpose of this paper is achieved while respecting this constraint.

Constitutive relation \eqref{eq:le-roux-rajagopal} can be seen as a special case of general implicit constitutive relation \cite{rajagopal.kr:on*3,rajagopal.kr:on*4}
\begin{equation}
  \label{eq:general-implicit-relation}
  \tensorf{h} ({\traceless{\cstress}}, \ensuremath{\tensorq{D}}) = \tensorq{0},
\end{equation}
where $\tensorf{h}$ is a tensorial function. It is however not clear whether the fully implicit case brings any advantage and whether it is physically substantiated. On the other hand, the special case $\ensuremath{\tensorq{D}}({\traceless{\cstress}})$ has a sound thermodynamic meaning and practical advantages, see Sec.~\ref{sec.CR-thermodynamics} and \cite{perlacova.t.prusa.v:tensorial}.

Constitutive relation \eqref{eq:le-roux-rajagopal} can be also regarded as a consequence of gradient dynamics\footnote{Gradient dynamics is a cornerstone of the GENERIC framework \cite{GO,OG}, where it generates the dissipative evolution. It has strong connection to the principle of large deviations \cite{Mielke-Peletier}, and it can be related to other techniques associated with the maximization of entropy production \cite{SEA,Rajagopal2004}, where for instance quantum dissipation is addressed. However, we prefer using gradient dynamics, since it seems to be more straightforward in the case of non-quadratic dissipation potentials \cite{janecka.a.pavelka.m:gradient}. Note that an alternative non-potential version of GENERIC is preferred by others \cite{HCO,Hutter2013}. Although the latter version is more general as it allows for antisymmetric dissipative coupling, it is not clear whether such a coupling is necessary and correct \cite{Miroslav-WhyGeneric}. We prefer the potential version of GENERIC because reversibility and irreversibility are then clearly distinguished \cite{PRE2014} and because we can then perform Legendre transformations, which are crucial for the manuscript at hand.}
\begin{subequations}
  \begin{equation}
    \label{eq.gradient}
    \ensuremath{\tensorq{D}} = \pd{\Xi^*}{{\traceless{\cstress}}},
  \end{equation}
  where the shear rate is a derivative of the dissipation potential $\Xi^*$ with respect to the shear stress.\footnote{The derivative is to be interpreted as the functional derivative in general. It will be also denoted simply by a subscript, i.e. $\Xi^*_{{\traceless{\cstress}}}$.} Alternatively, one can use the conjugate form
  \begin{equation}
    \label{eq.cgradient}
    {\traceless{\cstress}} = \pd{\Xi}{\ensuremath{\tensorq{D}}},
  \end{equation}
\end{subequations}
where the dissipation potential $\Xi$ depends on the shear rate, and it is conjugate to $\Xi^*$ by means of the Legendre transformation, see Sec.~\ref{sec.conjugate} for more details. The former case, Eq.~\eqref{eq.gradient}, is covered by dissipation potential
\begin{equation}
  \label{eq.theta}
  \Xi^* \left( {\traceless{\cstress}} \right) =
  \begin{cases}
    \frac{\alpha}{2 \beta (n + 1)} \left[ (1 + \beta \absnorm{{\traceless{\cstress}}}^2)^{n+1} - 1 \right] + \frac{\gamma}{2} \absnorm{{\traceless{\cstress}}}^2, & n \neq -1, \\
    \frac{\alpha}{2 \beta} \ln \left(1 + \beta \absnorm{{\traceless{\cstress}}}^2 \right) + \frac{\gamma}{2} \absnorm{{\traceless{\cstress}}}^2, & n = -1,
  \end{cases}
\end{equation}
whose derivative gives Eq.~\eqref{eq:le-roux-rajagopal} in the sense of Eq.~\eqref{eq.gradient}. This dissipation potential is plotted in Fig.~\ref{fig:dissipation-potential}. The latter case, Eq.~\eqref{eq.cgradient}, is covered by the dissipation potential $\Xi$ specified in Sec.~\ref{sec.conjugate}.

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.6\textwidth]{dissipation-potential.pdf}
  \caption{Dissipation potential $\Xi^{*} \left( {\traceless{\cstress}} \right)$ corresponding to the constitutive relation~\eqref{eq:le-roux-rajagopal} with parameter values $\alpha = 1$, $\beta = 1$, $\gamma = 0.02$ and $n = -2$. The segment $C-D$ is not convex, thereby the dissipation potential is not CR-stable in this region. The dashed tangent line $B-E$ indicates the convex hull of the potential and coincides with the Maxwell lever rule construction. The segments $B-C$ and $D-E$ correspond to the CR-metastable regions (compare with Fig.~\ref{fig:flow-curve}).}
  \label{fig:dissipation-potential}
\end{figure}

The important feature of the dissipation potential $\Xi^*$, Eq.~\eqref{eq.theta}, is that it is convex in the vicinity of zero while non-convex for higher shear stresses. Convexity is then regained for even higher shear stresses. By analogy with equilibrium thermodynamics \cite{Callen}, it can be expected that the convex parts generate CR-stable or CR-metastable evolution while the non-convex parts lead to CR-instability. This is indeed so as discussed in Sec.~\ref{sec.CR-thermodynamics} and Sec.~\ref{sec.stability}.

\subsection{Second law of thermodynamics}
\textbf{On the second law of thermodynamics and convexity}: Assume for a moment that the dissipation potential $\Xi(\gradsym)$ is smooth and convex. The potential is zero at the origin an thus is positive everywhere except in the thermodynamic equilibrium. Due to the convexity, the origin is also the only point where derivative of the dissipation potential is zero, $\dcstresssymb = \Xi_{\gradsym}|_{\gradsym=0} = \tensorq{0}$.
The conjugate dissipation potential $\Xi^*(\dcstresssymb)$, which is formed by Legendre transformation, is also convex and zero at the origin. Moreover, it holds that
\begin{equation}
\langle \gradsym, \dcstresssymb\rangle = \underbrace{\Xi(\gradsym)}_{\geq 0} + \underbrace{\Xi^*(\dcstresssymb)}_{\geq 0}\geq 0,
\end{equation}
see e.g. p. 267 in \cite{Roubicek}. Entropy production is the product of thermodynamic forces and fluxes as usually in  \cite{dGM} or in the GENERIC evolution equation for momentum \eqref{eq.evo.uu.GENERIC}, and this last inequality represents the positive sign of entropy production. Convexity of the dissipation potential thus implies compatibility with the second law of thermodynamics.

\textbf{On the second law and non-convexity}: Convexity of the dissipation potential is a sufficient condition for compatibility with the second law. However, it is not a necessary condition. The second law is fulfilled if and only if the dissipation potential generates non-negative entropy production,
\begin{equation}\label{eq.2nd}
  \langle \gradsym, \Xi_\gradsym\rangle \geq 0 \qquad\mbox{or}\qquad \langle \dcstresssymb, \Xi^*_\dcstresssymb\rangle \geq 0.
\end{equation}
Practically always in non-equilibrium thermodynamics this condition is required to hold locally (not only after integration over the whole system). Moreover, let us consider only the case where the dissipation potential depends only on a norm of the thermodynamic force or flux (as always in this paper). The local version of condition \eqref{eq.2nd} is then
\begin{equation}
  0 \leq \dcstresssymb : \Xi^*_\dcstresssymb = \absnorm{\dcstresssymb} \frac{\dcstresssymb}{\absnorm{\dcstresssymb}} : \Xi^*_{\absnorm{\dcstresssymb}} \pd{\absnorm{\dcstresssymb}}{\dcstresssymb} =
  \absnorm{\dcstresssymb} \frac{\dcstresssymb}{\absnorm{\dcstresssymb}} : \Xi^*_{\absnorm{\dcstresssymb}} \frac{\dcstresssymb}{\absnorm{\dcstresssymb}}
  = \absnorm{\dcstresssymb} \Xi^*_{\absnorm{\dcstresssymb}},
\end{equation}
which is satisfied for instance if the function $\Xi^*(|\dcstresssymb|)$ is monotonous, see Fig.~\ref{fig:dissipation-potential}. In this sense, a non-convex dissipation potential which is monotonous as a function of $\absnorm{\dcstresssymb}$ fulfills the second law of thermodynamics.

A generalization to dissipation potentials that are not necessarily simple functions of a norm of the forces or fluxes was found for instance in \cite{Mira-Milan} (Eqs. (28)-(30) therein).

It should be also borne in mind that condition \eqref{eq.2nd} expresses non-negativity of entropy production only if the thermodynamic potential whose derivatives generate the thermodynamic forces is convex. Therefore, we always assume that entropy is a concave function of the state variables and that the thermodynamic potential is convex.

We believe that the experimental observations can be understood also without non-convex dissipation potentials, but on more detailed levels of description. The non-convex dissipation potential should be a result of a reduction from a more detailed level of description equipped with a convex dissipation potential. It is not clear, however, what is the more detailed level of description. Perhaps a further extension in the sense of internal variables (that can be often seen as a CR-extension) could generate the same results using a convex dissipation potential. But unless such extension is known, let us proceed with a non-convex dissipation potential.

\subsection{Approach to equilibrium and existence}
Convexity of the dissipation potential near equilibrium, which can be interpreted as approximately quadratic behavior, is moreover in tight relation with stability of the equilibrium state. Apart from the convexity of the dissipation potential, convexity of the energy plays a crucial role in stability analysis. When these two requirements are satisfied, the equilibrium state is stable with respect to small fluctuations, as can be seen for example from analysis of the spectrum of the linearized GENERIC near equilibrium, see~\cite{PRE2014}, or as discussed in the paper by Matolcsi~\cite{Matolcsi-Lyapunov}. Thermodynamic equilibrium is stable if the dissipation potential is convex and the entropy concave near equilibrium.

On the other hand, although convexity of the dissipation potential (also monotonicity of the derivative) implies increase of the entropy by condition \eqref{eq.2nd}, it does not make up for an exact mathematical proof of the existence of solutions and of the approach to equilibrium from arbitrarily far-from-equilibrium states. It is, however, an important and essential ingredient from both the physical and the mathematical points of view. This can be well illustrated on the case of the Boltzmann equation where the exact mathematics has been done.

Boltzmann equation is an illustration of GENERIC with the non-quadratic but convex dissipation potential, see~\cite{GO}. The proof of the global existence of solutions has been done by Di Perna and Lions~\cite{Lions}, for which the latter author received the Fields Medal. The Boltzmann $H$-theorem is not sufficient but plays an essential role in the proof. The approach to the equilibrium state was also shown by Desvilletes and Villani \cite{Desvilettes-Villani}, where again the Boltzmann $H$-theorem does not suffice but plays an essential role.

Summing up, the increase of the entropy as implied by monotonicity of the dissipation potential does not, indeed, imply immediately all the mathematical rigor in the proofs of the existence of solutions and the approach to the equilibrium from any far-from-equilibrium state, but it is an essential ingredient in such proofs (that have been done only for a few mesoscopic time evolution equations so far).

\subsection{Shear and vorticity banding}

Let us now return to Fig. \ref{fig.exp}. It is possible to interpret the $\ensuremath{\tensorq{D}}({\traceless{\cstress}})$ dependence in Fig.~\ref{fig.exp} (open triangles) in terms of gradient dynamics. For each ${\traceless{\cstress}}$ one has a unique value of $\ensuremath{\tensorq{D}}$, which is the slope of the function $\Xi^*$ at the particular value of ${\traceless{\cstress}}$. In the next section, we try to answer how to interpret the dependence ${\traceless{\cstress}}(\ensuremath{\tensorq{D}})$.

  Before going on with the details on the conjugate dissipation potential, let us briefly comment on the relation of this manuscript to the gradient and vorticity banding phenomena, excellently reviewed in~\cite{Fielding,Olmsted}. Shear banding is usually experimentally observed in the cylindrical Couette flow (as in this paper), where there might develop instabilities in form of (visible) bands. If the bands develop in the radial direction, the phenomenon is referred to as gradient banding while if they develop in the direction of the rotational axis, the phenomenon is referred to as vorticity banding.

  Gradient banding usually leads to stress--strain curves with multiple shear rates corresponding to a single shear stress (''horizontal'' S-curve) while vorticity banding leads to curves where multiple shear stresses correspond to a single shear rate (''vertical'' S-curve). Although such distinction between the two phenomena is typical, it is not always true \cite{Olmsted}. When crossing the region of unstable shear rates in the case of gradient banding, the region of higher shear rates is usually characterized by lower effective viscosity, i.e. the fluid is shear thinning. On the other hand, vorticity banding often corresponds to increased viscosity---shear thickening fluids. However, there have been reported experiments where such distinction was not valid \cite{Olmsted}.

  Since it is difficult to characterize gradient/vorticity banding in the terms of horizontal/vertical S-curves or shear thinning/thickening fluids, we shall discuss the S-curves themselves. There are several models leading to horizontal S-curves for the stress--strain dependency (multiple shear rates for a single shear stress). For instance, the Johnson--Segalman model \cite{Johnson-Segalman} or rather its weakly non-local variant \cite{Olmsted-Radulescu}, thermodynamics with internal variables \cite{Verhas,Asszonyi2015}, viscoelastic models undergoing scission of polymeric chains \cite{Miroslav-micelles} or viscoelastic models with relaxation of internal structure \cite{Miroslav-suspensions,Germann-Beris} indicate the horizontal S-shape of the stress--strain relation. The literature is rather vast, and we do not attempt to make the list of references complete. Interested reader could consult, for example, recent special issue of the Journal of Rheology \cite{Fielding-special}. Most of the models usually exhibit shear thinning, and it is believed that shear thickening should be captured by models with non-linear dependence of the constitutive relation on the shear stress \cite{Radulescu}, which is a motivation for the constitutive relation~\eqref{eq:le-roux-rajagopal}. This paper thus aims at modeling of the vertical S-shaped stress--strain curve with shear thickening, typically observed in vorticity banding.

\section{Conjugate representation}
\label{sec.conjugate}

The purpose of this section is to present details on the Legendre duality between $\Xi$ and $\Xi^*$. A consequence of the non-convexity of $\Xi^*$ is that $\Xi$ is multivalued, i.e. there are multiple shear stresses corresponding to certain shear rates. The shear rates that correspond to multiple stresses then constitute the CR-metastable and CR-unstable regions of the dynamics, which are discussed in Sec.~\ref{sec.stability}. In fact, Sec.~\ref{sec.stability} contains the essence of this paper, and reading it before the following more technical parts could be advisable. Note also that Legendre transformation can be regarded as constrained extremization of the generating potential, which often turns out to be a practical method, see e.g. \cite{Junker2014}.

The balance of linear momentum reads
\begin{equation}
  \label{eq.u.evo}
  \pd{\mathbf{u}}{t} = \underbrace{-\nabla p - \divergence\left(\frac{\mathbf{u}\otimes\mathbf{u}}{\rho}\right)}_{=\mbox{reversible terms}} + \divergence({\traceless{\cstress}}),
\end{equation}
where $\mathbf{u}$ is the density of momentum (per volume). The irreversible part of this evolution equation is in the form of divergence of the shear stress ${\traceless{\cstress}}$ (thermodynamic flux). Shear rate $\ensuremath{\tensorq{D}}$ is thus interpreted as a thermodynamic force. The reversible terms, see Appendix \ref{sec.rev} for the definition of reversibility and irreversibility, are the gradient of the pressure and the divergence of convective terms, and their Hamiltonian nature is discussed in App.~\ref{sec.hydro}.

Only the field of momentum density plays the role of state variable, since we assume that the fluid is incompressible ($\divergence \ensuremath{\vec{v}} = 0$) and isothermal. When carrying out the CR-extension below (Sec. \ref{sec.promotion}), also the field of conjugate stress tensor becomes a state variable.

\begin{subequations}\label{eq.LT.Xic.Xi}
Equation \eqref{eq.gradient} can be then rewritten as
\begin{equation}
  \label{eq.Xic.JJ}
  \left. \Xi^*_\dcstresssymb \right|_{\tilde{\cstress}_\delta(\gradsym)} = \gradsym,
\end{equation}
and can be regarded as an equation for $\dcstresssymb$. Denoting the solutions to this equation by $\tilde{\cstress}_\delta (\gradsym)$, the equation can be further interpreted as
\begin{equation}
  \label{eq.Xic.Xi}
  \pd{}{\dcstresssymb} \left( -\Xi^*(\dcstresssymb) + \langle \dcstresssymb, \gradsym \rangle \right) \Big|_{\tilde{\cstress}_\delta (\gradsym)} = \tensorq{0},
\end{equation}
where $\langle\bullet,\bullet\rangle$ is the $L^2$ scalar product---spatial integral of product of the two arguments. The last equation leads to the formulation of the dissipation potential $\Xi$ conjugate to $\Xi^*$
\begin{equation}
  \label{eq.Xi}
  \Xi(\gradsym) =_\mathrm{def} \left( -\Xi^*(\dcstresssymb) + \langle \dcstresssymb,\gradsym \rangle \right) \Big|_{\tilde{\cstress}_\delta (\gradsym)}.
\end{equation}
Note that we do not require the solution to be unique, which means that we do not interpret $\tilde{\cstress}_\delta (\gradsym)$ as a function, but rather as a graph in the $(\gradsym,\dcstresssymb)$ plane or as a collection of functions. Therefore, dissipation potential $\Xi$ is multivalued as can be seen from Fig.~\ref{fig:conjugate-dissipation-potential}. Although the dissipation potential is multivalued, it is uniquely determined and the generalized (multivalued) Legendre transformation is invertible\footnote{Invertible in the sense that the backward Legendre transformation restores the original function (or rather graph).} \cite{dorst.l.boomgaard.r:analytical,dorst.l.boomgaard.r:morphological}.
\end{subequations}

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.6\textwidth]{conjugate-dissipation-potential.pdf}
  \caption{Conjugate dissipation potential $\Xi \left( \ensuremath{\tensorq{D}} \right)$ corresponding to the constitutive relation~\eqref{eq:le-roux-rajagopal} is obtained from the dissipation potential $\Xi^* \left( {\traceless{\cstress}} \right)$ (Fig.~\ref{fig:dissipation-potential}) using the generalized Legendre transformation~\eqref{eq.Xi}. In regions II and III, the potential is multivalued. The dashed line $C-D$ is the CR-unstable branch, lines $B-C$ and $D-E$ are CR-metastable branches, and the curve $A-B,E-F$ is globally CR-stable and corresponds to the convex hull of the dissipation potential $\Xi^*$. See Sec.~\ref{sec.stability} for detailed explanation.}
  \label{fig:conjugate-dissipation-potential}
\end{figure}

\begin{subequations}\label{eq.LT.Xi.Xic}
Using Eq.~\eqref{eq.Xic.JJ}, we obtain from Eq.~\eqref{eq.Xi} that
\begin{equation}
  \label{eq.Xi.XX}
  \Xi_{\gradsym} = \tilde{\cstress}_\delta (\gradsym),
\end{equation}
solution of which are denoted by $\tilde{\gradsym}(\dcstresssymb)$. These solutions can be regarded as the same graph in the $(\gradsym,\dcstresssymb)$ plane as solutions $\tilde{\cstress}_\delta (\gradsym)$. An example of such a graph, which expresses simply the constitutive relation is plotted in Fig.~\ref{fig:flow-curve}.
Equation~\eqref{eq.Xi.XX} can be rewritten as
\begin{equation}
  \pd{}{\gradsym} \left(-\Xi(\gradsym) + \langle\gradsym,\dcstresssymb\rangle\right)\Big|_{\tilde{\gradsym}(\dcstresssymb)} = \tensorq{0},
\end{equation}
 which finally leads back to the original dissipation potential,
\begin{equation}
  \Xi^*(\dcstresssymb) = -\Xi \left(\tilde{\gradsym}(\dcstresssymb)\right) + \left\langle\tilde{\gradsym}(\dcstresssymb), \dcstresssymb\right\rangle,
\end{equation}
derivative of which gives Eq.~\eqref{eq.Xic.JJ}. We have thus performed Legendre transformation from $\Xi^*$ to $\Xi$ and back. Interested reader can consult App.~\ref{sec.Legendres} for other forms of Legendre-like transformations.
\end{subequations}

The graph $\tilde{\cstress}_\delta (\gradsym)$ or $\tilde{\gradsym}(\dcstresssymb)$ is plotted in Fig.~\ref{fig:flow-curve}, where also the lever rule is indicated. Of course, on one hand, we have
\begin{equation}
  \int_{|{\traceless{\cstress}}|(B)}^{|{\traceless{\cstress}}|(E)} |\ensuremath{\tensorq{D}}(T)| \,\mathrm{d} T = \int_{|{\traceless{\cstress}}|(B)}^{|{\traceless{\cstress}}|(E)} \Xi^*_{|{\traceless{\cstress}}|(T)} \,\mathrm{d} T = \Xi^*(|{\traceless{\cstress}}(E)|)-\Xi^*(|{\traceless{\cstress}}(B)|),
\end{equation}
while on the other hand from Fig.~\ref{fig:dissipation-potential} it is clear that
\begin{equation}
  \frac{\Xi^*({\traceless{\cstress}}(E))-\Xi^*({\traceless{\cstress}}(B))}{|{\traceless{\cstress}}(E)|-|{\traceless{\cstress}}(B)|} = \Xi^*_{|{\traceless{\cstress}}(E)|}= \Xi^*_{|{\traceless{\cstress}}(B)|}.
\end{equation}
Combining these two results leads to
\begin{equation}
  \int_{|{\traceless{\cstress}}|(B)}^{|{\traceless{\cstress}}|(E)} |\ensuremath{\tensorq{D}}(T)| \,\mathrm{d} T =(|{\traceless{\cstress}}|(E) - |{\traceless{\cstress}}|(B))\:\underbrace{|\ensuremath{\tensorq{D}}(|{\traceless{\cstress}}(E)|)|}_{=|\ensuremath{\tensorq{D}}(|{\traceless{\cstress}}|(B)|)|},
\end{equation}
which means that the shaded areas in Fig.~\ref{fig:flow-curve} are equal---the Maxwell lever rule  known from equilibrium thermodynamics \cite{Callen}.

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.6\textwidth]{flow-curve.pdf}
  \caption{Flow curve (the constitutive relation) given by Eq.~\eqref{eq:le-roux-rajagopal} with parameter values $\alpha = 1$, $\beta = 1$, $\gamma = 0.02$ and $n = -2$. The line $B-E$ is constructed using the Maxwell lever rule, so that the shaded areas are equal. Since the segment $C-D$ corresponds to the region where the dissipation potential $\Xi^*$ loses convexity, it is CR-unstable, as commented in Sec.~\ref{sec.stability}. The lines $B-C$ and $D-E$ are CR-metastable---they coincide with the lower CR-stable branches of the dissipation potential $\Xi$, see Fig~\ref{fig:conjugate-dissipation-potential}. Note that the vertical axis expresses the strain rate, in contrast to the usual habit, and the S-curve in this figure is thus referred to as a vertical. }
  \label{fig:flow-curve}
\end{figure}

The Maxwell lever rule plays an analogical role in equilibrium phase transitions. The hysteresis is observed if fluctuations are sufficiently small. If they are large enough, the system evolves along the straight Maxwell line. The same is true in the case of the non-convex dissipation potential.

The Legendre duality between $\Xi$ as a function of forces and $\Xi^*$ as a function of fluxes is often advantageous because one of the forms can be more accessible. The former choice is used for example in chemical kinetics \cite{Grchemkin} while the latter form for example used in plasticity modeling \cite{Mielke2003}.

\textbf{Stability of constitutive relation}: At many places within this paper, we discuss CR-stability of constitutive relations generated by non-convex dissipation potentials. By CR-stability of constitutive relations we mean whether the constitutive relation is CR-stable with respect to perturbations. If the dissipation potential is convex everywhere, the dual (Legendre transformed) dissipation potential is convex as well, and for each value of the thermodynamic force there is only one value of the flux. On the other hand, if the potential is non-convex, the dual dissipation potential is multi-valued and for some values of thermodynamic forces (or fluxes) there are several possible values of fluxes (or forces). The graph of the constitutive relation then has several branches as in Fig.~\ref{fig:conjugate-dissipation-potential}. In such regions where multiple admissible fluxes (or forces) are present, perturbation of the flux can cause evolution of the system towards another branch of the graph of the constitutive relation. The constitutive relation is said to be CR-unstable if any perturbation leads to selection of a different branch, CR-metastable if only large enough perturbation switches the system to a different branch, and CR-stable if the current branch of the constitutive relation is not changed by any perturbation of the flux.

\section{Extended hydrodynamics (CR-extension)}
\label{sec.CR-thermodynamics}

  In this section the set of hydrodynamic state variables\footnote{The state variables are $(\rho,\mathbf{u},s)$ in the energetic representation or $(\rho, \mathbf{u}, e)$ in the entropic representation.} is expanded by promoting shear stress to an independent state variable, which is the reason for the name extended hydrodynamics. The idea of extending the set of state variables is ubiquitous in non-equilibrium thermodynamics, for example Extended Irreversible Thermodynamics (EIT) \cite{Jou-EIT} or thermodynamics with internal variables \cite{Van-Berezovski}. For simplicity, only the isothermal setting is considered in this paper. The isothermal hydrodynamic evolution is thus enriched to a form of extended hydrodynamics. The extra state variable has its own evolution equation, which will be shown to be the derivative of a \textbf{multiscale thermodynamic Lyapunov} (MTL) function with respect to the conjugate of the extra state variable (Sec.~\ref{sec.fast}). The extra state variable then evolves towards the maxima of the Lyapunov function. The figures in Sec.~\ref{sec.stability} show how the Lyapunov function depends on the shear rate and towards what values the shear stress evolves. This will give a graphical representation of the CR-stability analysis of the system. Motivation for reading this rather technical section can thus be found in Sec.~\ref{sec.stability}.

Dissipation potential $\Xi^*$ is convex near equilibrium and for high shear stresses, but it loses convexity in the intermediate range of the shear stress. Convexity leads to existence of a Lyapunov functional, see e.g. \cite{GO}, and Braun--Le Chatelier principle \cite{Braun-LeChat}, that guarantees stability of the evolution. The former paper is crucial for understanding everything that follows. Instability can be thus anticipated in the non-convex regions. Let us now have a look at a non-equilibrium thermodynamic description of stability, metastability and instability. The analysis is formulated within a recent framework developed by Grmela~\cite{MG-CR}, the CR-thermodynamics (CR stands for constitutive relation), which is a form of mesoscopic multiscale thermodynamics.\footnote{We use a slightly different notation compared with the original paper \cite{MG-CR}. Instead of thermodynamic potentials, we use the negatives of the potentials, and we call them MTL functions. The reason is that we can then always use the same form of Legendre transformation. Another difference is that the extra variable is denoted by $\cstress_\delta^*$ instead of $\dcstresssymb$ so that the usual meaning of flux is kept.}

\subsection{GENERIC structure of hydrodynamics}

Evolution of momentum density is governed by Eq.~\eqref{eq.u.evo}, which can be rewritten as (within the GENERIC framework)
\begin{equation}
  \label{eq.evo.uu.GENERIC}
  \pd{\mathbf{u}}{t} = \left( \pd{\mathbf{u}}{t} \right)_{rev} + \left. \Xi_{\mathbf{u}^*} \right|_{\mathbf{u}^* = \Psi^{(CH\rightarrow E)}_\mathbf{u}}.
\end{equation}
The reversible part is Hamiltonian, see in Appendix \ref{sec.hydro}, and it expresses the Euler equations of compressible hydrodynamics. The MTL function (the negative of the thermodynamic potential from \cite{MG-CR}) driving classical hydrodynamics (CH) to equilibrium (E) is
\begin{equation}
  \Psi^{(CH\rightarrow E)} = \int \mathrm{d}\mathbf{r}\, S(\rho,\mathbf{u},e) - \frac{1}{T_0} E(\rho,\mathbf{u},e) +\frac{\mu_0}{T_0} M(\rho,\mathbf{u},e),
\end{equation}
where $e$ is total energy density (per volume).
Total entropy is
\begin{equation}\label{eq.S.hydro}
  S = \int \mathrm{d}\mathbf{r}\, s\left(\rho, e - \frac{\mathbf{u}^2}{2\rho} \right),
\end{equation}
where $s(\rho,\varepsilon)$ denotes local equilibrium entropy density (which depends on density and internal energy density $\varepsilon=e-\frac{\mathbf{u}^2}{2\rho}$). Total energy is $E=\int\mathrm{d}\mathbf{r}\,e(\mathbf{r})$, and $M=\int\mathrm{d}\mathbf{r}\, \rho$ is the total mass. Constant $T_0$ is the temperature an isolated system would have after relaxation to equilibrium and constant $\mu_0$ the equilibrium chemical potential. Local temperature is defined as the derivative of entropy with respect to the total energy density, i.e. $T^{-1}(\mathbf{r}) =_\mathrm{def} S_{e(\mathbf{r})}$. Derivative of the MTL function with respect to momentum is then
\begin{equation}
  \mathbf{u}^* = \left( \pd{\Psi^{(CH\rightarrow E)}}{\mathbf{u}} \right)_{\rho,s} = S_\mathbf{u} = -\frac{1}{T}  \frac{\mathbf{u}}{\rho}.
\end{equation}
Assuming isothermal conditions, the conjugate momentum, $\mathbf{u}^*$, is proportional to velocity $\vecv = \mathbf{u}/\rho$.

The dissipation potential is often in the form
\begin{equation}
  \Xi = \Xi \left( \Gamma(\mathbf{u}^*) \right),
\end{equation}
$\Gamma$ being a linear operator. In our case
\begin{equation}
  \label{eq:operator-gamma}
  \Gamma(\bullet) = -\frac{1}{2} \left(\nabla \bullet +(\nabla \bullet)^\top \right),
\end{equation}
which leads to shear rate $\Gamma(\mathbf{u}^*)_{\mathbf{u}^*=\Psi_\mathbf{u}^{(CH\rightarrow E)}} = \frac{1}{T_0} \ensuremath{\tensorq{D}}$. The constant prefactor $T_0$ is going to be ignored for the sake of simplicity. In a slightly overloaded notation, we can then write
\begin{equation}
  \Xi_{\mathbf{u}^*} = \left\langle \Xi_{\gradsym},\pd{\gradsym}{\mathbf{u}^*} \right\rangle
  =\left\langle \Xi_{\gradsym}, \pd{\Gamma(\mathbf{u}^*)}{\mathbf{u}^*}\right\rangle
  =\transpose{\Gamma} \left(\Xi_{\gradsym}\right),
\end{equation}
where the adjoint operator $\transpose{\Gamma}$ is defined by
\begin{equation}
  \int \mathrm{d} \mathbf{r}\, \Gamma(\mathbf{a}) \cdot \mathbf{b} = \int \mathrm{d}\mathbf{r}\, \mathbf{a} \cdot \transpose{\Gamma} (\mathbf{b}), \qquad \forall \mbox{ symmetric tensor fields } \mathbf{a}, \mathbf{b}.
\end{equation}
This also specifies the scalar product $\langle\bullet,\bullet\rangle$ of two functions as the spatial integral of product of the functions ($L^2$ scalar product). In our case, we have $\transpose{\Gamma} (\bullet) = \divergence (\bullet)$. Eventually, the evolution equation for momentum density becomes
\begin{equation}\label{eq.evo.uu.Gamma}
  \pd{\mathbf{u}}{t} = \left( \pd{\mathbf{u}}{t} \right)_{rev} + \transpose{\Gamma} (\dcstresssymb),
\end{equation}
where $\dcstresssymb = \Xi_\gradsym$. Note that the shear stress ${\traceless{\cstress}}$ is a thermodynamic flux because it is conjugate to the force ($\gradsym$), which comes from differentiation of the MTL function (or entropy). Equation~\eqref{eq.evo.uu.Gamma}, which is the standard evolution equation for momentum density in fluid dynamics, is thus fully compatible with the GENERIC expression \eqref{eq.evo.uu.GENERIC}.

\subsection{Promoting conjugate stress to a new state variable}
\label{sec.promotion}

The goal is now to promote the flux to an independent state variable. But $\dcstresssymb$ can not be promoted itself because derivative of a functional with respect to $\dcstresssymb$ would then have to be on right hand side of evolution equation for $\mathbf{u}$ instead of $\dcstresssymb$ itself. Rather, a conjugate flux $\cstress_\delta^*$ is the new variable, that is coupled with momentum density. Since $\dcstresssymb$ is interpreted the shear stress, $\cstress_\delta^*$ is interpreted as conjugate shear stress. As usual in non-equilibrium thermodynamics \cite{RedExt}, the dissipation can be moved onto a higher (more detailed) level of description, i.e. into the evolution equation for the extra state variable. Dissipative evolution of momentum density can be then recovered by letting the extra variable relax. The reversible coupling between the original state variables and the extra state variable is required to be antisymmetric in order to fulfill the Onsager--Casimir reciprocal relations \cite{Casimir1945,PRE2014,HCO} and thus we write (still following \cite{MG-CR}) evolution equations for extended hydrodynamics (EH) as follows
\begin{subequations}
  \label{eq.EH.evo}
  \begin{eqnarray}
    \label{eq.u.evo.EH}
    \pd{\mathbf{u}}{t} &=& \left( \pd{\mathbf{u}}{t} \right)_{rev} + \transpose{\Gamma} (\cstress_\delta^{**}),\\
    \label{eq.JJc.evo}
    \pd{\cstress_\delta^*}{t} &=& - \Gamma(\mathbf{u}^*) + \Xi^*_{\cstress_\delta^{**}}.
  \end{eqnarray}
\end{subequations}
To close the evolution equations, the conjugates are evaluated as the corresponding derivatives of a new extended MTL function $\Psi^{(EH \to E)}$, $\mathbf{u}^* = \Psi^{(EH \rightarrow E)}_{\mathbf{u}}$ and $\cstress_\delta^{**} = \dcstresssymb = \Psi^{(EH \rightarrow E)}_{\cstress_\delta^*}$. Note that $\cstress_\delta^{**}$ was identified with $\dcstresssymb$ because the generalized Legendre transformation is idempotent (second application of the transformation restores the original functional), i.e. $\cstress_\delta^{**} = \dcstresssymb$. The MTL function $\Psi^{(EH \rightarrow E)}$, which drives evolution of the extended hydrodynamics towards equilibrium, has yet to be specified in order to define the conjugate stress tensor completely. Evolution equations \eqref{eq.EH.evo} represent evolution on the level of extended hydrodynamics.

Parity with respect to the time-reversal \cite{PRE2014} of $\dcstresssymb$ is now even, which is in contrast to the odd parity of $\dcstresssymb$ before promoting $\cstress_\delta^*$ to a state variable. The physical reason is that we have now a more detailed description of the shear stress, which can be regarded for example as the stress formulated by kinetic theory (App.~\ref{sec.KT}). Similarly, $\cstress_\delta^*$ and $\cstress_\delta^{**}$ are even with respect to time reversal, since MTL functions and dissipation potentials, with respect to which the fluxes are conjugate, are always even. Therefore, evolution equation \eqref{eq.u.evo.EH} is now reversible (in contrast to Eq.~\eqref{eq.u.evo}) while equation \eqref{eq.JJc.evo} contains a reversible term (first) and an irreversible term (second). By moving to a more detailed level of description, parity of the quantity promoted to an extra state variable may change because the quantity becomes independent and is no longer enslaved by the constitutive relation valid before the promotion. As a consequence, irreversible terms may become reversible.

\subsection{Specification of the MTL function}
The conjugate stress will be now specified more precisely. The relation between $\dcstresssymb$ and $\cstress_\delta^*$ is a Legendre transformation,
\begin{equation}
  \label{eq.LT.JJ.JJc}
  \dcstresssymb = \Psi^{(EH\rightarrow E)}_{\cstress_\delta^*}, \quad \text{ or } \quad
  \pd{}{\cstress_\delta^*} \left( - \Psi^{(EH\rightarrow E)} + \int \mathrm{d} \mathbf{r}\, \dcstresssymb : \cstress_\delta^* \right) = \tensorq{0},
\end{equation}
which is carried out with respect to the MTL function
\begin{equation}\label{eq.Psi.EH.E}
  \Psi^{(EH \rightarrow E)} =_\mathrm{def} S^{(EH\rightarrow E)} - \frac{1}{T_0}E^{(EH\rightarrow E)} + \frac{\mu_0}{T_0}M^{(EH\rightarrow E)}.
\end{equation}
This MTL function drives the extended hydrodynamics towards equilibrium.  Note that all quantities with $(EH\rightarrow E)$ superscript depend on the hydrodynamic fields $\mathbf{x}^{(CH)} =_\mathrm{def} (\rho,\mathbf{u},e)$ and the conjugate stress tensor $\cstress_\delta^*$. Quantities $E^{(EH\rightarrow E)}$, $S^{(EH\rightarrow E)}$ and $M^{(EH\rightarrow E)}$, which have so far been left undetermined, are the building blocks of the corresponding MTL function, and they can reflect the particular nature of the physical system. Alternatively, the extra evolution equation can be interpreted as a numerical scheme leading to constitutive relation \eqref{eq:le-roux-rajagopal}.

Let us consider at least one possible choice of $\Psi^{(EH\rightarrow E)}$. Within the kinetic theory (App.~\ref{sec.KT}), kinetic energy of particles is equal to the sum of kinetic energy of the overall motion $\mathbf{u}^2/(2\rho)$ and pressure, and the deviatoric part of stress tensor thus does not contribute to energy (when neglecting long-range pair interactions). On the other hand, additional knowledge expressed as the extra state variable reduces entropy, and the entropy on the $(EH)$ level is thus lower than the entropy on the level of classical hydrodynamics $(CH)$. Similarly as in Extended Irreversible Thermodynamics \cite{Jou-EIT}, the first approximation of the $(EH)$ entropy is
\begin{equation}
  \label{eq.EH.S}
  S^{(EH)} = S^{(CH)} -\int \mathrm{d}\mathbf{r}\, \frac{1}{2}\beta \absnorm{\cstress_\delta^*}^2 = \int\mathrm{d}\mathbf{r}\, s\left(\rho,e-\frac{\mathbf{u}^2}{2\rho}\right) - \frac{1}{2}\beta \absnorm{\cstress_\delta^*}^2.
\end{equation}
Coefficient $\beta$ is positive so that the entropy is concave. Energy and mass on the $(EH)$ level are the same as on the $(CH)$ level, thus the MTL function driving evolution from $(EH)$ to equilibrium is
\begin{equation}
  \Psi^{(EH\rightarrow E)} = \Psi^{(CH\rightarrow E)} -\int\mathrm{d}\mathbf{r}\, \frac{1}{2} \beta\absnorm{\cstress_\delta^*}^2,
\end{equation}
from which it follows that
\begin{equation}\label{eq.T.Tc}
  \dcstresssymb = \Psi^{(EH\rightarrow E)}_{\cstress_\delta^*} = - \beta \cstress_\delta^*.
\end{equation}
In summary, for the particular choice of the entropy \eqref{eq.EH.S}, the conjugate stress tensor $\cstress_\delta^*$ is proportional to the standard shear stress tensor ${\traceless{\cstress}}$   and the evolution equations~\eqref{eq.EH.evo} for the constitutive relation~\eqref{eq:le-roux-rajagopal} and the operator $\Gamma$ defined by~\eqref{eq:operator-gamma} read
\begin{subequations}

  \begin{align}
    \pd{\ensuremath{\vec{u}}}{t} &= \left( \pd{\ensuremath{\vec{u}}}{t} \right)_{rev} + \divergence {\traceless{\cstress}}, \\
    - \frac{1}{\beta} \pd{{\traceless{\cstress}}}{t} &= - \ensuremath{\tensorq{D}} + \left[ \alpha (1 + \beta \absnorm{{\traceless{\cstress}}}^2)^n + \gamma \right] {\traceless{\cstress}}. \label{eq:26b}
  \end{align}
\end{subequations}
The second equation can be apparently interpreted as the time derivative regularization of the original constitutive relation.

\subsection{Duality between $\dcstresssymb$ and $\cstress_\delta^*$}
\label{sec.EH.duality}

The inverse Legendre transformation from $\cstress_\delta^*$ to $\dcstresssymb$ follows. The conjugate MTL function $\Psi^{*,(EH\rightarrow E)}$ can be constructed as
\begin{equation}
  \Psi^{*,(EH\rightarrow E)} (\mathbf{u}^*,\dcstresssymb) = - \Psi^{(EH\rightarrow E)} \left( \mathbf{u}, \cstress_\delta^*(\dcstresssymb) \right) + \int \mathrm{d}\mathbf{r}\, \cstress_\delta^*(\dcstresssymb) : \dcstresssymb
\end{equation}
where $\cstress_\delta^*(\dcstresssymb)$ is the solution to \eqref{eq.LT.JJ.JJc}, i.e. Eq. \eqref{eq.T.Tc} for the particular choices of entropy and energy. Conversely, the Legendre transformation
\begin{equation}
  \label{eq.LT.JJc.JJ}
  \cstress_\delta^* = \Psi^{*,(EH\rightarrow E)}_{\dcstresssymb}, \quad \text{ or } \quad
  \frac{\partial}{\partial \dcstresssymb}\left(-\Psi^{*,(EH\rightarrow E)} + \int \mathrm{d} \mathbf{r}\, \dcstresssymb : \cstress_\delta^* \right) = \tensorq{0}
\end{equation}
leads back to the original MTL function,
\begin{equation}
  \Psi^{(EH\rightarrow E)} (\mathbf{u}, \cstress_\delta^*) = -\Psi^{*,(EH\rightarrow E)} \left( \mathbf{u}^*, \dcstresssymb(\cstress_\delta^*) \right) + \int \mathrm{d} \mathbf{r}\, \dcstresssymb(\cstress_\delta^*):\cstress_\delta^*,
\end{equation}
where $\dcstresssymb(\cstress_\delta^*)$ are the solutions to Eq.~\eqref{eq.LT.JJc.JJ}. The Legendre duality between $\dcstresssymb$ and $\cstress_\delta^*$ has thus been established.

\subsection{Fast reducing evolution}
\label{sec.fast}

Now we are coming to the highlight of the thermodynamic extension. Let us now identify the evolution of $\cstress_\delta^*$ that leads back to classical hydrodynamics, i.e. the reducing (fast) evolution leading to the reduced (hydrodynamic) evolution of variables $\mathbf{x}$. Equation \eqref{eq.JJc.evo} can be rewritten as
\begin{equation}
  \label{eq.reducing}
  \pd{\cstress_\delta^*}{t} = - \pd{}{\dcstresssymb} \left( \left\langle \Gamma(\mathbf{u}^*), \dcstresssymb \right\rangle - \Xi^* \right) \Big|_{\left(\mathbf{u}^* = \Psi^{(EH \rightarrow E)}_{\mathbf{u}}, \dcstresssymb= \Psi^{(EH \rightarrow E)}_{\cstress_\delta^*}\right)} = - \Psi^{(EH\rightarrow CH)}_{\dcstresssymb},
\end{equation}
where
\begin{equation}
  \label{eq.Psi.EH.CH}
  \Psi^{(EH\rightarrow CH)}(\gradsym,\dcstresssymb) =_\mathrm{def} -\Xi^*(\dcstresssymb) + \int \mathrm{d} \mathbf{r}\, \gradsym : \dcstresssymb,
\end{equation}
with $\gradsym = \Gamma(\mathbf{u}^*)$, is the MTL function driving evolution from extended hydrodynamics to classical hydrodynamics. This MTL function will be referred to as the \textbf{reducing MTL function}. Equation \eqref{eq.reducing} then means that conjugate shear stress $\cstress_\delta^*$ evolves in the direction of steepest ascent of $\Psi^{(EH\rightarrow CH)}(\gradsym,\dcstresssymb)$. MTL function \eqref{eq.Psi.EH.CH} generates Legendre transformation from $\Xi^*$ to $\Xi$ as
\begin{equation}
  \Psi^{(EH\rightarrow CH)}_\dcstresssymb = \tensorq{0},
\end{equation}
which is equivalent to Eq.~\eqref{eq.Xic.Xi}. The constitutive relation is thus recovered after the extra variable ($\cstress_\delta^*$) has relaxed to a state enslaved by the slow variable ($\mathbf{u}$).

The reducing evolution \eqref{eq.reducing} can be alternatively rewritten as
\begin{equation}\label{eq.reducing.J}
  \pd{\dcstresssymb}{t} = - \left\langle \frac{\partial^2 \Psi^{(EH\rightarrow E)}}{\partial \cstress_\delta^* \partial \cstress_\delta^*}, \Psi^{(EH\rightarrow CH)}_{\dcstresssymb} \right\rangle,
\end{equation}
where Eq.~\eqref{eq.LT.JJ.JJc} was used, and compatibility of the reducing evolution \eqref{eq.reducing} with the second law of thermodynamics is guaranteed by
\begin{equation}
  \dot{\Psi}^{(EH\rightarrow CH)} = -\int \mathrm{d} \mathbf{r}\, \Psi^{(EH\rightarrow CH)}_{\dcstresssymb} \Psi^{(EH\rightarrow E)}_{\cstress_\delta^* \cstress_\delta^*} \Psi^{(EH\rightarrow CH)}_{\dcstresssymb} \geq 0,
\end{equation}
for $\Psi^{(EH\rightarrow E)}$ concave. Note that only derivatives with respect to the fast variable are taken into account because only evolution of the fast variable is considered in the reducing evolution. The MTL function $\Psi^{(EH\rightarrow E)}$ has to be thus chosen concave, and the reducing evolution is compatible with the second law of thermodynamics.

On the other hand, dissipation potentials $\Xi$ and $\Xi^*$ were not convex and thus
\begin{equation}
  \dot{\Psi}^{(CH\rightarrow E)} = \int \mathrm{d} \mathbf{r}\, \Psi^{(CH\rightarrow E)}_{\mathbf{u}} \left. \Xi^*_{\mathbf{u}^*} \right|_{\mathbf{u}^* = \Psi^{(CH\rightarrow E)}}
\end{equation}
does not need to have definite sign and evolution on the hydrodynamic level can have apparent negative entropy production. Indeed, if shear rate and shear stress were spatially homogeneous with only for example $xy$ component nonzero, the slope of $\Xi^*$ could be negative for positive $\tensorc{D}_{xy}$. The way out of such a problem is the lift to extended hydrodynamics, where entropy production keeps its positivity. In other words, apparent negative entropy production would indicate evolution of instability (phase transition), and can be regarded as an invitation to a higher (more detailed) level of description.

After having demonstrated how the constitutive relation can be interpreted as a consequence of a more detailed evolution, study of behavior of the more detailed evolution gives information about CR-stability of the evolution with respect to the constitutive relation itself. The constitutive relation is stable if any perturbation decreases in time. It is metastable if only small perturbations decrease in time, and it is unstable if no perturbation decreases in time. The relation between CR-stability and the usual mathematical concept of stability is discussed in \ref{sec.stab-CRstab}.

\subsection{Alternative forms of the extension}

Evolution equations \eqref{eq.EH.evo} are relatively simple and useful for qualitative analysis of CR-stability carried out in Sec. \ref{sec.stability}, but they are not objective due to the presence of partial time derivative of the conjugate stress instead of an objective derivative. Moreover, Jacobi identity is not fulfilled for the system of equations (checked by program \cite{Kroeger2010}). Hence, the evolution equations have to be regarded rather as a toy model useful for explaining construction of the extension and the CR-stability implications.

The extension can be made more physically sound by choosing a better operator $\Gamma$. Such a choice, however, must be accompanied by an another choice of the quantity to be promoted to an extra state variable. For example, the choice
\begin{equation}
  \label{eq.Gamma.c}
  \divergence {\tensorq{T}}= \Gamma^T(\mathbf{c}^*) = T_0 c_{jk}\partial_i c^*_{jk} -T_0 \partial_k\left( c_{kj} (c^*_{ij}+c^*_{ji})\right)
\end{equation}
with $\mathbf{c}$ being the new variable (conformation tensor) leads to the adjoint operator
\begin{equation}
  \Gamma(\mathbf{u}^*) = -T_0\partial_k(c_{ij} u^*_k) + T_0 c_{kj}\partial_k u^*_i + c_{ki}\partial_k u^*_j.
\end{equation}
It should be stressed that this extension is considered in the energetic representation, where entropy density is among the state variables instead of the total energy density. Let us refer to this extension as to the conformation tensor (CT) level.

The evolution equation for the extra variable is then
\begin{equation}
  \label{eq.evo.c}
  \frac{\partial c_{ij}}{\partial t} = -\Gamma(\mathbf{u}^*) + \frac{\partial \Xi^*}{\partial c^*_{ij}}
  =\partial_k(c_{ij} T_0 u^*_k) - c_{kj}T_0 \partial_k u^*_i - c_{ki}T_0 \partial_k u^*_j + \frac{\partial \Xi^*}{\partial c^*_{ij}},
\end{equation}
which is analogical to Eq. \eqref{eq.EH.evo}. The reversible part is the upper-convected derivative of $\mathbf{c}$ (with an extra compressibility term, which is sometimes referred as to the Truesdell rate). Indeed, in the energetic representation total entropy depends only on the entropy density and derivative of MTL function \eqref{eq.Psi.EH.E} becomes
\begin{equation}
  \Psi^{(CT\rightarrow E)}_{\mathbf{u}} = -\frac{1}{T_0} E_\mathbf{u} = -\frac{1}{T_0}\frac{\mathbf{u}}{\rho},
\end{equation}
and Eq. \eqref{eq.evo.c} then becomes
\begin{equation}
  \frac{\partial c_{ij}}{\partial t}
  + \partial_k(c_{ij} v_k) - c_{kj} \partial_k v_i - c_{ki} \partial_k v_j =  \frac{\partial \Xi^*}{\partial c^*_{ij}},
\end{equation}
where the left hand side is the upper-convective derivative of $\mathbf{c}$. Note that velocity is given by $\mathbf{u}/\rho$. Such an extended evolution also fulfills the Jacobi identity, is thus Hamiltonian, see e.g. \cite{Gn1}, and is also objective. Objectivity seems to be in tight relation with validity of Jacobi identity. Evolution equation \eqref{eq.evo.c} is the analogy of Eq.~\eqref{eq.EH.evo} and this whole paragraph is the analogy of Sec.~\ref{sec.promotion}.

Conjugate dissipation potential $\Xi^*$ is defined through Legendre transformation
\begin{subequations}
  \begin{equation}
    \pd{}{\mathbf{u}^*} \left(-\Xi(\mathbf{u}^*) + \langle \mathbf{u}^*, \Gamma^T(\mathbf{c}^*)\rangle\right) = \tensorq{0} \quad\Rightarrow\quad \mathbf{u}^* = \mathbf{u}^*(\mathbf{c}^*),
  \end{equation}
  and
  \begin{equation}
    \Xi^*(\mathbf{c}^*) = - \Xi \left( \mathbf{u}^*(\mathbf{c}^*) \right) + \langle \Gamma(\mathbf{u}^*(\mathbf{c}^*)), \mathbf{c}^*\rangle,
  \end{equation}
\end{subequations}
which is analogical to the transformations carried out in Sec. \ref{sec.conjugate}.

The relation between $\mathbf{c}$ and $\mathbf{c}^*$ can be obtained as in Sec.~\ref{sec.EH.duality}, i.e. through Legendre transformation with respect to a MTL function $\Psi^{(CT\rightarrow E)}$. The MTL function is constructed from entropy, energy and total mass (as in Eq.~\eqref{eq.Psi.EH.E}), see e.g. \cite{Sarti-Marrucci,Gn1}. A particular choice of entropy could be
\begin{equation}
  \label{eq.CT.S}
  S (\rho, \mathbf{u}, e, \mathbf{c})^{(CT)} = \int \mathrm{d}\mathbf{r}\, s \left( \rho,e-\frac{\mathbf{u}^2}{2\rho} - H \Tr(\mathbf{c})\right) + \frac{1}{2}k_B\frac{\rho}{m} \ln \det \frac{\mathbf{c}}{Q} - \frac{1}{2}\beta \partial_i c_{jk} \partial_i c_{jk},
\end{equation}
where $H$, $Q$ and $\beta$ are some material-dependent parameters and the last term expresses weakly non-local effects as in \cite{Vitek-stress-diffusion}. The weakly non-local effects could play a role when the conformation tensor should exhibit steep spatial variations, as is usual in shear banding. Entropy \eqref{eq.CT.S} is just one possible choice made here for being explicit. However, the concrete form of entropy depends on the class of materials considered.

Finally, Eq.~\eqref{eq.evo.c} can be rewritten as
\begin{equation}
  \pd{\mathbf{c}}{t} = - \pd{}{\mathbf{c}^*} \Psi^{(CT \rightarrow CH)},
\end{equation}
where the reducing MTL function is
\begin{equation}
  \Psi^{(CT\rightarrow CH)} = -\Xi^*(\mathbf{c}^*) + \langle \Gamma(\mathbf{u}^*), \mathbf{c}^*\rangle.
\end{equation}
This is completely analogical to Sec.~\ref{sec.fast}, and analysis of the second law of thermodynamics could proceed as in that section.

This alternative extension brings this paper closer to the common modeling of shear banding, where viscoelastic models are usually employed (see the discussion in Sec.~\ref{sec.CR}). Similarly, by choosing a different operator $\Gamma$, one can derive evolution equations for Reynolds stress or weakly non-local vorticity, that can be found in \cite{Miroslav-turbulence,PhysicaD-2015}. In the following section, however, we stick to the extension where conjugate stress tensor is the extra state variable for simplicity.

\section{Interpretation of CR-stability}
\subsection{Thermodynamic interpretation}\label{sec.stability}

The reducing MTL function $\Psi^{(EH\rightarrow CH)}$ is plotted in Fig.~\ref{fig:slope-transform} for several values of shear rates $\gradsym$.

\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.6\textwidth]{slope-transform.pdf}
  \caption{Reducing MTL function $\Psi^{(EH\rightarrow CH)}$, Eq.~\eqref{eq.Psi.EH.CH}. For low shear rates $\absnorm{\ensuremath{\tensorq{D}}}$ (thermodynamic force $\gradsym$) the function has only one local maximum, which is located near the origin ($\absnorm{{\traceless{\cstress}}} = \unit[0]{Pa}$). As the shear rate increases, a second local maximum develops in regions with higher shear stress ${\traceless{\cstress}}$. As $\absnorm{\ensuremath{\tensorq{D}}}$ further increases, the second maximum becomes higher than the first maximum. Finally, for even higher shear rates, the first maximum disappears. Local maxima are indicated by circles ($\circ$) while local minima by squares ($\square$).}
  \label{fig:slope-transform}
\end{figure}

The curves plotted in the figure correspond to the following different regimes of the reducing evolution, indicated also in Fig.~\ref{fig:conjugate-dissipation-potential}.
\begin{enumerate}[I]
\item For low shear rates $\absnorm{\ensuremath{\tensorq{D}}}$ (Regime I when comparing to Fig. \ref{fig:conjugate-dissipation-potential}) there is only one local maximum, which is near the origin ($\absnorm{{\traceless{\cstress}}} = \unit[0]{Pa}$), as the curve for $\absnorm{\ensuremath{\tensorq{D}}} = \unit[0.05]{s^{-1}}$ indicates. Bearing in mind that $\Psi^{(EH\rightarrow E)}$ has to be chosen concave, Eq.~\eqref{eq.reducing.J} tells that shear stress evolves in the direction of derivative of the reducing MTL function. In other words, if the derivative is positive, ${\traceless{\cstress}}$ increases, and if it is negative, ${\traceless{\cstress}}$ decreases. Thus, if ${\traceless{\cstress}}$ is initially to the left of the local maximum, where the derivative is positive, it tends to the local maximum. If it is to the right of the maximum, where the derivative is negative, it again tends to the maximum. ${\traceless{\cstress}}$ is thus attracted to the local (and global) maximum. The evolution towards this maximum is thus CR-stable.
\item For slightly higher shear rates (Regime II), another local maximum appears in higher shear stresses (curve for $\absnorm{\ensuremath{\tensorq{D}}} = \unit[0.12]{s^{-1}}$). Similarly as in Regime I (low shear rates), this new local maximum attracts ${\traceless{\cstress}}$ from its neighborhood. However, since the first maximum is still higher, settling of ${\traceless{\cstress}}$ in the second maximum can be seen as a CR-metastable state. After a sufficient fluctuation, ${\traceless{\cstress}}$ will end up in the first (higher) maximum. Moreover, a local minimum appears between the local maxima. The local minimum is an CR-unstable state, since ${\traceless{\cstress}}$ is attracted to the local maxima regardless how close to the local minimum it is. That can be seen from signs of derivative of the reducing MTL function to the left and to the right of the local minimum.
\item For even higher shear rates (Regime III and curve $\absnorm{\ensuremath{\tensorq{D}}} = \unit[0.2]{s^{-1}}$), the second maximum becomes higher than the first one, and the first maximum thus becomes CR-metastable while the second becomes CR-stable. The local minimum is still located between the maxima, and although it is a stationary point ($\dot{{\tensorq{T}}}_\delta = \tensorq{0}$ in the local minimum), such a state is CR-unstable, since any arbitrarily small fluctuation in the stress would make the stress go towards one of the maxima. Note that the value of the shear rate between Regimes II and III, when the two local maxima have the same height, corresponds to the shear rate selected by the Maxwell lever rule (see Sec. \ref{sec.CR}).
\item Finally, (Regime IV) the first local maximum disappears (as well as the local minimum) for even higher shear rates (curve $\absnorm{\ensuremath{\tensorq{D}}} = \unit[0.35]{s^{-1}}$), and the second maximum becomes the only CR-stable state.
\end{enumerate}

The described dynamics of the stress tensor can be seen also in Fig.~\ref{fig:conjugate-dissipation-potential}. Regime I (low shear rates) corresponds to the part of curve $A-B$ left of point $D$, where only one ${\traceless{\cstress}} = \Xi_\ensuremath{\tensorq{D}}$ is possible. Regime II (one CR-stable state, one CR-unstable state and one CR-metastable state) is indicated in that figure as well. The upper solid curve is the CR-stable state, the lower solid curve is the CR-metastable state and the dashed curve is the CR-unstable state. Regime III (one CR-metastable state, one CR-unstable state and one CR-stable state) is described analogically. Finally, Regime IV has only one state, which is CR-stable. The reducing evolution \eqref{eq.reducing.J} is thus able to distinguish between CR-stable, CR-metastable and CR-instable states in Fig.~\ref{fig:conjugate-dissipation-potential}, and the analysis is similar to the treatment of phase transitions in equilibrium thermodynamics \cite{Landau5,Callen}. We are thus coming to a kind of \textbf{dissipative phase transitions}.

Let us now return to the experimental results in Fig.~\ref{fig:conjugate-dissipation-potential} and to the second question raised in the introduction. When slowly raising the shear rate, we start in Regime I, where shear stress is determined uniquely. After passing to Regime II, an another possible state appears, but shear stress evolves continuously due to carefulness of the experimentalists (smallness of fluctuations). Passing also through Regime III, point $C$ (boundary between Regimes III and IV) is reached, where the CR-metastable state becomes CR-unstable and shear stress jumps to the CR-stable state. That is why \textbf{shear stress behaves discontinuously}.

Moreover, it could be expected that when going back from high shear rates to the small ones, shear stress would fall as low as to point $D$ in Fig.~\ref{fig.exp}, from which it would jump to the CR-stable state on branch $A$. One should thus observe \textbf{hysteresis} indicated in that figure. We are not aware of such a backward experiment having been carried out, but modeling of hysteresis by means of introducing a new evolution equation is usual in literature \cite{kubin.lp.poirier.jp:relaxation}.

On the other hand, when slowly varying ${\traceless{\cstress}}$, there is always a unique value of $\ensuremath{\tensorq{D}}$, determined by slope of dissipation potential $\Xi^*$ in Fig.~\ref{fig:dissipation-potential}, and the slope varies continuously. That is why \textbf{no jump in shear rate} is observed in the experiment.

  \subsection{Stability of the regularized constitutive relation}
  \label{sec.stab-CRstab}

  The same dynamical behavior can be concluded directly from the stability analysis of the regularized constitutive relation \eqref{eq:26b}. Its one-dimensional counterpart with a general conjugate dissipation potential $\Xi^*$ reads
  \begin{equation}
    \label{eq:45}
    \frac{1}{\beta} \pd{\absnorm{{\traceless{\cstress}}}}{t} = \absnorm{\ensuremath{\tensorq{D}}} - \Xi^*_{\absnorm{{\traceless{\cstress}}}} \left(\absnorm{{\traceless{\cstress}}} \right).
  \end{equation}

In the shear-rate-controlled experiment we assume that $\absnorm{\ensuremath{\tensorq{D}}} = \absnorm{\ensuremath{\tensorq{D}}} (t)$ is given. Equation \eqref{eq:45} is then a non-autonomous ordinary differential equation with stationary points $\absnorm{\ensuremath{\tensorq{D}}} - \Xi^*_{\absnorm{{\traceless{\cstress}}}} \left(\absnorm{{\traceless{\cstress}}} \right)$. It can equivalently written as an extended autonomous system
  \begin{subequations}
    \begin{align}
      \pd{\absnorm{{\traceless{\cstress}}}}{t} &= \beta \left[ \absnorm{\ensuremath{\tensorq{D}}} (\tau) - \Xi^*_{\absnorm{{\traceless{\cstress}}}} \left(\absnorm{{\traceless{\cstress}}} \right) \right], \\
      \pd{\tau}{t} &= 1.
    \end{align}
  \end{subequations}
  From this system, we can deduce the phase field in the $\absnorm{\ensuremath{\tensorq{D}}}$--$\absnorm{{\traceless{\cstress}}}$ plane. In Figure~\ref{fig:phase-field}, we show the phase field for the particular choice of the dual dissipation potential \eqref{eq.theta}.
  \begin{figure}[!htbp]
    \centering
    \includegraphics[width=0.8\textwidth]{phase-field.pdf}
    \caption{Constitutive relation generated by the dual dissipation potential \eqref{eq.theta} with parameters $\alpha = 1$, $\beta = 1$, $\gamma = 0.12$ and $n = -2$, and the corresponding phase field. Vectors are rescaled to have uniform size.}
    \label{fig:phase-field}
  \end{figure}
  We can easily see that the solutions are quickly (depending on the parameter $\beta$) attracted towards the increasing parts of the constitutive relation and repelled from the decreasing part. In other words, the decreasing part is unstable in the usual mathematical sense while the increasing parts are asymptotically stable. This leads to the jumps (and possibly hysteresis) in the values of the stress $\absnorm{{\traceless{\cstress}}}$ when changing $\absnorm{\ensuremath{\tensorq{D}}}$.

 On the other hand, if we prescribe $\absnorm{{\traceless{\cstress}}} = \absnorm{{\traceless{\cstress}}} (t)$ (shear stress controlled experiment) and we further assume that its time derivative is small enough, \eqref{eq:45} reduces back to the original constitutive relation. In this case, there are no jumps.

At least in this case the CR-stability analysis leads to the same conclusion as the usual stability analysis for ordinary differential equations provided the evolution of the extra state variable (here conjugate stress) is assumed much faster than evolution of the remaining variables.

\subsection{Linearized hydrodynamic stability}
\label{sec:linearized-stability}

The CR-instability of the decreasing part of the constitutive relation can be also verified by means of the linearized hydrodynamic stability. Perturbing a known simple shear base flow, the aim is to study the temporal evolution of the disturbance. Linearizing the governing equations for the disturbance yields a generalized eigenvalue problem. Since constitutive relation \eqref{eq:le-roux-rajagopal} does not have to be invertible, the resulting problem might be in the form of a system of equations, unlike its classical counterpart---the single Orr--Sommerfeld equation, see Orr \cite{orr.w:stability*1,orr.w:stability} and Sommerfeld \cite{sommerfeld.a:beitrag}. Using a variant of a pseudospectral collocation method, e.g. Canuto and Hussaini \cite{canuto.c.hussaini.my.ea:spectral*1,canuto.c.hussaini.my.ea:spectral}, we are able to solve the problem and numerically confirm that once the flow is governed by a non-monotone constitutive relation, the perturbation grows in time and the flow is unstable. The derivation of the problem and description of the numerical method is rather lengthy and it will be subject of a future publication.

The hydrodynamic stability is usually studied for a Poiseuille or Couette flow of a fluid in some archetypal geometry like a two-dimensional channel or a pipe with prescribed boundary conditions. In the classical setting, this yields a single differential eigenvalue equation, while in our case, it leads to a system of equations. In either case, it has to be solved numerically and the spectral collocation method has proven effective, see \cite{Orszag}. This is different from, for example, Ván \cite{Van-stability} who could use the Routh--Hurwitz criteria, since he studied the stability of the equilibrium and he considered the perturbation in the form of an unbounded plane wave, resulting in a third-order polynomial. We cannot make use of the mentioned criteria as we want to study stability of non-equilibrium flows with imposed perturbations also satisfying the boundary conditions, and consequently we do not obtain just a polynomial. Anyway, determining the stability of the equilibrium can be also done from the kinetic energy of the disturbance.

\section{Braun--Le Chatelier principle}

Braun--Le Chatelier principle in equilibrium thermodynamics is the statement that a perturbation of a thermodynamic force acting on a system triggers such internal processes inside the system that the effect of the perturbation is reduced after the processes have relaxed \cite{LeChatelier}. The principle can be derived from convexity of the energy functional \cite{Landau5}, and it is equivalent with stability of the system. Analogical derivation starting from convexity of a dissipation potential leads to an analogical result, Braun--Le Chatelier principle for dissipative thermodynamics \cite{Braun-LeChat},
\begin{equation}
  \left( \pd{\absnorm{{\traceless{\cstress}}}}{\absnorm{\ensuremath{\tensorq{D}}}} \right)_{J_{int}=0} <
  \left( \pd{\absnorm{{\traceless{\cstress}}}}{\absnorm{\ensuremath{\tensorq{D}}}} \right)_{X_{int}},
\end{equation}
where $J_{int}$ is a rate of change of an internal variable and $X_{int}$ is the thermodynamic force corresponding on the internal variable. In other words, apparent viscosity (the derivative of stress with respect to strain) is lower after the internal variable has relaxed ($J_{int}=0$) than just after the perturbation (change of shear rate).

Therefore, Braun--Le Chatelier principle, which is a consequence of relaxation of an internal variable and convexity of a dissipation potential, leads to shear thinning. This is in agreement with that shear thinning is usually observed in the gradient banding experiments, which are well described by a viscoelastic model. The internal variable dynamics of which leads to the viscoelastic behavior relaxes and causes shear thinning in agreement with the Braun--Le Chatelier principle.

On the other hand, shear thickening is more difficult to explain, which supports the use of non-linear or even non-convex dissipation potentials.

\section{Non-smooth dissipation potentials}

This section stands a little apart from the rest of the paper, as it addresses yield stress behavior instead of S-curves. It is included because the authors believe that discussing non-smoothness of dissipation potentials is a natural continuation of the discussion of non-convexity.

Inspired by~\cite{Rajagopal2004}, dissipation potential exhibiting yield stress (no shear rate for sufficiently low shear stress) is for example
\begin{equation}
  \label{eq.Bingham.Xi}
  \Xi(\ensuremath{\tensorq{D}}) = \int \mathrm{d}\mathbf{r}\, \alpha \absnorm{\ensuremath{\tensorq{D}}} + \frac{1}{2}\beta \absnorm{\ensuremath{\tensorq{D}}}^2,
\end{equation}
where $\alpha$ and $\beta$ are positive constants, see Fig.~\ref{fig.Bingham.Xi}. The yield stress (given by the coefficient $\alpha$) is caused by the fact that for low shear stresses (low slopes of the dissipation potential) there is no corresponding $\ensuremath{\tensorq{D}}$. More exactly, the only shear rate corresponding to low shear stresses is zero.

\begin{figure}[!htbp]
  \centering
  \begin{minipage}[t]{0.48\textwidth}
  \includegraphics[width=\textwidth]{bingham-dissipation-potential.pdf}
  \caption{Dissipation potential \eqref{eq.Bingham.Xi} exhibiting yield stress. The potential is not differentiable at $\absnorm{\ensuremath{\tensorq{D}}} = 0$. Numerical values of $\alpha$ and $\beta$ are taken as 1.}
  \label{fig.Bingham.Xi}
  \end{minipage}
  \hfill
  \begin{minipage}[t]{0.48\textwidth}
  \includegraphics[width=\textwidth]{bingham-conjugate-dissipation-potential.pdf}
  \caption{Conjugate dissipation potential $\Xi^*$ obtained by Legendre transformation using subdifferentials from~\eqref{eq.Bingham.Xi}. Note that it is smooth everywhere in contrast to $\Xi$. Conjugate potential is flat near the origin and grows for higher shear stresses.}
  \label{fig.Bingham.Xic}
  \end{minipage}
\end{figure}

Dissipation potential \eqref{eq.Bingham.Xi} is not smooth at $\absnorm{\ensuremath{\tensorq{D}}} = 0$, where it has no derivative. Therefore, Legendre transformation to $\Xi$ has to be generalized by means of subdifferentials, see e.g. \cite{Roubicek},
\begin{equation}
  \label{eq.LT.subdif}
  {\traceless{\cstress}} \in \partial \Xi,
\end{equation}
where $\partial \Xi$ is the set of slopes of all hyperplanes touching $\Xi$ at a point $\ensuremath{\tensorq{D}}$. When constructing the subdifferential we assume that $\Xi$ is convex and continuous in a neighborhood of the point of non-smoothness. Continuity then guarantees that the set is non-empty (a consequence of Hahn--Banach theorem as commented in \cite{Roubicek}). Non-smoothness in non-convex regions can be handled as in Remark 5.8 in \cite{Roubicek}. When $\Xi$ is smooth, i.e. it has standard functional derivative, inclusion \eqref{eq.LT.subdif} becomes equality and the Legendre transformation restores its usual form.

Inclusion \eqref{eq.LT.subdif} can be interpreted as the following problem: For a given ${\traceless{\cstress}}$ find all values of $\ensuremath{\tensorq{D}}$ such that the inclusion is fulfilled. Solution to the problem is a graph $\ensuremath{\tensorq{D}}({\traceless{\cstress}})$. It can be seen in Fig.~\ref{fig.Bingham.Xi} that for sufficiently high shear stresses, there is always a non-zero $\ensuremath{\tensorq{D}}$ for which $\Xi_{\ensuremath{\tensorq{D}}}={\traceless{\cstress}}$. However, only $\partial \Xi$ at $\absnorm{\ensuremath{\tensorq{D}}} = 0$ contains the low shear stresses, which means that when the shear stress is sufficiently low, only $\absnorm{\ensuremath{\tensorq{D}}} = 0$ is the solution to the problem. In other words, there is no shear rate for sufficiently low shear stresses, but when shear stress crosses a threshold, shear rates starts to increase. We will return to this point graphically later on.

The conjugate dissipation potential is defined as
\begin{equation}
  \Xi^*({\traceless{\cstress}}) = -\Xi \left( \ensuremath{\tensorq{D}}({\traceless{\cstress}}) \right) + \int\mathrm{d}\mathbf{r}\, \ensuremath{\tensorq{D}}({\traceless{\cstress}}):{\traceless{\cstress}},
\end{equation}
and it is depicted in Fig.~\ref{fig.Bingham.Xic}. Since the conjugate dissipation potential is smooth, the backward Legendre transformation is carried out with derivatives as usually.

Let us now have a look at these dissipation potentials from the perspective of multiscale mesoscopic thermodynamics as in Sec. \ref{sec.CR-thermodynamics}. The reducing MTL function is again $-\Xi^*({\traceless{\cstress}}) + \int\mathrm{d}\mathbf{r}\, {\traceless{\cstress}} :\ensuremath{\tensorq{D}}$ and it is plotted in Fig.~\ref{fig.Bingham.Massieu}.
\begin{figure}[!htbp]
  \centering
  \includegraphics[width=0.6\textwidth]{bingham-slope-transform.pdf}
  \caption{Reducing MTL function for the yield stress model with parameters $\alpha=\beta=1$. Shear stress is attracted to the local maxima. For $\absnorm{\ensuremath{\tensorq{D}}} = 0$ the region near $\absnorm{{\traceless{\cstress}}} = 0$ is flat, which means that the shear stress can end up anywhere between $0$ and $1$. For higher shear rates a local maximum develops for non-zero stress.}
  \label{fig.Bingham.Massieu}
\end{figure}
The shear stress again tends to the maximum of the MTL function for the given shear rate. When $\absnorm{\ensuremath{\tensorq{D}}} = 0$, norm of shear stress can end up anywhere between 0 and 1 as the function is flat in that region. When shear rate is higher, the maximum is located to the right of $\absnorm{{\traceless{\cstress}}} = 1$. In other words, if there is non-zero shear rate, shear stress has to be higher than the threshold $\absnorm{{\traceless{\cstress}}} = 1$. The material thus exhibits yield stress behavior and its irreversible evolution is CR-stable.

  \section{Critical heat flux}
  \label{sec.CHF}

  In 1934, Nukiyama \cite{nukiyama.s:maximum} experimentally showed that by controlling the heat flux from a hot platinum wire to boiling water, the temperature of the wire jumps between the heating branch and the cooling branch and exhibits a hysteresis loop. Then in 1937, Drew and Mueller \cite{drew.tb.mueller.ac:boiling} were able to control the temperature and observed that the heat flux initially grows with the temperature, then decreases until it finally grows again. The boiling curve, heat flux versus temperature difference (thermodynamic force) dependence, is then non-monotonous, see also Stosic \cite{Stosic} and the book by Bergman and co-authors \cite{bergman.tl:fundamentals}. In Fig.~\ref{fig:boiling-curve}, we have plotted the experimental data by Nukiyama \cite{nukiyama.s:maximum} and fitted them with a relation of the form \eqref{eq:le-roux-rajagopal}.

  \begin{figure}[!htbp]
    \centering
    \includegraphics[width=0.6\textwidth]{nukiyama.pdf}
    \caption{Boiling curve for a platinum wire of diameter $d = \unit[0.14]{mm}$ and boiling water $T = \unit[100]{^{\circ} C}$ from Nukiyama \cite{nukiyama.s:maximum}. The data were fitted with a relation of the form \eqref{eq:le-roux-rajagopal} with $\alpha = 0.98$, $\beta = 0.00026$, $\gamma = 0.03$ and $n = -1.41$.}
    \label{fig:boiling-curve}
  \end{figure}

  As the flux-force relation can be described with an analogue of \eqref{eq:le-roux-rajagopal}, it can be likewise generated by a non-convex dissipation potential. Such a formulation would then lead to implications on the CR-stability of the constitutive relation and the experimentally observed hysteresis. For instance, the CR-unstable part is inaccessible when controlling the heat flux while careful temperature variation makes it available. This is analogous to the stress/shear rate behavior described in this paper. We hope to address this issue more explicitly in the near future.

\section{Conclusion}
\label{sec:conclusion}

Constitutive relation~\eqref{eq:le-roux-rajagopal}, that fits well the considered experimental data (Fig.~\ref{fig.exp}), can be formulated by means of dissipation potential $\Xi^*({\traceless{\cstress}})$, defined in Eq.~\eqref{eq.theta}. This dissipation potential is convex near equilibrium, but loses convexity for higher shear stresses. Due to the loss of convexity the Legendre-conjugate dissipation potential $\Xi(\ensuremath{\tensorq{D}})$ is multivalued. This reflects the experimental observation that varying shear stress leads to continuous evolution of shear rate (value of $\ensuremath{\tensorq{D}} = \Xi^*_{{\traceless{\cstress}}}$ is unique) while varying shear rate leads to a jump (and possibly hysteresis) in shear stress (value of ${\traceless{\cstress}} = \Xi_{\ensuremath{\tensorq{D}}}$ is not unique). The phenomenon of critical heat flux can be addressed analogically.

By extending classical hydrodynamics to extended hydrodynamics as in Sec.~\ref{sec.CR-thermodynamics}, one can show which parts of the multivalued dissipation potential $\Xi$ are CR-stable (with respect to perturbations of the constitutive relation), which are CR-metastable and which CR-unstable. It is then possible to identify between which branches of the dissipation potential the shear stress will jump, see Sec.~\ref{sec.stability}. The findings are compatible with the experimental data.

\section*{Acknowledgment}
We are grateful to Vít Průša, Václav Klika and On\v{r}ej Sou\v{c}ek for their invaluable comments when discussing the manuscript and for pointing out further research directions (e.g. hysteresis modeling in rheology). We also want to express our gratitude to Gian Paolo Beretta for drawing our attention to the Critical heat flux.

Most ideas in this work are either inspired or brought up by Miroslav Grmela, founder of our thermodynamic tribe.

This work was supported by Czech Science Foundation, project no.  17-15498Y.
This work has been supported by Charles University Research program No.
UNCE/SCI/023.

Adam Jane\v{c}ka acknowledges the support of Project No. LL1202 in the programme ERC-CZ funded by the Ministry of Education, Youth and Sports of the Czech Republic and the support of Project 260 449/2017 ``Student research in the field of physics didactics and mathematical and computer modelling''.

\bibliographystyle{spmpsci}