Mathematical Aspects in Non-equilibrium Thermodynamics

Edited by Róbert Kovács, Patrizia Rogolino and Francesco Oliveri

www.mdpi.com/journal/symmetry

## **Mathematical Aspects in Non-Equilibrium Thermodynamics**

## **Mathematical Aspects in Non-Equilibrium Thermodynamics**

Editors

**R ´obert Kov´acs Patrizia Rogolino Francesco Oliveri**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Robert Kov ´ acs ´ Wigner Research Institute for Physics Hungary

Patrizia Rogolino Universita di Messina ` Italy

Francesco Oliveri Universita di Messina ` Italy

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Symmetry* (ISSN 2073-8994) (available at: https://www.mdpi.com/journal/symmetry/special issues/Mathematical Aspects Non-equilibrium Thermodynamics).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-7762-3 (Hbk) ISBN 978-3-0365-7763-0 (PDF)**

© 2023 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editors**

#### **R ´obert Kov ´acs**

Robert Kov ´ acs has a degree in mechanical engineering, graduated in 2015, and then obtained the ´ Ph.D in 2017 in Budapest University of Technology and Economics (BME), in collaboration with the Hungarian Academy of Sciences. He is currently working on continuum thermodynamics, both on experimental and theoretical problems, including the development of numerical schemes, analytical solutions of evolution equations, and other mathematical aspects. He has been a visting researcher at Northeastern University, MA, USA; at Messina University, Italy; and at Universite du Qu ´ ebec ´ a` Chicoutimi, Quebec, Canada. ´

#### **Patrizia Rogolino**

Patrizia Rogolino is Professor of mathematical physics at the University of Messina. Her research interests include irreversible thermodynamics of continuous media, nonlinear theories of heat conduction, and thermoelectric phenomena in nanosystems. She is the author of several scientific papers on peer reviewed journals.

#### **Francesco Oliveri**

Francesco Oliveri is full professor of mathematical physics at the University of Messina. His research interests include mathematical apparatus of quantum mechanics for the modelization of macroscopic systems, Lie symmetries of differential equations, thermodynamics of continuous media with nonlocal constitutive relations, computer algebra. He is the author of numerous scientific papers on these topics in peer-reviewed journals as well author and editor of some volumes.

### *Editorial* **Mathematical Aspects in Non-Equilibrium Thermodynamics**

**Róbert Kovács 1,2,3,\*, Patrizia Rogolino <sup>4</sup> and Francesco Oliveri <sup>4</sup>**


#### **Introduction**

Prof. Csaba Asszonyi, D.Sc. (1941–2022): The present Special Issue is dedicated to the memory of our beloved, respected friend, colleague and teacher, the late Professor Csaba Asszonyi.

The research field of professor Asszonyi was continuum mechanics and irreversible thermodynamics. He played a pioneering role in establishing the thermodynamical background of rock rheology and engineering rock mechanics. He was educated as a mechanical engineer, started his career as research engineer, and then performed coordinated mining research in Hungary. Later on, he went into the industry and became a company group leader. Seventeen years ago, he established the Montavid Thermodynamic Research Group.

His thinking focused on thermodynamic concepts, connecting an application-oriented, engineering attitude with deep theoretical ideas. He developed several industrial applications of thermodynamic rheology. His contributions included the extension of linear viscoelasticity with internal variables and the unification of classical rheological bodies in a thermodynamic framework. He was the author of more than two hundred articles, dozens of patents, and ten books. He refused honours and distinctions, and only at the end of his life became the honorary president of the Society for the Unity of Science and Technology.

The Special Issue "Mathematical Aspects in Non-equilibrium Thermodynamics" consists of five original research papers. Although the current topic has a long history, there are still numerous open questions regarding the structure of evolution equations, the corresponding thermodynamically compatible initial and boundary conditions, and also their relation to experimental and practical aspects. These five papers actually cover various recent and relevant topics such as optimization, finite time thermodynamics, the role of the second law in continuum physics, multi-component mixtures, and boundary conditions. We hope that this Special Issue will be able to play a role in further progress to come in the future.

In the paper "The Role of the Second Law of Thermodynamics in Continuum Physics: A Muschik and Ehrentraut Theorem Revisited" by V. A. Cimmelli and P. Rogolino [1], the authors revisited the second law of thermodynamics and how the entropy inequality plays a crucial role in the derivation of evolution equations, also providing local and global formulations of the second law. The classical results of Muschik and Ehrentraut are reformulated in the present modern mathematical context of second law, thus highlighting a few geometric aspects as an outcome. They also emphasized that the non-equilibrium concept of temperature and entropy far from equilibrium is not necessarily identical to the one close to equilibrium, and how these notions need further investigation.

The paper "Integrability of the Multi-Species TASEP with Species-Dependent Rates", written by Eunghyun Lee [2], is related to totally asymmetrical simple exclusion processes; it

**Citation:** Kovács, R.; Rogolino, P.; Oliveri, F. Mathematical Aspects in Non-Equilibrium Thermodynamics. *Symmetry* **2023**, *15*, 929. https:// doi.org/10.3390/sym15040929

Received: 2 March 2023 Accepted: 10 April 2023 Published: 17 April 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

was shown that the Bethe ansatz method is applicable to processes in the present *N*-particle mixtures with species-dependent rates, providing transition probabilities for all possible compositions of species. Despite the limitations detailed in [2], this study serves as a basis for future investigation to see if the methods can be used to study the species inhomogeneity of other multi-species models.

In the paper "Shock Structure and Relaxation in the Multi-Component Mixture of Euler Fluids" by Madjarevic et al. [3], an important benchmark study related to shock structures is presented. Here, the authors utilized a multi-component mixture of Euler fluids, whose evolution equations possess a hyperbolic structure, originating from extended thermodynamics. The present study is concerned with a three-component mixture of polyatomic gases inheriting the kinetic theory formulation for the phenomenological coefficients. The the quantitative characteristics of the shock profiles, such as the temperature overshoot, the shock thickness, and the resulting relaxation times were investigated, thus providing a deeper insight into a complex, coupled phenomenon.

The paper "Cyclic Control Optimization Algorithm for Stirling Engines" by Raphael Paul and Karl Heinz Hoffmann [4] deals with an optimization problem related to nonequilibrium Stirling engines. The authors focused their attention on the optimization of both the power and efficiency, using an indirect iterative gradient algorithm. The problem formulation led to a particular Hamiltonian system, describing attractive and repulsive limit cycles, with periodic boundary conditions. They provided detailed insight into the problem formulation and optimization algorithm, and therefore their results are of high importance in dealing with similar optimization tasks for other thermodynamic cycles.

The last paper of the present Special Issue, titled "Recent Advances on Boundary Conditions for Equations in Nonequilibrium Thermodynamics" [5], is written by Wen-An Yong and Yizhou Zhou. They focused on linearized systems obeying the hyperbolic structure originating from extended thermodynamics and reviewed the possible (proper) boundary conditions in the light of uniform and generalized Kreiss conditions. The structural stability of the studied PDEs was also satisfied. As these conditions are strongly related to the suitability of hyperbolic equations, the present results could serve as a future basis for the comparison of various thermodynamic approaches and provide hints to extend the present formalism to nonlinear problems.

**Acknowledgments:** We would like to express our gratitude to the Editorial Board of Symmetry for they helpful attitude and also to the Authors who made this Special Issue successful.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **The Role of the Second Law of Thermodynamics in Continuum Physics: A Muschik and Ehrentraut Theorem Revisited**

**Vito Antonio Cimmelli 1,\*,† and Patrizia Rogolino 2,†**


**Abstract:** In continuum physics, constitutive equations model the material properties of physical systems. In those equations, material symmetry is taken into account by applying suitable representation theorems for symmetric and/or isotropic functions. Such mathematical representations must be in accordance with the second law of thermodynamics, which imposes that, in any thermodynamic process, the entropy production must be nonnegative. This requirement is fulfilled by assigning the constitutive equations in a form that guaranties that second law of thermodynamics is satisfied along arbitrary processes. Such an approach, in practice regards the second law of thermodynamics as a restriction on the constitutive equations, which must guarantee that any solution of the balance laws also satisfy the entropy inequality. This is a useful operative assumption, but not a consequence of general physical laws. Indeed, a different point of view, which regards the second law of thermodynamics as a restriction on the thermodynamic processes, i.e., on the solutions of the system of balance laws, is possible. This is tantamount to assuming that there are solutions of the balance laws that satisfy the entropy inequality, and solutions that do not satisfy it. In order to decide what is the correct approach, Muschik and Ehrentraut in 1996, postulated an amendment to the second law, which makes explicit the evident (but rather hidden) assumption that, in any point of the body, the entropy production is zero if, and only if, this point is a thermodynamic equilibrium. Then they proved that, given the amendment, the second law of thermodynamics is necessarily a restriction on the constitutive equations and not on the thermodynamic processes. In the present paper, we revisit their proof, lighting up some geometric aspects that were hidden in therein. Moreover, we propose an alternative formulation of the second law of thermodynamics, which incorporates the amendment. In this way we make this important result more intuitive and easily accessible to a wider audience.

**Keywords:** second law of thermodynamics; dissipation principle; state space; balance laws; entropy inequality

#### **1. Introduction**

Let *B* be a continuous body undergoing a thermomechanical transformation, whose evolution in the spacetime is ruled by the system of balance laws

$$
\mathcal{U}\_{\beta,t} + \mathcal{U}\_{\beta,j} v\_j + \Phi\_{k,k}^{\beta} = r\_{\beta}, \quad \beta = 1 \dots \omega,\tag{1}
$$

with *vj* as the components of the velocity field on *B* entering the total time derivative, Φ*β <sup>k</sup>* as the components of the flux of *Uβ*, and *r<sup>β</sup>* as the production of *U<sup>β</sup>* (for the sake of simplicity, we assume that the supplies are zero). Moreover, the symbols *f*,*<sup>t</sup>* and *f*,*<sup>j</sup>* mean the partial derivatives of function *f* with respect to time and to the spatial coordinate *xj*, *j* = 1, 2, 3, respectively.

**Citation:** Cimmelli, V.A.; Rogolino, P. The Role of the Second Law of Thermodynamics in Continuum Physics: A Muschik and Ehrentraut Theorem Revisited. *Symmetry* **2022**, *14*, 763. https://doi.org/10.3390/ sym14040763

Academic Editors: Alessandro Sarracino and Kazuharu Bamba

Received: 2 February 2022 Accepted: 1 April 2022 Published: 7 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

For instance, in classical rational thermodynamics, Equation (1) is the balance of mass, linear momentum, angular momentum, and energy [1,2],

$$
\rho\_{,t} + \rho\_{,i}v\_{i} + \rho v\_{j,\bar{j}} = 0,\tag{2}
$$

$$
\rho v\_{i,t} + \rho v\_{i,j} v\_j - T\_{i j,j} = \rho b\_{i\prime} \tag{3}
$$

$$T\_{\text{ij}} = T\_{\text{ji}} \tag{4}$$

$$
\rho \varepsilon\_{,t} + \rho \varepsilon\_{,j} \upsilon\_{j} - T\_{ij} \upsilon\_{i,j} + q\_{j,j} = \rho r\_{,} \tag{5}
$$

where *ρ* denotes the mass density, *ε* the specific internal energy, *r* the energy supply, *vi*, *qi*, and *bi* the components of velocity, heat flux, and body force, respectively, and *Tij* are the components of the Cauchy stress [1,2], while in the extended non-equilibrium thermodynamic theories, taking the fluxes as independent variables, the set of field equations includes the balance laws for the independent fluxes [3–7].

We suppose that the fields *<sup>U</sup>β*, the fluxes <sup>Φ</sup>*<sup>β</sup> <sup>k</sup>* , and the productions *r<sup>β</sup>* depend on *ω* unknown fields *<sup>z</sup>α*(*xj*, *<sup>t</sup>*) and on their spatial derivatives *<sup>z</sup>α*,*j*(*xj*, *<sup>t</sup>*). Then, suitable constitutive equations must be assigned for them. Such equations must account for material symmetry, namely, the invariance of certain physical properties with respect to a group of coordinate transformations (symmetry group). Thus, suitable representation theorems for symmetric and/or isotropic functions must be applied [8]. Through the constitutive equations, material symmetry reflects itself in the Equation (1) field, which inherit from them particular mathematical properties. Meaningful examples of such reciprocal influences are illustrated in [9], where some important equations of continuum physics are studied by the Lie symmetry analysis of differential equations. On the other hand, the mathematical representation of material symmetry cannot neglect the constraint dictated by the second law of thermodynamics, which imposes the local entropy production

$$
\sigma^{(s)} = \rho \mathbf{s}\_{\prime 1} + \rho \mathbf{s}\_{\prime \dot{\jmath}} \mathbf{v}\_{\dot{\jmath}} + J\_{k,k} - \rho (r/\vartheta), \tag{6}
$$

where *s* is the specific entropy, *Jk* are the components of the entropy flux, and *ϑ* the absolute temperature, are nonnegative, namely,

$$
\sigma^{(s)} \ge 0,\tag{7}
$$

whatever the thermodynamic process is [1,2]. Indeed, the unilateral differential constraint (7) can be interpreted either as a restriction on the solutions of Equation (1), or as a restriction on the constitutive equations that characterize the fields *<sup>U</sup>β*, <sup>Φ</sup>*<sup>β</sup> <sup>k</sup>* , *rβ*, *s*, and *Jk*. In the first case, it leads to the following assumption:

**Assumption 1.** *Among the mathematical solutions of Equation* (1)*, we must find those that are physically realizable by substituting them into the Inequality* (7) *and checking the sign.*

In the second case, instead, it leads to a different assumption, namely:

**Assumption 2.** *The constitutive equations for <sup>U</sup>β,* <sup>Φ</sup>*<sup>β</sup> <sup>k</sup> , rβ, s, and Jk must be assigned in such a way that the constraint* (7) *is satisfied along arbitrary thermodynamic processes.*

Modern constitutive theories of continuum thermodynamics are based on the second statement, which was formulated for the first time by Coleman and Noll [1] in 1963, and is universally known as the entropy (or dissipation) principle [10]. On the other hand, since the determination of solutions of Equation (1) is, in general, very difficult, the Coleman– Noll approach is the most convenient one used for determining the consequences of the second law. Two rigorous mathematical procedures for analyzing the restrictions dictated by the second law on the constitutive equations are based on such an assumption, namely, the Coleman–Noll and Liu procedures [1,11,12].

In order to investigate if the entropy principle, as formulated by Coleman and Noll, is just an operative assumption or a consequence of a general physical law, Muschik and Ehrentraut in 1996 proved their celebrated theorem in [13], which we revisit here within a geometric perspective. Since the theorem has a complex mathematical formulation, herein we limit ourselves to provide a sketch of the result, referring the reader to reference [13] for details. Their starting point is the following concept of equilibrium for a thermodynamic system.

*A thermodynamic system is said to be in equilibrium (stable or metastable), if isolation of the system from its environment does not affect the observables of the system (or in other words, if the state of the isolated system is the same as the state of the system prior to insulation)*.

Then, they define the process direction vectors as those vectors that are tangent to the curve representative of a process in the state space. Moreover, they show that in a fixed point of such a curve, the entropy production is linear in the process direction vectors *yα*. The vectors *<sup>y</sup><sup>α</sup>* are said irreversible, if *<sup>σ</sup>*(*s*)(*yα*) > 0, reversible, if *<sup>σ</sup>*(*s*)(*yα*) = 0, non-real, if *<sup>σ</sup>*(*s*)(*yα*) < 0. In particular, a vector *<sup>y</sup><sup>α</sup>* such that *<sup>σ</sup>*(*s*)(*yα*) = 0, is called the reversible process direction.

At this point, Muschik and Ehrentraut prove their fundamental.

**Proposition 1.** *If in a point of non-equilibrium of the curve representative of the process in the state space there exist both irreversible and non-real process directions, then, in that point, it is possible to obtain a reversible process direction as a linear combination of the latter.*

Such a result is unphysical, since in a non-equilibrium point the entropy production can be zero. In order to overcome such a discrepancy, Muschik and Ehrentraut add to the classical formulation of the second law, represented by the Inequality (7), the following amendment.

#### *Except in equilibria, reversible process directions in state space do not exist*.

The consequences of the amendment are severe, because it implies that, in a given point, the process directions are either all irreversible (or reversible, in the limit of quasistatic transformation), or all non-real. On the other hand, since non-real process directions do not exist in nature, we must choose the first option. As a consequence, non-real process directions that must be excluded by the second law do not exist. Moreover, since the point on the curve of the process is arbitrary, we are led to conclude that there are no non-real thermodynamic processes to be forbidden by the second law and, hence, the last necessarily restricts the constitutive equations and not the thermodynamic processes. In this way, the classical Coleman–Noll approach follows by a rigorous proof.

It is worth observing that Proposition 1 and its consequences are not aimed at giving more deep physical insight on some basic concepts of thermodynamics, such as equilibrium and reversibility, since in [13] those concepts are the classical ones. Indeed, Proposition 1 aims at answering the following question, which is fundamental from the methodological point of view: is the Coleman–Noll interpretation of the second law, as a constraint on the constitutive equations, a mathematical consequence of the basic laws of thermodynamics or an additional operative assumption that does not follow from those laws? The answer, given by Proposition 1, is that such an interpretation is a consequence of the second law, provided the amendment is postulated. In the absence of the amendment, the Coleman– Noll approach would be an arbitrary assumption and, hence, it could be questioned. The present paper is motivated by the observation that the important result illustrated above can be put in a more general and accessible form within a geometric framework. To achieve that task, we use the results in references [14,15], where a geometric perspective on non-equilibrium thermodynamics has been given. The chosen state space is different with respect to that considered reference [13], because we do not include in it the time derivatives. In this way, the constitutive equations we deal with are suitable to satisfy the

principle of material indifference [2]. After defining the space of the higher derivatives, we introduce the definitions of the real, ideal, and over-ideal vector of the higher derivatives. For thermodynamic processes, we give the definitions of irreversible, reversible, and over-reversible process, by analyzing the properties of its representative curve in the fiber bundle of the configuration spaces. Once the geometric framework is complete, we reformulate the second law of thermodynamics, both locally and globally, in space and time, in order to encompass the amendment. In this way, we are able to prove a new formulation of the Muschik and Ehrentraut theorem. In the discussion in Section 5, we highlight some of the main advantages of the geometric approach presented here. We underline that our generalized formulation of the second law seems to be more adapt for some recent thermodynamic theories that analyze real transformations, which occur in a finite time, and not by quasi-static transformations, which require an infinite time to be realized [16–18]. Moreover, we will see that the use of the fiber bundle allows to provide, in a very intuitive manner, the mathematical definition of reversible and irreversible processes, and to generalize, in a natural way, local results to global ones. Finally, our approach leads to a physically sound interpretation of the principle of local equilibrium [19], which provides the theoretical justification of the definition of entropy and temperature outside the equilibrium.

The paper runs as follows.

In Section 2, we construct a new thermodynamic framework for non-equilibrium processes. In Section 3, we present a new formulation, both locally and globally, of the second law of thermodynamics. In Section 4, we prove the Muschik and Ehrentraut theorem. In Section 5, we discuss our results and present some open problems that will be considered in future research.

#### **2. The Thermodynamic Framework**

In this section, we construct a geometric framework where our main results can be formulated. To this end, we begin by providing some basic definitions.

**Definition 1.** *The space Ct of the configurations at the instant t is represented by a ω-dimensional vector space spanned by the solutions <sup>z</sup>α*(*xj*, *<sup>t</sup>*) *of Equation* (1) *with a structure of a finitedimensional manifold.*

We assume that the total configuration space is given by the disjointed union

$$\mathcal{C} = \bigcup\_{t \in [0, \infty]} \{t\} \times \mathbb{C}\_{t\prime} \tag{8}$$

with a given natural structure of a fiber bundle over the real line R where time flows [14,15].

**Definition 2.** C *is called the configuration bundle.*

Under the natural assumption that *Ct* does not vary in time, namely, *Ct* <sup>=</sup> *<sup>C</sup>* <sup>∀</sup>*t*, then C has the topology of the Cartesian product

$$
\mathcal{C} = \mathbb{R} \times \mathbb{C}.\tag{9}
$$

**Definition 3.** *A vector valued function <sup>π</sup>* : *<sup>t</sup>* <sup>∈</sup> [*τ*0, *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>*] <sup>⊆</sup> <sup>R</sup> <sup>→</sup> *<sup>z</sup>α*(*xj*, *<sup>t</sup>*) ∈ C *is said to be a thermodynamic process π of duration τ. Moreover, π* = *π*(*t*) *is the parametric equation of the curve* Γ *representative of π in* C*.*

**Definition 4.** *For <sup>t</sup>*<sup>0</sup> <sup>∈</sup> [*τ*0, *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>*]*, a vector valued function <sup>p</sup>* : *<sup>t</sup>* <sup>∈</sup> [*t*0, *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>*] <sup>⊆</sup> <sup>R</sup> <sup>→</sup> *<sup>z</sup>α*(*xj*, *<sup>t</sup>*) ∈ C *is said to be a restricted thermodynamic process <sup>p</sup> of initial point <sup>t</sup>*<sup>0</sup> *and duration <sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>* <sup>−</sup> *<sup>t</sup>*0*, reference [14]. Moreover, <sup>p</sup>* <sup>=</sup> *<sup>p</sup>*(*t*) *is the parametric equation of the curve <sup>γ</sup> representative of p in* C*.*

**Remark 1.** *For <sup>t</sup>*<sup>0</sup> = *<sup>τ</sup>*<sup>0</sup> *we have <sup>p</sup>*(*t*) = *<sup>π</sup>*(*t*)*, for <sup>t</sup>*<sup>0</sup> = *<sup>τ</sup>*<sup>0</sup> + *<sup>τ</sup>, <sup>p</sup>*(*t*) *is the process of duration* 0*, i.e., the null process.*

As said in Section 2, in order to find the fields *<sup>z</sup>α*(*xj*, *<sup>t</sup>*), i.e., to solve the system (1), for the quantities *Uβ*, Φ*<sup>α</sup> <sup>k</sup>* and *r<sup>α</sup>* constitutive equations must be assigned on a suitable state space.

**Definition 5.** *The* 4*ω-dimensional vector space with the structure of a finite-dimensional manifold*

$$\Sigma\_t = \left\{ z\_a(x\_{\dot{j}\prime}, t), \ z\_{a, \dot{j}}(x\_{\dot{j}\prime}, t) \right\}. \tag{10}$$

*for any value of the time variable t, it represents a local in the time state space and it is called state space at the instant t.*

**Definition 6.** *The disjoint union*

$$\mathcal{S} = \bigcup\_{t \in [0, \infty]} \{t\} \times \Sigma\_{\mathfrak{k}} \tag{11}$$

*with a given natural structure of a fiber bundle over the real line* R *where time flows, it represents the total configuration space and it is said to be the thermodynamic bundle.*

Again, under the natural assumption that <sup>Σ</sup>*<sup>t</sup>* does not vary in time, namely, <sup>Σ</sup>*<sup>t</sup>* <sup>=</sup> <sup>Σ</sup> <sup>∀</sup>*t*, then S has the topology of the Cartesian product

$$\mathcal{S} = \mathbb{R} \times \Sigma. \tag{12}$$

Of course,

$$\mathcal{C}\_t \subset \Sigma\_t, \quad \mathcal{C} \subset \mathcal{S}. \tag{13}$$

The balance Equation (1) on the local in time state space Σ*<sup>t</sup>* reads

$$\frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{a}} z\_{a,t} + \frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{a,j}} z\_{a,jt} + \frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{a}} z\_{a,j} v\_{j} + \frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{a,k}} z\_{a,kj} v\_{j} + \frac{\partial \Phi\_{k}^{\mathcal{S}}}{\partial z\_{a}} z\_{a,k} + \frac{\partial \Phi\_{k}^{\mathcal{S}}}{\partial z\_{a,j}} z\_{a,jk} = r\_{\beta}. \tag{14}$$

In Equation (14), we may individuate the 10*ω* higher derivatives *<sup>z</sup>α*,*t*, *<sup>z</sup>α*,*jt*, *<sup>z</sup>α*,*jk* , which are the space and time derivatives of the elements of Σ*t*.

**Definition 7.** *The local in time* 10*ω-dimensional vector space*

$$H\_t = \left\{ z\_{a,t}(\mathbf{x}\_{j},t), z\_{a,jt}(\mathbf{x}\_{j},t), z\_{a,jk}(\mathbf{x}\_{j},t) \right\},\tag{15}$$

*and the fiber bundle*

$$\mathcal{H} = \bigcup\_{t \in [0, \infty]} \{t\} \times H\_{\mathfrak{l}}.\tag{16}$$

*represent the space of the higher derivatives at time t and its fiber bundle, respectively.*

**Definition 8.** *A point P* ∈ *B is said to be in a state of equilibrium at the instant t, if*

$$z\_{\alpha,t}(P,t) + z\_{\alpha,j}(P,t)v\_j = 0 \,\,\forall \, t\_r$$

*and*

$$z\_{\alpha,jt}(P,t) + z\_{\alpha,jk}(P,t)v\_k = 0 \,\,\forall \,\, t.$$

The Definition 8 leads to

**Definition 9.** *The equilibrium subspace of Ht and its fiber bundle are given by*

$$\mathcal{E}\_t = \left\{ z\_{a,jk}(x\_{\bar{\jmath}}, t) \right\}\_{\prime} \tag{17}$$

*and*

$$\hat{\mathcal{E}} = \bigcup\_{t \in [0, \infty]} \{t\} \times \hat{E}\_t. \tag{18}$$

Analogously, the entropy inequality on the state space reads

$$
\rho \frac{\partial \mathbf{s}}{\partial z\_{\mathbf{z}}} z\_{\mathbf{z},t} + \rho \frac{\partial \mathbf{s}}{\partial z\_{\mathbf{z},j}} z\_{\mathbf{z},jt} + \rho \frac{\partial \mathbf{s}}{\partial z\_{\mathbf{z}}} z\_{\mathbf{z},j} v\_{j} + \rho \frac{\partial \mathbf{s}}{\partial z\_{\mathbf{z},k}} z\_{\mathbf{z},kj} v\_{j} + \frac{\partial f\_{k}}{\partial z\_{\mathbf{z}}} z\_{\mathbf{z},k} + \frac{\partial f\_{k}}{\partial z\_{\mathbf{z},j}} z\_{\mathbf{z},jk} \ge 0. \tag{19}
$$

**Definition 10.** *The local in time* 10*ω-dimensional vector space at time t*

$$\mathcal{W}\_t = \left\{ z\_{\mathfrak{a},t}(\mathbf{x}\_{j\prime},t), z\_{\mathfrak{a},j\dagger}(\mathbf{x}\_{j\prime},t), z\_{\mathfrak{a},jk}(\mathbf{x}\_{j\prime},t) \right\},\tag{20}$$

*and the fiber bundle*

$$\mathcal{W} = \bigcup\_{t \in [0, \infty]} \{t\} \times \mathcal{W}\_t. \tag{21}$$

*define the vector space and the fiber bundle of the higher derivatives, respectively, whose state vectors satisfy the entropy inequality. Moreover, the equilibrium subspace of Wt and its fiber bundle are given by*

$$E\_1 = \left\{ z\_{u,jk}(x\_{j'}, t) \right\}\_{\prime} \tag{22}$$

*and*

$$\mathcal{E} = \bigcup\_{t \in [0, \infty]} \{t\} \times E\_t. \tag{23}$$

**Remark 2.** *We defined two different spaces of the higher derivatives, one for the balance equations and another for the entropy inequality, related to the fundamental focus of the present investigation, namely, to determine the conditions, if any, under which all of the solutions of the balance laws are also solutions of the entropy inequality. This will be discussed in detail in the next section.*

The relations in Equations (14) and (19) can be arranged as follows

$$\frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{a}} z\_{a,t} + \frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{a,j}} z\_{a,jt} + \left(\frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{a,k}} v\_{j} + \frac{\partial \Phi^{\beta}\_{j}}{\partial z\_{a,k}}\right) z\_{a,kj} = r\_{\beta} - \frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{a}} z\_{a,j} v\_{j} - \frac{\partial \Phi^{\beta}\_{j}}{\partial z\_{a}} z\_{a,j}.\tag{24}$$

$$
\rho \frac{\partial \mathbf{s}}{\partial z\_a} z\_{a,t} + \rho \frac{\partial \mathbf{s}}{\partial z\_{a,j}} z\_{a,jt} + \left(\rho \frac{\partial \mathbf{s}}{\partial z\_{a,k}} v\_i + \frac{\partial f\_i}{\partial z\_{a,k}}\right) z\_{a,ki} \ge -\rho \frac{\partial \mathbf{s}}{\partial z\_a} z\_{a,j} v\_j - \frac{\partial f\_i}{\partial z\_a} z\_{a,i} . \tag{25}$$

Let us now define the 10*ω* × 1 column vector function

$$\mathbf{y}\_{\mathbf{a}} \equiv \begin{pmatrix} z\_{\mathbf{a}, \mathbf{t}\_{\prime}} \ z\_{\mathbf{a}, \mathbf{j}, \mathbf{t}\_{\prime}} \ z\_{\mathbf{a}, \mathbf{k} \mathbf{j}} \end{pmatrix}^{T}, \tag{2.6}$$

the *ω* × 1 column vector

$$\mathbf{C}\_{\beta} \equiv r\_{\beta} - \frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{a}} z\_{a,j} \upsilon\_{j} - \frac{\partial \Phi\_{j}^{\beta}}{\partial z\_{a}} z\_{a,j}, \quad \beta = 1 \ldots \omega,\tag{27}$$

and the *ω* × 10*ω* matrix

$$A\_{\beta\mathfrak{a}} \equiv \left[ \frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{\mathfrak{a}}}, \ \ \ \frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{\mathfrak{a},j}}, \ \left( \ \frac{\partial \mathcal{U}\_{\beta}}{\partial z\_{\mathfrak{a},k}} v\_{\mathfrak{f}} + \frac{\partial \Phi\_{j}^{\beta}}{\partial z\_{\mathfrak{a},k}} \right) \right], \quad j, k = 1, 2, 3,\tag{28}$$

with *C<sup>β</sup>* and *Aβα* defined on S. In this way, the balance Equation (24) can be rearranged as

$$A\_{\beta\mathfrak{a}}(\mathcal{S})y\_{\mathfrak{a}} = \mathbb{C}\_{\beta}(\mathcal{S}).\tag{29}$$

Analogously, after defining the 10*ω* × 1 column vector function

$$B\_{\mathbf{a}}(\mathcal{S}) \equiv \left(\rho \frac{\partial \mathbf{s}}{\partial z\_{a}}, \ \ \rho \frac{\partial \mathbf{s}}{\partial z\_{a,j}}, \ \left(\rho \frac{\partial \mathbf{s}}{\partial z\_{a,k}} v\_{i} + \frac{\partial f\_{i}}{\partial z\_{a,k}}\right)\right)^{T},\tag{30}$$

and the scalar function

$$D(\mathcal{S}) \equiv -\rho \frac{\partial \mathbf{s}}{\partial z\_{\alpha}} z\_{\alpha,j} v\_j - \frac{\partial f\_i}{\partial z\_{\alpha}} z\_{\alpha,i} \tag{31}$$

we can write the Inequality (25) as

$$B\_{\mathfrak{A}}(\mathcal{S})y\_{\mathfrak{A}} \ge D(\mathcal{S}).\tag{32}$$

**Remark 3.** *It is worth observing that the higher derivatives entering the system* (29) *are elements of Ht, while those entering the Inequality* (32) *are elements of Wt.*

From now on, we pursue our analysis under the hypothesis that *B* occupies the whole space. Then, for arbitrary *<sup>t</sup>*<sup>0</sup> <sup>∈</sup> [*τ*0, *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>*], we consider the restricted process *<sup>p</sup>* of the initial instant *<sup>t</sup>*<sup>0</sup> and duration *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>* <sup>−</sup> *<sup>t</sup>*0, and suppose that it corresponds to the solution of the Cauchy problem for the system (29) with initial conditions

$$z\_{\mathfrak{a}}(\mathbf{x}\_{\mathfrak{j}}, \mathfrak{t}\_{\mathfrak{l}}) = z\_{\mathfrak{a}} \mathfrak{o}(\mathfrak{x}\_{\mathfrak{j}}), \quad \forall P \in \mathbb{C}. \tag{33}$$

If *Aβα* and *C<sup>β</sup>* are regular, and *Aβα* is invertible, the theorem of Cauchy–Kovalevskaya ensures that the Cauchy problem (29) and (33) has a unique solution continuously depending on the initial data (33), reference [20]. However, such a solution does not necessarily correspond to a thermodynamic process that is physically realizable, since the physically admissible solutions of (29) and (33) are only those solutions, which additionally satisfy the unilateral differential constraint (32). On the other hand, the problems (29) and (33) are very difficult to solve, in general, so that finding a solution, and verifying ex post if it also satisfies (32), does not seem to be a convenient procedure. For that reason, Coleman and Noll [1] in 1963 postulated the constitutive principle referred in Section 2, reference [10]. It is important to investigate if the Coleman and Noll postulate is a consequence of a general physical law or an arbitrary, although very useful, assumption, as observed by Muschik and Ehrentraut [13]. Such a study will be carried out in the following sections.

#### **3. Local and Global Formulations of the Second Law of Thermodynamics**

Let us consider now a fixed point *P*<sup>0</sup> ∈ *B* whose vector position will be indicated by **<sup>x</sup>**0, a fixed instant of time *<sup>t</sup>*<sup>0</sup> <sup>∈</sup> [*τ*0, *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>*]. We note that, whatever is *<sup>t</sup>*0, it can can ever be considered as the initial time of a restricted process of duration *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>* <sup>−</sup> *<sup>t</sup>*0. Moreover, let <sup>Σ</sup>0, *<sup>H</sup>*0, and *<sup>E</sup>*<sup>0</sup> the vector spaces <sup>Σ</sup>*t*(*P*0, *<sup>t</sup>*0), *Ht*(*P*0, *<sup>t</sup>*0), and *Et*(*P*0, *<sup>t</sup>*0). When evaluated in (*P*0, *<sup>t</sup>*0), the balance Equation (29) and the entropy Inequality (32) transform in the algebraic relations:

$$A\_{\beta\mathfrak{a}}(\Sigma\_0)y\_{\mathfrak{a}} = \mathbb{C}\_{\beta}(\Sigma\_0),\tag{34}$$

$$B\_{\mathfrak{a}}(\Sigma\_0)y\_{\mathfrak{a}} \ge D(\Sigma\_0). \tag{35}$$

In this way, we can regard the *<sup>ω</sup>* <sup>×</sup> <sup>10</sup>*<sup>ω</sup>* matrix *<sup>A</sup>βα*(Σ0) as a linear morphism from *H*<sup>0</sup> to the *ω*-dimensional Euclidean vector space defined on Σ0. Analogously, the vector *<sup>B</sup>α*(Σ0) can be regarded as a linear application from *<sup>H</sup>*<sup>0</sup> in R, so that *<sup>B</sup>α*(Σ0) belongs to the dual space *H*∗ <sup>0</sup> of *H*0. It is worth observing that, since *Aβα* is supposed to be invertible (otherwise the Cauchy problem (29) and (33) would not admit a unique solution), the algebraic relations (34) allow to determine *ω* of the 10*ω* components of *yα*. Moreover, by spatial derivation of the initial conditions (33), we have

$$z\_{\mathfrak{a},jk}(\mathfrak{x}\_{\mathfrak{j}},t\_0) = z\_{\mathfrak{a},0,jk}(\mathfrak{x}\_{\mathfrak{j}}),\tag{36}$$

which, once evaluated in *P*0, allow to determine 6*ω* components of *yα*. It is worth observing that, since the initial conditions can be assigned arbitrarily, such 6*ω* quantities can assume arbitrary values. Moreover, there are further 3*ω* components of the vector *y<sup>α</sup>* that remain completely arbitrary, since the system (34) and the initial relations (36) allow to determine only 7*ω* of the 10*ω* components of *yα*. It is not guaranteed that the Inequality (35) is satisfied whatever is *y<sup>α</sup>* ∈ *H*0. Thus, we define the space *W*<sup>0</sup> ⊆ *H*<sup>0</sup> constituted by the vectors of *H*0, which satisfy both Equation (34) and the Inequality (35).

**Remark 4.** *It is worth observing that, although it is not guaranteed that the Inequality* (35) *is satisfied, whatever is y<sup>α</sup>* ∈ *H*0*, at this stage, we do not have elements to exclude such a possibility. In other words, we do not have elements to decide if, actually, W*<sup>0</sup> *is a proper subspace of H*<sup>0</sup> *or if it coincides with H*0*.*

In order to decide if *<sup>W</sup>*<sup>0</sup> <sup>⊂</sup> *<sup>H</sup>*0, or *<sup>W</sup>*<sup>0</sup> <sup>=</sup> *<sup>H</sup>*0, we follow the way paved by Muschik and Ehrentraut [13], who observed that such a decision cannot be ensued solely by the second law of thermodynamics, because such a law does not contain information regarding Equation (34) or the initial conditions (36). In order to fill this gap, Muschik and Ehrentraut completed the information contained in the Inequality (35) by an amendment that clarifies how the reversible transformations can be realized from the operative point of view. Here, we followed their strategy, but we propose a more general approach that includes the amendment into a new formulation of the second law. To achieve that task, we needed some preliminary definitions. To this end, we observed that, in the real world, reversible thermodynamic transformations do not exist, but they are approximated by very slow (quasi-static) transformations in which, in any point *P*<sup>0</sup> ∈ *B*, the system is close to the thermodynamic equilibrium. From an ideal point of view, a quasi-static transformation requires an infinite time to occur, and in any point of the system, the value of the state variable is constant in time.

**Remark 5.** *Alternative formulations of the thermodynamic laws that consider realistic transformations occurring in a finite time are proposed within the framework of finite time thermodynamics [16–18].*

As far as the thermodynamic framework is concerned, if *B* undergoes a quasi-static transformation, along with Muschik and Ehrentraut [13], we say that, in any point (*P*0, *<sup>t</sup>*0), the vectors of the higher derivatives are elements of *E*0. Such an observation suggests the following definitions.

**Definition 11.** *A vector y<sup>α</sup>* ∈ *H*<sup>0</sup> *is said:*


We can establish the following, owing to the definitions above:

**Postulate 1** (**A local formulation of the second law of thermodynamics).** *Let B be a body, and let the couple* (*P*0, *<sup>t</sup>*0) *represent an arbitrary point of <sup>B</sup> at an arbitrary instant <sup>t</sup>*<sup>0</sup> <sup>∈</sup> [*τ*0, *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>*]*. Suppose B undergoes an arbitrary thermodynamic process of initial instant t*<sup>0</sup> *and duration <sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>* <sup>−</sup> *<sup>t</sup>*0*. Then, the local space of the higher derivatives <sup>W</sup>*<sup>0</sup> *does not contain over-ideal vectors. Moreover, a vector y<sup>α</sup>* <sup>∈</sup> *<sup>W</sup>*<sup>0</sup> *is ideal if, and only if,* (*P*0, *<sup>t</sup>*0) *is in thermodynamic equilibrium.*

The postulate, the above traduces the experimental evidence, in that, in a thermodynamic process, the entropy production cannot be negative at any point *P*<sup>0</sup> of *B* at any instant *t*0. Moreover, it also expresses the further experimental fact, which is often tacit in the formulations of the second law of thermodynamics—that the entropy production can be zero only in the points of *B* that are in equilibrium. In particular, if the point *P*<sup>0</sup> ∈ *B* at the instant *<sup>t</sup>*<sup>0</sup> is in a thermodynamic equilibrium, then *<sup>W</sup>*<sup>0</sup> <sup>=</sup> *<sup>E</sup>*0. Hence, *<sup>y</sup><sup>α</sup>* <sup>∈</sup> *<sup>E</sup>*<sup>0</sup> is ideal if, and only if, (*P*0, *<sup>t</sup>*0) is in a thermodynamic equilibrium. Finally, a vector *<sup>y</sup><sup>α</sup>* <sup>∈</sup> *<sup>W</sup>*<sup>0</sup> is either real or over-ideal if, and only if, (*P*0, *<sup>t</sup>*0) is not in a thermodynamic equilibrium.

**Remark 6.** *We should note that the local formulation of the second law of thermodynamics prohibits that over-ideal vectors be in W*0*, but it does not prevent that they be in H*0*. Whether H*<sup>0</sup> *contains over-ideal vectors (or not) is the focus of the present investigation.*

**Definition 12.** *Let B be a body undergoing an arbitrary thermodynamic process p of the initial instant <sup>t</sup>*<sup>0</sup> *and duration <sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>* <sup>−</sup> *<sup>t</sup>*0*, and let <sup>γ</sup> be the curve representative of the process in* <sup>C</sup>*. The process p is said to be:*


The definitions above allow us to enunciate the following:

**Postulate 2** (**The global formulation of the second law of thermodynamics).** *Over-reversible processes do not occur in nature. Moreover, a thermodynamic process is reversible if, and only if, any point P* ∈ *B, at any instant t, is in a thermodynamic equilibrium.*

The previous formulations (local and global) of the second law of thermodynamics include information not present in the classical ones—that the reversible transformations are necessarily quasi-static and, hence, they need an infinite time to occur. Thus, they represent ideal processes, which, in nature, are approximated by very slow transformations. Here, we take into account such a situation by admitting that this is possible if, and only if, at any point of the body, at any instant *t*, there is a state of thermodynamic equilibrium. On the other hand, this implies that at any point of *γ*, the vector of the higher derivatives *<sup>y</sup>α*(*P*, *<sup>t</sup>*) lie in the equilibrium bundle <sup>E</sup>, and is ideal.

#### **4. The Muschik and Ehrentraut Theorem Revisited**

In this section, we present a novel formulation of the Muschik and Ehrentraut theorem proved in reference [13]. To this end, we use the thermodynamic framework and the generalized formulations of the second law established above.

**Theorem 1.** *Let <sup>B</sup> be a body, and let the couple* (*P*0, *<sup>t</sup>*0) *represent an arbitrary point of <sup>B</sup> at an arbitrary instant t*<sup>0</sup> <sup>∈</sup> [*τ*0, *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>*]*. Then, H*<sup>0</sup> *= W*0*.*

**Proof.** To prove the theorem, it is enough to demonstrate that the vectors of *H*<sup>0</sup> are all (and only) the vectors of *<sup>W</sup>*0. To this end, we observe that, in the generic point (*P*0, *<sup>t</sup>*0), fixed values of *<sup>A</sup>βα*(Σ0), *<sup>C</sup>β*(Σ0), *<sup>B</sup>α*(Σ0), and *<sup>D</sup>*(Σ0) correspond to infinite vectors *<sup>y</sup>α*(*P*0, *<sup>t</sup>*0), because only *<sup>ω</sup>* components of *<sup>y</sup>α*(*P*0, *<sup>t</sup>*0) are determined by the balance equations while the remaining 9*ω* are completely arbitrary (see discussion in Section 3). Moreover, if all the *y<sup>α</sup>* in *H*<sup>0</sup> are over-ideal, the vector space *W*<sup>0</sup> would be empty, because the second law of thermodynamics prohibits that it contains over-ideal vectors. As a consequence, in (*P*0, *<sup>t</sup>*0), no process would be possible. On the other hand, since (*P*0, *<sup>t</sup>*0) is arbitrary, no thermodynamic transformation could occur in *<sup>B</sup>* in the interval of time [*τ*0, *<sup>τ</sup>*<sup>0</sup> + *<sup>τ</sup>*]. So, in (*P*0, *<sup>t</sup>*0), the space *<sup>H</sup>*<sup>0</sup> contains, in principle, both real/ideal vectors and over-ideal ones.

Let us suppose that in (*P*0, *<sup>t</sup>*0), the space *<sup>H</sup>*<sup>0</sup> contains an ideal vector *<sup>y</sup>*<sup>1</sup> *<sup>α</sup>* and an overideal vector *y*<sup>2</sup> *<sup>α</sup>*. Since the existence of *y*<sup>1</sup> *<sup>α</sup>* is possible if, and only if, (*P*0, *<sup>t</sup>*0) is in the thermodynamic equilibrium, while *y*<sup>2</sup> *<sup>α</sup>* exists if, and only if (*P*0, *<sup>t</sup>*0) is not in the thermodynamic equilibrium, such a situation is impossible to be realized.

Analogously, let us suppose that *y*<sup>1</sup> *<sup>α</sup>* is ideal and *y*<sup>2</sup> *<sup>α</sup>* is real. Again, such a situation is impossible, because it would require (*P*0, *<sup>t</sup>*0) to be in equilibrium and not in equilibrium.

Finally, let *y*<sup>1</sup> *<sup>α</sup>* be a real vector, and *y*<sup>2</sup> *<sup>α</sup>* be an over-ideal one. Such a situation is possible, in principle, provided (*P*0, *<sup>t</sup>*0) is not in equilibrium.

In such a case, due to the local formulation of the second law, neither *y*<sup>1</sup> *<sup>α</sup>* nor *y*<sup>2</sup> *<sup>α</sup>* are elements of *E*0.

Let us consider the linear combination *y*<sup>3</sup> *<sup>α</sup>* = *<sup>λ</sup>y*<sup>1</sup> *<sup>α</sup>* + (<sup>1</sup> <sup>−</sup> *<sup>λ</sup>*)*y*<sup>2</sup> *<sup>α</sup>*, with *<sup>λ</sup>* <sup>∈</sup>]0, 1[. Since *<sup>y</sup>*<sup>1</sup> *α* and *y*<sup>2</sup> *<sup>α</sup>* are in *H*0, they satisfy the following equations

$$A\_{\beta a}(\Sigma\_0) y\_a^1 = C\_{\beta}(\Sigma\_0), \tag{37}$$

$$A\_{\beta a}(\Sigma\_0) y\_a^2 = \mathbb{C}\_{\beta}(\Sigma\_0). \tag{38}$$

The combination of Equation (37) multiplied by *λ* and Equation (38) multiplied by (<sup>1</sup> <sup>−</sup> *<sup>λ</sup>*) leads to

$$A\_{\beta a}(\Sigma\_0) y\_a^3 = \mathbb{C}\_{\beta}(\Sigma\_0),\tag{39}$$

namely, *y*<sup>3</sup> *<sup>α</sup>* is also a solution of Equation (34), i.e., it is in *H*0. On the other hand, the local entropy production corresponding to *y*<sup>3</sup> *<sup>α</sup>* can be written as

$$\sigma^3 = \lambda \left[ B\_a(\Sigma\_0) y\_a^1 - D(\Sigma\_0) \right] + (1 - \lambda) \left[ B\_a(\Sigma\_0) y\_a^2 - D(\Sigma\_0) \right] = \tag{40}$$

$$= B\_a(\Sigma\_0) \left[ \lambda y\_a^1 + (1 - \lambda) y\_a^2 \right] - D(\Sigma\_0).$$

Since *λ* is arbitrary in ]0, 1[, nothing prevents choosing it as

$$
\lambda = \frac{D(\Sigma\_0) - B\_a(\Sigma\_0) y\_a^2}{B\_a(\Sigma\_0)[y\_a^1 - y\_a^2]},
\tag{41}
$$

because, as it is easily seen, the right-hand side of Equation (41) is in the interval ]0, 1[. In fact, being *y*<sup>2</sup> *<sup>α</sup>* over-ideal we have *<sup>D</sup>*(Σ0) <sup>−</sup> *<sup>B</sup>α*(Σ0)*y*<sup>2</sup> *<sup>α</sup>* > 0. Moreover, being *y*<sup>1</sup> *<sup>α</sup>* real, we have *<sup>B</sup>α*(Σ0) *y*1 *<sup>α</sup>* <sup>−</sup> *<sup>y</sup>*<sup>2</sup> *α* <sup>&</sup>gt; *<sup>D</sup>*(Σ0) <sup>−</sup> *<sup>D</sup>*(Σ0), namely, *<sup>B</sup>α*(Σ0) *y*1 *<sup>α</sup>* <sup>−</sup> *<sup>y</sup>*<sup>2</sup> *α* > 0. Hence, *λ* > 0. Moreover, being *y*<sup>1</sup> *<sup>α</sup>* real, we also have *<sup>B</sup>α*(Σ0) *y*1 *<sup>α</sup>* <sup>−</sup> *<sup>y</sup>*<sup>2</sup> *α* <sup>&</sup>gt; *<sup>D</sup>*(Σ0) <sup>−</sup> *<sup>B</sup>α*(Σ0)*y*<sup>2</sup> *<sup>α</sup>* and, hence, *λ* < 1.

Consequently, the right-hand side of Equation (40) vanishes, so that *y*<sup>3</sup> *<sup>α</sup>* is in *E*0. However, this is impossible, otherwise (*P*0, *<sup>t</sup>*0) would be in a thermodynamic equilibrium. Thus, it is forbidden that in (*P*0, *<sup>t</sup>*0) there are both real and over-ideal vectors that are solutions of the local balance laws (34).

Furthermore, suppose that both *y*<sup>1</sup> *<sup>α</sup>* and *y*<sup>2</sup> *<sup>α</sup>* are real. Then, it is easy to verify by direct calculations that *λ* can be taken, such that *σ*<sup>3</sup> > 0.

Finally, if (*P*0, *<sup>t</sup>*0) is a point of equilibrium, then the entropy production related to *<sup>y</sup>*<sup>1</sup> *α* and *y*<sup>2</sup> *<sup>α</sup>* vanishes, so that by Equation (40), it follows that *σ*<sup>3</sup> is zero.

The considerations above show the impossibility that, in a point *P*<sup>0</sup> of *B*, at a given instant *t*0, the solutions of Equation (34) can be of different types. Moreover, they cannot be over-ideal only, because this contradicts the local form of the second law of thermodynamics. Thus, *<sup>H</sup>*<sup>0</sup> may contain either only real vectors, and in such a case, (*P*0, *<sup>t</sup>*0) is a point of non-equilibrium, or only ideal vectors, and in such a case, (*P*0, *<sup>t</sup>*0) is a point of equilibrium. This conclusion proves the theorem.

### **Corollary 2.** H *=* W*.*

**Proof.** This corollary is an immediate consequence of the arbitrariness of the initial instant *t*0, and of the point *P*0. In particular, whatever *t*<sup>0</sup> is, we can ever consider it as the initial instant of the restricted process of duration *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>* <sup>−</sup> *<sup>t</sup>*0, so that *<sup>H</sup>*<sup>0</sup> <sup>≡</sup> *Ht*(*P*0, *<sup>t</sup>*0) has dimension 10*ω*. Moreover, only 7*ω* of components of the vectors of *H*<sup>0</sup> can be determined by the algebraic relations (34) and (36) while the further 3*ω* components are completely arbitrary. Thus, to *H*<sup>0</sup> can be applied to the conclusions established in Theorem 1. This is enough to prove that, for any *<sup>t</sup>* <sup>∈</sup> [*τ*0, *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>*], the space of the higher derivatives *Ht* contains only real or ideal vectors.

**Remark 7.** *The Corollary <sup>2</sup> also implies* <sup>E</sup><sup>ˆ</sup> *<sup>=</sup>* <sup>E</sup>*.*

**Corollary 3.** *The unilateral differential constraint* (32) *is a restriction on the constitutive quantities Uβ, rβ, s, and Jk, and not on the thermodynamic processes p.*

**Proof.** In fact, any process *<sup>p</sup>* : *<sup>t</sup>* <sup>∈</sup> [*τ*0, *<sup>τ</sup>*<sup>0</sup> <sup>+</sup> *<sup>τ</sup>*] <sup>⊆</sup> <sup>R</sup> <sup>→</sup> *<sup>z</sup>α*(*xj*, *<sup>t</sup>*) ∈ C, where *<sup>z</sup>α*(*xj*, *<sup>t</sup>*) is a solution of the balance laws (29), can only be either irreversible or reversible but not over-reversible, because otherwise its representative curve *γ* would contain at least an over-ideal point against Corollary 2. On the other hand, such a property of the solutions of the system of balance laws is not guaranteed, whatever *Aβα* and *C<sup>β</sup>* are, and for arbitrary *s* and *Jk*, because, given the state space, only particular forms of those functions defined on it lead to a nonnegative entropy production. Then, the role of the unilateral differential constraint in Equation (32) is just to select such forms.

#### **5. Discussion**

Exploitation of the second law of thermodynamics is based on the assumption that the constitutive equations modeling material properties must be assigned in such a way that all solutions of the field equations satisfy the entropy inequality. An alternative interpretation of the second law consequences is that we must exclude from the set of solutions of the balance equations ones that do not guarantee a nonnegative entropy production. The problem of the choice between the two interpretations above was solved in 1996 by Muschik and Ehrentraut [13], by postulating an amendment to the second law, which states that, at a fixed instant of time and in any point of the body, the entropy production is zero if, and only if, this point is in a thermodynamic equilibrium. Muschik and Ehrentraut proved that, under the validity of the amendment, the constitutive equations and not the thermodynamic processes are restricted by the second law of thermodynamics. Such a result provides the theoretical basis to the rigorous mathematical procedures for the exploitation of the second law proposed by Coleman and Noll in 1963 [1], and by Liu in 1972 [12].

In the present paper, we revisited the Muschik and Ehrentraut theorem, highlighting some geometric aspects that were hidden in reference [13]. Moreover, we proposed a generalized formulation of the second law of thermodynamics that incorporates the amendment. Progresses in the analysis of the entropy principle achieved in this way are the following:

1. In the analysis of the efficiency of the thermodynamic systems, the concept of the Carnot cycle plays a fundamental role. Carnot was the first to prove that the efficiency of a quasi-static transformation is maximum for a Carnot cycle. However, since quasi-static transformations require an infinite time, the Carnot efficiency is not suitable to describe the performance of real heat engines, which produce irreversible transformations taking over in a finite time. Thus, for real processes, it is important to investigate how much the efficiency deteriorates when the cycle is operated in a finite time. This is the task of finite time thermodynamics, a modern non-equilibrium theory, which has been developed in the last four decades [16–18]. Since the global formulation of the second law of thermodynamics proposed here explicitly mentions the fact that a reversible transformation lies in the equilibrium bundle E, it seems to be more appropriate for finite time thermodynamics with respect to the classical formulation, which does not contain such information.


**Figure 1.** A schematized explanation of the analysis of the entropy principle developed in the present paper, and of its relation to the Muschik and Ehrentraut theory. An Amendment to second law of thermodynamics has been proposed by Muschik and Ehrentraut in 1996 [13]; rigorous procedures for the exploitation of the entropy principle have been formulated by Coleman and Noll in 1963 [1], and by Liu in 1972 [12]; the general tenets of Finite Time Thermodynamics have been summarized by Andresen, Salomon and Berry in 1984 [17].

In reference [21], the results in [13] were extended in order to encompass non-regular processes. Such an investigation deserves consideration, since thermodynamic processes involving discontinuous solutions are frequent in physics. In future research studies, we will aim to reanalyze the results in [21], within the mathematical framework presented here.

**Author Contributions:** Conceptualization, V.A.C. and P.R.; formal analysis, V.A.C. and P.R.; investigation, V.A.C. and P.R.; writing—original draft, V.A.C. and P.R.; writing—review and editing, V.A.C. and P.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the University of Basilicata (RIL 2020), the University of Messina (FFABR Unime 2020), and the Italian National Group of Mathematical Physics (GNFM-INdAM).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Patrizia Rogolino thanks the University of Messina and the Italian National Group of Mathematical Physics (GNFM-INdAM) for financial support. Vito Antonio Cimmelli thanks the University of Basilicata and the Italian National Group of Mathematical Physics (GNFM-INdAM) for financial support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Integrability of the Multi-Species TASEP with Species-Dependent Rates**

**Eunghyun Lee**

Department of Mathematics, Nazarbayev University, Nur-Sultan 010000, Kazakhstan; eunghyun.lee@nu.edu.kz

**Abstract:** Assume that each species *l* has its own jump rate *bl* in the multi-species totally asymmetric simple exclusion process. We show that this model is *integrable* in the sense that the Bethe ansatz method is applicable to obtain the transition probabilities for all possible *N*-particle systems with up to *N* different species.

**Keywords:** TASEP; ASEP; multi-species; Bethe ansatz

#### **1. Introduction**

The multi-species asymmetric simple exclusion process on Z is a generalization of the asymmetric simple exclusion process (ASEP) on Z in the sense that each particle may belong to a different species labelled by an integer *l* ∈ {1, 2, ···}. Each particle jumps to the right by one step with the probability *p* or to the left by one step with the probability *<sup>q</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>p</sup>* after a waiting time that is exponentially distributed with rate 1. If a particle belonging to *l* attempts to jump to the site occupied by a particle belonging to *l* ≥ *l*, the jump is prohibited; however, if a particle belonging to *l* attempts to jump to the site occupied by a particle belonging to *l* < *l* , then the jump occurs by interchanging positions.

The transition probabilities and some determinantal formulas for the multi-species ASEP and its special cases were found in [1–5]. For certain special initial conditions with a single second class particle, some distributions and their asymptotics were studied in [6,7]. More recently, asymptotic behaviors of the second class particles were studied using the color-position symmetry, see [8]. In fact, the multi-species asymmetric simple exclusion process can be considered in a more general context—that is, the coloured six vertex model [9].

Another direction of generalizing the ASEP and other models studied in the integrable probability is to make the jump rates inhomogeneous. It is known that the Bethe ansatz method is still applicable to some single-species model with inhomogeneous jump rates. The basic idea of using the Bethe ansatz in the ASEP is that the generator of the ASEP is a similarity transformation of the XXZ quantum spin system. Considering that the Bethe ansatz is a method to find eigenvalues and eigenvectors of a certain class of quantum spin systems, we use the Bethe ansatz to find the solution of the forward equation of a certain class of Markov processes, that is, the transition probabilities of the processes.

Of course, for some particle models, the Bethe ansatz method cannot be used. For the background of Bethe ansatz, see [10,11]. It is known that the Bethe ansatz is applicable to some generalization of the ASEP. For example, the transition probability and the current distribution of the totally asymmetric simple exclusion process (TASEP) with particledependent rates were studied in [12], and the transition probabilities and some asymptotic results for the *q*-deformed totally asymmetric zero range process with site-dependent rates were studied in [13,14].

In this paper, we consider the multi-species totally asymmetric simple exclusion processes with *N* particles in which particles move to the right and each species *l* is allowed to have its own rate *bl*. Following the notations used in [4], let *<sup>X</sup>* = (*x*1, ··· , *xN*) <sup>∈</sup> <sup>Z</sup>*<sup>N</sup>* with *<sup>x</sup>*<sup>1</sup> <sup>&</sup>lt; ··· <sup>&</sup>lt; *xN* represent the positions of particles, and let *<sup>π</sup>* <sup>=</sup> *<sup>π</sup>*(1)*π*(2) ··· *<sup>π</sup>*(*N*)

**Citation:** Lee, E. Integrability of the Multi-Species TASEP with Species-Dependent Rates. *Symmetry* **2021**, *13*, 1578. https://doi.org/10.3390/ sym13091578

Academic Editors: Róbert Kovács, Patrizia Rogolino, Francesco Oliveri and Jaume Giné

Received: 28 July 2021 Accepted: 21 August 2021 Published: 27 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

be a permutation of a multi-set <sup>M</sup> = [*i*1, ··· , *iN*] with elements taken from {1, ··· , *<sup>N</sup>*} and *π*(*i*) representing the species of the *i th* leftmost particle. Then, the state of an *N*-particle system is denoted by

$$(X,\pi) = (\mathfrak{x}\_1, \dots, \mathfrak{x}\_{N}, \pi(\mathbf{1}), \cdot, \dots, \pi(N)).$$

Let us write *<sup>P</sup>*(*Y*,*ν*)(*X*, *<sup>π</sup>*; *<sup>t</sup>*) for the transition probability from (*Y*, *<sup>ν</sup>*) at *<sup>t</sup>* <sup>=</sup> 0 to (*X*, *<sup>π</sup>*) at a later time *<sup>t</sup>*. For fixed *<sup>X</sup>* and *<sup>Y</sup>*, *<sup>P</sup>*(*Y*,*ν*)(*X*, *<sup>π</sup>*; *<sup>t</sup>*) is regarded as a matrix element of an *<sup>N</sup><sup>N</sup>* <sup>×</sup> *<sup>N</sup><sup>N</sup>* matrix denoted by **<sup>P</sup>***Y*(*X*; *<sup>t</sup>*) whose columns and rows are labelled by *<sup>ν</sup>*, *<sup>π</sup>* <sup>=</sup> <sup>1</sup> ··· 11, 1 ··· 12, ··· , *<sup>N</sup>* ··· *<sup>N</sup>*, respectively. Throughout this paper, given an *<sup>N</sup><sup>n</sup>* <sup>×</sup> *<sup>N</sup><sup>n</sup>* matrix, we assume that its rows (*i*<sup>1</sup> ··· *in*) and columns (*j*<sup>1</sup> ··· *jn*) are labelled by 1 ··· 1, ··· , *n* ··· *n* and that these labels are listed lexicographically, unless stated otherwise. The main result of this paper is that the multi-species TASEP with species-dependent rates is an integrable model, and we provide a formula analogous to (2.12) in [4] using the Bethe ansatz method.

#### *Statement of the Results*

We first introduce a few objects to state the main theorem. Define an *<sup>N</sup>*<sup>2</sup> <sup>×</sup> *<sup>N</sup>*<sup>2</sup> matrix **R***βα* = *Rij*,*kl* with

$$R\_{ij,kl} = \begin{cases} S\_{\beta a}(i) & \text{if } ij = kl \text{ with } i \le j; \\ -1 & \text{if } ij = kl \text{ with } i > j; \\ T\_{\beta a}(i) & \text{if } ij = lk \text{ with } i < j; \\ 0 & \text{for all other cases,} \end{cases} \tag{1}$$

where

$$S\_{\not\beta a}(i) = -\frac{1 - b\_i \mathbb{I}\_{\beta}}{1 - b\_i \mathbb{I}\_{\beta}}, \quad T\_{\not\alpha}(i) = \frac{b\_i(\mathbb{I}\_{\beta} - \mathbb{I}\_{\alpha})}{1 - b\_i \mathbb{I}\_{\alpha}}, \quad \mathbb{I}\_{\alpha} \mathbb{I}\_{\beta} \in \mathbb{C}. \tag{2}$$

**Remark 1.** *The form of the matrix (1) was obtained by induction via similar arguments to Sections 2.1 and 2.2 in [3], which treats a special case, and the motivation of (2) is given in Section 2.1. Finding the form of (1) with (2) is the key idea of this paper.*

Let T*<sup>l</sup>* be the simple transposition that interchanges the number at the *l*th slot and the number at the (*<sup>l</sup>* <sup>+</sup> <sup>1</sup>)st slot. If <sup>T</sup>*<sup>l</sup>* maps a permutation ( ··· *αβ* ···) to ( ··· *βα* ···), we write <sup>T</sup>*<sup>l</sup>* <sup>=</sup> <sup>T</sup>*l*(*β*, *<sup>α</sup>*) when necessary. Corresponding to a simple transposition <sup>T</sup>*l*(*β*, *<sup>α</sup>*), we define *<sup>N</sup><sup>N</sup>* <sup>×</sup> *<sup>N</sup><sup>N</sup>* matrix **<sup>T</sup>***l*(*β*, *<sup>α</sup>*) by the tensor product of matrices,

$$\mathbf{T}\_{l}(\boldsymbol{\beta},a) = \underbrace{\mathbf{I}\_{N} \otimes \dots \otimes \mathbf{I}\_{N}}\_{(l-1)\text{ times}} \otimes \mathbf{R}\_{\text{fa}} \otimes \underbrace{\mathbf{I}\_{N} \otimes \dots \otimes \mathbf{I}\_{N}}\_{(N-l-1)\text{ times}}\tag{3}$$

where **I***<sup>N</sup>* is the *N* × *N* identity matrix. For a permutation *σ* in the symmetric group S*<sup>N</sup>* written

$$
\sigma = \mathcal{T}\_{\mathbf{i}\_{\mathbf{j}}} \cdots \mathcal{T}\_{\mathbf{i}\_{\mathbf{i}}} = \mathcal{T}\_{\mathbf{i}\_{\mathbf{j}}}(\boldsymbol{\beta}\_{\mathbf{j}}, \boldsymbol{\alpha}\_{\mathbf{j}}) \cdots \mathcal{T}\_{\mathbf{i}\_{\mathbf{i}}}(\boldsymbol{\beta}\_{\mathbf{i}}, \boldsymbol{\alpha}\_{\mathbf{i}}) \tag{4}
$$

for some *i*1, ··· , *ij* ∈ {1, ··· , *N* − 1}, we define

$$\mathbf{A}\_{\sigma} = \mathbf{T}\_{i\_j}(\beta\_{j'}\boldsymbol{\alpha}\_j) \cdot \cdots \cdot \mathbf{T}\_{i\_1}(\beta\_1, \boldsymbol{\alpha}\_1). \tag{5}$$

Here, **A***<sup>σ</sup>* is well defined, that is, **A***<sup>σ</sup>* is unique regardless of the representation of *σ* by simple transpositions. This *well-definedness* is due to the following lemma.

**Lemma 1.** *The following consistency relations are satisfied.*


The relations in Lemma <sup>1</sup> with *bl* <sup>=</sup> 1 for all *<sup>l</sup>* are already known for the multi-species ASEP.

**Remark 2.** *The definitions of* <sup>T</sup>*<sup>l</sup> and* **<sup>A</sup>***<sup>σ</sup> are motivated by the arguments for <sup>N</sup>* <sup>=</sup> 2, 3 *in Sections 2.1 and 2.2 in [3], which treats a special case.*

Let **<sup>J</sup>**(*t*) be the *<sup>N</sup><sup>N</sup>* <sup>×</sup> *<sup>N</sup><sup>N</sup>* diagonal matrix whose (*π*, *<sup>π</sup>*)-element is given by *<sup>e</sup>εππ<sup>t</sup>* where

$$\varepsilon\_{\pi\pi\pi} = \varepsilon\_{\pi\pi\pi} (\xi\_1 \wedge \dots \wedge \xi\_N) = \frac{1}{\tilde{\xi}\_1} + \dots + \frac{1}{\tilde{\xi}\_N} - \left(b\_{\pi(1)} + \dots + b\_{\pi(N)}\right), \quad \xi\_1 \wedge \dots \wedge \xi\_N \in \mathbb{C}\_{\succ}$$

let **<sup>D</sup>**(*x*1, ··· , *xN*) be the *<sup>N</sup><sup>N</sup>* <sup>×</sup> *<sup>N</sup><sup>N</sup>* diagonal matrix whose (*π*, *<sup>π</sup>*)-element is given by *b x*1 *<sup>π</sup>*(1) ··· *b xN π*(*N*) , and let **D** (*y*1, ··· , *yN*) be the *<sup>N</sup><sup>N</sup>* <sup>×</sup> *<sup>N</sup><sup>N</sup>* diagonal matrix whose (*ν*, *<sup>ν</sup>*) element is given by *b* −*y*<sup>1</sup> *<sup>ν</sup>*(1) ··· *b* −*yN <sup>ν</sup>*(*N*) where *yi*s are the initial positions. In the next theorem, the integral of a matrix implies that the integral is taken element-wise, and − # implies <sup>1</sup> 2*πi* # .

**Theorem 1.** *Let* **A***σ be given as in (5) and c be a positively oriented circle centred at the origin with a radius less than bl for all <sup>l</sup> in the complex plane* C*. Then, the matrix of the transition probabilities of the multi-species TASEP with species-dependent rates is*

$$\mathbf{P}\_{Y}(X;t) = \sum\_{\sigma \in \mathfrak{S}\_{N}} \int\_{\mathcal{L}} \cdots \int\_{\mathcal{L}} \mathbf{J}(t) \mathbf{D}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N}) \mathbf{A}\_{\sigma} \mathbf{D}'(y\_{1}, \cdots, y\_{N}) \prod\_{l=1}^{N} \left(\tilde{\xi}\_{\sigma(l)}^{\mathbf{x}\_{l} - y\_{\sigma(l)} - 1}\right) d\tilde{\xi}\_{1} \cdots d\tilde{\xi}\_{N}. \tag{6}$$

**Remark 3.** *If <sup>ν</sup>* <sup>=</sup> <sup>12</sup> ··· *N, in other words, all <sup>N</sup> particles belong to different species, and the species are initially arranged in ascending order, then the system is the same as the TASEP with particle-dependent rates studied in [12]. Hence, the transition probability <sup>P</sup>*(*Y*, 12 ··· *<sup>N</sup>*)(*X*, 12 ··· *<sup>N</sup>*; *<sup>t</sup>*) *can be expressed as a determinant (see Theorem 1 in [12]).*

**Remark 4.** *Theorem <sup>1</sup> partially extends (2.12) in [4]. In other words, (6) with bl* = <sup>1</sup> *is equal to (2.12) in [4] with p* = 1*.*

The proofs of Lemma 1 and Theorem 1 are given in the next section.

#### **2. Proof of Theorem 1**

In order to prove that the (*π*, *<sup>ν</sup>*)-element of the right-hand side of (6) is *<sup>P</sup>*(*Y*,*ν*)(*X*, *<sup>π</sup>*; *<sup>t</sup>*), we should show that the (*π*, *ν*)-element satisfies its forward equation and the initial condition *<sup>P</sup>*(*Y*,*ν*)(*X*, *<sup>π</sup>*; 0) = *<sup>δ</sup>X*,*Y*.

#### *2.1. Forward Equations*

We first study the two-particle systems, which will be building-blocks for the formulas for *<sup>N</sup>*-particle systems. When *<sup>x</sup>*<sup>1</sup> <sup>&</sup>lt; *<sup>x</sup>*<sup>2</sup> <sup>−</sup> 1, the forward equations of *<sup>P</sup>*(*Y*,*ν*)(*x*1, *<sup>x</sup>*2, *<sup>π</sup>*; *<sup>t</sup>*) are straightforward because two particles act as *free* particles. Hence, the forward equations of *<sup>P</sup>*(*Y*,*ν*)(*x*1, *<sup>x</sup>*2, *<sup>π</sup>*; *<sup>t</sup>*) are expressed as

$$\begin{aligned} \frac{d}{dt} \mathbf{P}\_{Y}(\mathbf{x}\_{1}, \mathbf{x}\_{2}; t) &= \underbrace{\begin{bmatrix} b\_{1} & 0 & 0 & 0 \\ 0 & b\_{1} & 0 & 0 \\ 0 & 0 & b\_{2} & 0 \\ 0 & 0 & 0 & b\_{2} \end{bmatrix}}\_{(r\_{1})} \mathbf{P}\_{Y}(\mathbf{x}\_{1} - 1, \mathbf{x}\_{2}; t) + \underbrace{\begin{bmatrix} b\_{1} & 0 & 0 & 0 \\ 0 & b\_{2} & 0 & 0 \\ 0 & 0 & b\_{1} & 0 \\ 0 & 0 & 0 & b\_{2} \end{bmatrix}}\_{(r\_{2})} \mathbf{P}\_{Y}(\mathbf{x}\_{1}, \mathbf{x}\_{2} - 1; t) \\ &- \begin{bmatrix} b\_{1} + b\_{1} & 0 & 0 & 0 \\ 0 & b\_{1} + b\_{2} & 0 & 0 \\ 0 & 0 & b\_{2} + b\_{1} & 0 \\ 0 & 0 & 0 & b\_{2} + b\_{2} \end{bmatrix} \mathbf{P}\_{Y}(\mathbf{x}\_{1}, \mathbf{x}\_{2}; t) \end{aligned} \tag{7}$$

where the derivative of the matrix **<sup>P</sup>***Y*(*x*1, *<sup>x</sup>*2; *<sup>t</sup>*) on the left-hand side implies the matrix of the derivatives of elements of **<sup>P</sup>***Y*(*x*1, *<sup>x</sup>*2; *<sup>t</sup>*). The matrices (*r*1) and (*r*2) account for the *probability current* in the states (*x*1, *<sup>x</sup>*2, *<sup>π</sup>*) by a particle's jump to the right next site, which is empty. On the other hand, when *<sup>x</sup>*<sup>1</sup> <sup>=</sup> *<sup>x</sup>*<sup>2</sup> <sup>−</sup> 1, if two particles belong to different species, two particles may swap their positions. For example, if the initial state is (*Y*, 12), the system cannot be at (*X*, 21) at any later time *<sup>t</sup>*. Hence, the forward equation of *<sup>P</sup>*(*Y*,12)(*x*1, *<sup>x</sup>*<sup>1</sup> <sup>+</sup> 1, 12; *t*) is

$$\frac{d}{dt}P\_{(Y,12)}(\mathbf{x}\_1,\mathbf{x}\_1+1,12;t) = b\_1 P\_{(Y,12)}(\mathbf{x}\_1-1,\mathbf{x}\_1+1,12;t) - b\_2 P\_{(Y,12)}(\mathbf{x}\_1,\mathbf{x}\_1+1,12;t).$$

On the other hand, *<sup>P</sup>*(*Y*,12)(*x*1, *<sup>x</sup>*<sup>1</sup> <sup>+</sup> 1, 21; *<sup>t</sup>*) = 0 for all *<sup>t</sup>*, because the model is totally asymmetric. If the initial state is (*Y*, 21), the forward equation of *<sup>P</sup>*(*Y*,21)(*x*1, *<sup>x</sup>*<sup>1</sup> <sup>+</sup> 1, 21; *<sup>t</sup>*) is

$$\frac{d}{dt}P\_{\{\mathbf{Y},21\}}(\mathbf{x}\_1,\mathbf{x}\_1+1,21;t) = b\_2 P\_{\{\mathbf{Y},21\}}(\mathbf{x}\_1-1,\mathbf{x}\_1+1,21;t) - (b\_2+b\_1) P\_{\{\mathbf{Y},21\}}(\mathbf{x}\_1,\mathbf{x}\_1+1,21;t)$$

and the forward equation of *<sup>P</sup>*(*Y*,21)(*x*1, *<sup>x</sup>*<sup>1</sup> <sup>+</sup> 1, 12; *<sup>t</sup>*) is

$$\begin{aligned} \frac{d}{dt}P\_{(Y,21)}(\mathbf{x}\_1, \mathbf{x}\_1 + 1, 12; t) &= b\_1 P\_{(Y,21)}(\mathbf{x}\_1 - 1, \mathbf{x}\_1 + 1, 12; t) + b\_2 P\_{(Y,21)}(\mathbf{x}\_1, \mathbf{x}\_1 + 1, 21; t) \\ &- b\_2 P\_{(Y,21)}(\mathbf{x}\_1, \mathbf{x}\_1 + 1, 12; t). \end{aligned}$$

Hence, the forward equations of *<sup>P</sup>*(*Y*,*ν*)(*x*1, *<sup>x</sup>*<sup>1</sup> <sup>+</sup> 1, *<sup>π</sup>*; *<sup>t</sup>*) are expressed as

$$\begin{aligned} \frac{d}{dt} \mathbf{P}\_{Y}(\mathbf{x}\_{1}, \mathbf{x}\_{1} + 1; t) &= \\\\ \begin{bmatrix} b\_{1} & 0 & 0 & 0 \\ 0 & b\_{1} & 0 & 0 \\ 0 & 0 & b\_{2} & 0 \\ 0 & 0 & 0 & b\_{2} \end{bmatrix} \mathbf{P}\_{Y}(\mathbf{x}\_{1} - 1, \mathbf{x}\_{1} + 1; t) + \underbrace{\begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & b\_{2} & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}}\_{(B)} \mathbf{P}\_{Y}(\mathbf{x}\_{1}, \mathbf{x}\_{1} + 1; t) \\\\ -\underbrace{\left( \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & b\_{2} & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}}\_{(C)} + \begin{bmatrix} b\_{1} & 0 & 0 & 0 \\ 0 & b\_{2} & 0 & 0 \\ 0 & 0 & b\_{1} & 0 \\ 0 & 0 & 0 & b\_{2} \end{bmatrix} \right) \mathbf{P}\_{Y}(\mathbf{x}\_{1}, \mathbf{x}\_{1} + 1; t). \end{aligned} \tag{8}$$

Here, the matrix (*B*) accounts for *probability current* going in the states (*x*1, *<sup>x</sup>*<sup>1</sup> + 1, 12) by the species-2 particle's jump from the state (*x*1, *<sup>x</sup>*<sup>1</sup> + 1, 21). Similarly, the matrix (*C*) accounts for *probability current* going out of the states (*x*1, *<sup>x</sup>*<sup>1</sup> + 1, 21) by species-2 particle's jump to the state (*x*1, *<sup>x</sup>*<sup>1</sup> + 1, 12). Equations (7) and (8) imply that, if **<sup>U</sup>**(*x*1, *<sup>x</sup>*2; *<sup>t</sup>*) is a <sup>4</sup> <sup>×</sup> 4 matrix whose elements are functions on <sup>Z</sup><sup>2</sup> <sup>×</sup> [0, <sup>∞</sup>), then the forward equation of *<sup>P</sup>*(*Y*,*ν*)(*x*1, *<sup>x</sup>*2, *<sup>π</sup>*; *<sup>t</sup>*) for any *<sup>x</sup>*<sup>1</sup> <sup>&</sup>lt; *<sup>x</sup>*<sup>2</sup> is in the form of the (*π*, *<sup>ν</sup>*)-element of

$$\frac{d}{dt}\mathbf{U}(\mathbf{x}\_{1},\mathbf{x}\_{2};t) = \begin{bmatrix} b\_{1} & 0 & 0 & 0 \\ 0 & b\_{1} & 0 & 0 \\ 0 & 0 & b\_{2} & 0 \\ 0 & 0 & 0 & b\_{2} \end{bmatrix} \mathbf{U}(\mathbf{x}\_{1}-1,\mathbf{x}\_{2};t) + \begin{bmatrix} b\_{1} & 0 & 0 & 0 \\ 0 & b\_{2} & 0 & 0 \\ 0 & 0 & b\_{1} & 0 \\ 0 & 0 & 0 & b\_{2} \end{bmatrix} \mathbf{U}(\mathbf{x}\_{1},\mathbf{x}\_{2}-1;t)$$
 
$$-\begin{bmatrix} b\_{1}+b\_{1} & 0 & 0 & 0 \\ 0 & b\_{1}+b\_{2} & 0 & 0 \\ 0 & 0 & b\_{2}+b\_{1} & 0 \\ 0 & 0 & 0 & b\_{2}+b\_{2} \end{bmatrix} \mathbf{U}(\mathbf{x}\_{1},\mathbf{x}\_{2};t)$$

subject to the (*π*, *ν*)-element of

$$
\begin{bmatrix} b\_1 & 0 & 0 & 0 \\ 0 & b\_2 & 0 & 0 \\ 0 & 0 & b\_1 & 0 \\ 0 & 0 & 0 & b\_2 \end{bmatrix} \mathbf{U}(\mathbf{x}\_1, \mathbf{x}\_1; t) = \begin{bmatrix} b\_1 & 0 & 0 & 0 \\ 0 & b\_1 & b\_2 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & b\_2 \end{bmatrix} \mathbf{U}(\mathbf{x}\_1, \mathbf{x}\_1 + 1; t).
$$

Now, we extend the argument for two-particle systems to *N*-particle systems. The matrices (*r*1) and (*r*2) in (7) for two-particle systems are generalized to

$$\mathbf{r}\_{i} = \underbrace{\mathbf{I}\_{N} \otimes \cdots \otimes \mathbf{I}\_{N}}\_{i-1} \otimes \mathbf{r}\_{i} \otimes \underbrace{\mathbf{I}\_{N} \otimes \cdots \otimes \mathbf{I}\_{N}}\_{N-i} \quad i = 1, \cdots, N$$

where **r** is the diagonal matrix,

$$\mathbf{r} = \begin{bmatrix} b\_1 \\ & \ddots \\ & & b\_N \end{bmatrix}.$$

The matrix (*B*) in (8) is generalized to an *<sup>N</sup>*<sup>2</sup> <sup>×</sup> *<sup>N</sup>*<sup>2</sup> matrix **<sup>B</sup>** <sup>=</sup> *Bij*,*kl* with

> *Bij*,*kl* = ) *bj* if *ij* = *lk* with *<sup>i</sup>* < *<sup>j</sup>*; 0 for all other cases,

and let

$$\mathbf{B}\_{i} = \underbrace{\mathbf{I}\_{N} \otimes \dots \otimes \mathbf{I}\_{N}}\_{i-1} \otimes \mathbf{B} \otimes \underbrace{\mathbf{I}\_{N} \otimes \dots \otimes \mathbf{I}\_{N}}\_{N-i-1}.$$

The matrix (*C*) in (8) is generalized to *<sup>N</sup>*<sup>2</sup> <sup>×</sup> *<sup>N</sup>*<sup>2</sup> matrix **<sup>C</sup>** <sup>=</sup> *Cij*,*kl* with

$$\mathbb{C}\_{ij,kl} = \begin{cases} b\_i & \text{if } ij = kl \text{ with } i < j; \\ 0 & \text{for all other cases,} \end{cases}$$

and let

$$\mathbf{C}\_{i} = \underbrace{\mathbf{I}\_{N} \otimes \dots \otimes \mathbf{I}\_{N}}\_{i-1} \otimes \mathbf{C} \otimes \underbrace{\mathbf{I}\_{N} \otimes \dots \otimes \mathbf{I}\_{N}}\_{N-i-1}.$$

All forward equations of *<sup>P</sup>*(*Y*,*ν*)(*x*1, ··· , *xN*, *<sup>π</sup>*; *<sup>t</sup>*) may be expressed as a matrix equation. For example, if *xi* <sup>&</sup>lt; *xi*+<sup>1</sup> <sup>−</sup> 1 for all *<sup>i</sup>*, then the forward equation of *<sup>P</sup>*(*Y*,*ν*)(*x*1, ··· , *xN*, *<sup>π</sup>*; *<sup>t</sup>*) is the (*π*, *ν*)-element of

$$\begin{split} \frac{d}{dt} \mathbf{P}\_{\mathbf{Y}} (\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N}; t) &= \mathbf{r}\_{1} \mathbf{P}\_{\mathbf{Y}} (\mathbf{x}\_{1} - 1, \mathbf{x}\_{2}, \cdots, \mathbf{x}\_{N}; t) + \cdots \\ &\quad + \mathbf{r}\_{N} \mathbf{P}\_{\mathbf{Y}} (\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N-1}, \mathbf{x}\_{N} - 1; t) - (\mathbf{r}\_{1} + \cdots + \mathbf{r}\_{N}) \mathbf{P}\_{\mathbf{Y}} (\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N}; t) \\ \text{and if } \mathbf{x} = \mathbf{x}\_{i} &= \mathbf{x}\_{i+1} - 1 \text{ and } \mathbf{x}\_{j} < \mathbf{x}\_{j+1} - 1 \text{ for all } j \neq i, \end{split} \tag{9}$$

$$\begin{split} \frac{d}{dt} \mathbf{P}\_{Y}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{i-1}, \mathbf{x}, \mathbf{x} + 1, \mathbf{x}\_{i}, \mathbf{x}\_{i+2}, \cdots, \mathbf{x}\_{Ni}; t) &= \\ \frac{1}{j+i+1} \mathbf{P}\_{Y} \mathbf{P}\_{Y}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{j-1}, \mathbf{x}\_{j} - 1, \mathbf{x}\_{j+1}, \cdots, \mathbf{x}\_{N}; t) + \mathbf{B}\_{i} \mathbf{P}\_{Y} (\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{i-1}, \mathbf{x}, \mathbf{x} + 1, \mathbf{x}\_{i+2}, \cdots, \mathbf{x}\_{N}; t) \\ - \sum\_{j \neq i} \mathbf{r}\_{j} \mathbf{P}\_{Y} (\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{j-1}, \mathbf{x}, \mathbf{x} + 1, \mathbf{x}\_{j+2}, \cdots, \mathbf{x}\_{N}; t) - \mathbf{C}\_{i} \mathbf{P}\_{Y} (\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{i-1}, \mathbf{x}, \mathbf{x} + 1, \mathbf{x}\_{i}, \mathbf{x}\_{i+2}, \cdots, \mathbf{x}\_{N}; t). \end{split} \tag{10}$$

For other configurations of (*x*1, ··· , *xN*), the form of the matrix of the forward equations may be different from (9) and (10). However, as in other Bethe ansatz applicable models, if **<sup>U</sup>**(*x*1, ··· , *xN*; *<sup>t</sup>*) = *<sup>U</sup>πν*(*x*1, ··· , *xN*; *<sup>t</sup>*) is an *<sup>N</sup><sup>N</sup>* <sup>×</sup> *<sup>N</sup><sup>N</sup>* matrix whose elements *<sup>U</sup>πν*(*x*1, ··· , *xN*; *<sup>t</sup>*) are functions on <sup>Z</sup>*<sup>N</sup>* <sup>×</sup> [0, <sup>∞</sup>), then the forward equation of *<sup>P</sup>*(*Y*,*ν*)(*x*1, ··· , *xN*, *<sup>π</sup>*; *<sup>t</sup>*) for any *<sup>x</sup>*<sup>1</sup> <sup>&</sup>lt; ··· <sup>&</sup>lt; *xN* is in the form of the (*π*, *<sup>ν</sup>*)-element of

$$\frac{d}{dt}\mathbf{U}(\mathbf{x}\_{1},\cdots,\mathbf{x}\_{N};t) = \sum\_{j=1}^{N} \mathbf{r}\_{j} \mathbf{U}(\mathbf{x}\_{1},\cdots,\mathbf{x}\_{j-1},\mathbf{x}\_{j} - 1,\mathbf{x}\_{j+1},\cdots,\mathbf{x}\_{N};t) - (\mathbf{r}\_{1} + \cdots + \mathbf{r}\_{N})\mathbf{U}(\mathbf{x}\_{1},\cdots,\mathbf{x}\_{N};t) \tag{11}$$

subject to the (*π*, *ν*)-element of

$$\begin{aligned} \mathbf{r}\_{i+1} \mathbf{U}(\mathbf{x}\_{1\prime} \cdots \mathbf{r}\_{i-1}, \mathbf{x}\_i, \mathbf{x}\_i, \mathbf{x}\_{i\prime}, \mathbf{x}\_{i+2\prime} \cdots \mathbf{r}\_{\prime\prime}, \mathbf{x}\_N; t) &= \\ (\mathbf{B}\_i + \mathbf{r}\_i - \mathbf{C}\_i) \mathbf{U}(\mathbf{x}\_{1\prime} \cdots \mathbf{r}\_{i-1}, \mathbf{x}\_i, \mathbf{x}\_i; \mathbf{x}\_i + 1, \mathbf{x}\_{i+2\prime} \cdots \mathbf{r}\_{\prime\prime}, \mathbf{x}\_N; t) \end{aligned} \tag{12}$$

for all *<sup>i</sup>* <sup>=</sup> 1, ··· , *<sup>N</sup>* <sup>−</sup> 1.

## *2.2. Solutions of the Forward Equations via Bethe Ansatz*

The (*π*, *ν*)-element of (11) is

*d dtUπν*(*x*1, ··· , *xN*; *<sup>t</sup>*) = *N* ∑ *j*=1 *<sup>b</sup>π*(*j*)*Uπν*(*x*1, ··· , *xj*−1, *xj* <sup>−</sup> 1, *xj*+1, ··· , *xN*; *<sup>t</sup>*) <sup>−</sup> *<sup>b</sup>π*(1) <sup>+</sup> ··· <sup>+</sup> *<sup>b</sup>π*(*N*) *<sup>U</sup>πν*(*x*1, ··· , *xN*; *<sup>t</sup>*).

Assume the separation of variables to write *<sup>U</sup>πν*(*x*1, ··· , *xN*; *<sup>t</sup>*) = *<sup>U</sup>πν*(*x*1, ··· , *xN*)*T*(*t*). Then, the equation of the spatial variables is

$$\begin{split} \varepsilon \mathbb{I} l\_{\pi \mathsf{V}}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N}) &= b\_{\pi(1)} l\_{\pi \mathsf{V}}(\mathbf{x}\_{1} - 1, \mathbf{x}\_{2}, \cdots, \mathbf{x}\_{N}) + \cdots \\ &+ b\_{\pi(N)} l\_{\pi \mathsf{V}}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N-1}, \mathbf{x}\_{N} - 1) - (b\_{\pi(1)} + \cdots + b\_{\pi(N)}) l\_{\pi \mathsf{V}}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N}) \end{split} \tag{13}$$

for some constant *ε* with respect to *t*, *x*1, ··· , *xN*. Then, we observe that, for any *σ* ∈ S*N*,

$$\prod\_{i=1}^{N} \left( b\_{\pi(i)} \xi\_{\sigma(i)} \right)^{x\_i} \quad \xi\_{1\prime} \cdots \cdot \xi\_N \in \mathbb{C}^N$$

solves (13) if and only if

$$\varepsilon = \frac{1}{\mathcal{J}\_1} + \dots + \frac{1}{\mathcal{J}\_N} - \left( b\_{\pi(1)} + \dots + b\_{\pi(N)} \right). \tag{14}$$

Based on the observation in the above, assume that the matrix **<sup>U</sup>**(*x*1, ··· , *xN*; *<sup>t</sup>*) is invertible and that it is decomposed as

$$\mathbf{U}(\mathbf{x}\_1 \cdot \cdots \cdot, \mathbf{x}\_N; t) = \mathbf{J}(t)\mathbf{U}(\mathbf{x}\_1 \cdot \cdots \cdot, \mathbf{x}\_N)$$

where **J**(*t*) = *Jπν*(*t*) is an *<sup>N</sup><sup>N</sup>* <sup>×</sup> *<sup>N</sup><sup>N</sup>* diagonal matrix where *<sup>J</sup>ππ*(*t*) are functions of time only. Hence, from (11), we obtain

$$\begin{split} \mathbf{J}(\mathbf{J}(t)^{-1}) \left( \frac{d}{dt} \mathbf{J}(t) \right) &= \left( \mathbf{r}\_1 \mathbf{U}(\mathbf{x}\_1 - 1, \mathbf{x}\_{2\prime}, \cdots, \mathbf{x}\_N) + \cdots + \mathbf{r}\_N \mathbf{U}(\mathbf{x}\_1, \cdots, \mathbf{x}\_{N-1}, \mathbf{x}\_N - 1) \right. \\ &\left. - (\mathbf{r}\_1 + \cdots + \mathbf{r}\_N) \mathbf{U}(\mathbf{x}\_1, \cdots, \mathbf{x}\_N) \right) \mathbf{U}(\mathbf{x}\_1, \cdots, \mathbf{x}\_N)^{-1} . \end{split} \tag{15}$$

Both sides of (15) must be a diagonal matrix **E** = *εππ* whose elements are some constants with respect to *t*, *x*1, ··· , *xN*. Thus, we obtain the matrix equation for spatial variables

$$\begin{aligned} \mathbf{E} \mathbf{U}(\mathbf{x}\_1, \dots, \mathbf{x}\_N) &= \mathbf{r}\_1 \mathbf{U}(\mathbf{x}\_1 - 1, \mathbf{x}\_2, \dots, \mathbf{x}\_N) + \dots + \mathbf{r}\_N \mathbf{U}(\mathbf{x}\_1, \dots, \mathbf{x}\_{N-1}, \mathbf{x}\_N - 1) \\ &- (\mathbf{r}\_1 + \dots + \mathbf{r}\_N) \mathbf{U}(\mathbf{x}\_1, \dots, \mathbf{x}\_N) \end{aligned} \tag{16}$$

and the matrix equation for the time variable

$$\frac{d}{dt}\mathbf{J}(t) = \mathbf{J}(t)\mathbf{E} = \mathbf{E}\mathbf{J}(t).$$

**Lemma 2.** *Let* **<sup>D</sup>**(*x*1, ··· , *xN*) = *<sup>D</sup>πν*(*x*1, ··· , *xN*) *be an N<sup>N</sup>* <sup>×</sup> *<sup>N</sup><sup>N</sup> diagonal matrix with*

$$D\_{\pi\pi}(\mathfrak{x}\_1\cdot\cdots\cdot,\mathfrak{x}\_N) = b^{\varkappa\_1}\_{\pi(1)}\cdot\cdots\cdot b^{\varkappa\_N}\_{\pi(N)}.$$

*Then, for any σ* ∈ S*N,*

$$\mathbf{U}(\mathbf{x}\_1, \dots, \mathbf{x}\_N) = \mathbf{D}(\mathbf{x}\_1, \dots, \mathbf{x}\_N) \mathbf{A} \prod\_{i=1}^N \xi\_{\sigma(i)}^{x\_i} \tag{17}$$

*where* **<sup>A</sup>** *is an arbitrary invertible <sup>N</sup><sup>N</sup>* <sup>×</sup> *<sup>N</sup><sup>N</sup> matrix whose elements are constants with respect to x*1, ··· , *xN is a solution of (16) if and only if εππ is given by*

$$\varepsilon\_{\pi\pi} = \frac{1}{\xi\_1} + \dots + \frac{1}{\xi\_N} - \left(b\_{\pi(1)} + \dots + b\_{\pi(N)}\right). \tag{18}$$

**Proof.** First, we observe that

$$\mathbf{r}\_i \mathbf{D}(\mathbf{x}\_1 \cdot \cdots \cdot \mathbf{x}\_{i-1}, \mathbf{x}\_i - 1, \mathbf{x}\_{i+1} \cdot \cdots \cdot \mathbf{x}\_N) = \mathbf{D}(\mathbf{x}\_1 \cdot \cdots \cdot \mathbf{x}\_N)$$

because **<sup>r</sup>***<sup>i</sup>* is a diagonal matrix whose (*π*, *<sup>π</sup>*)-element is *<sup>b</sup>π*(*i*). We observe that

$$\mathbf{D}^{-1}(\mathbf{x}\_1 \cdot \cdots \cdot \mathbf{x}\_N) = \mathbf{D}(\mathbf{x}\_1^{-1} \cdot \cdots \cdot \mathbf{x}\_N^{-1}) \,,$$

and

$$\mathbf{D}(\mathbf{x}\_{1\prime}\cdots\cdot,\mathbf{x}\_{N})\mathbf{D}(y\_{1\prime}\cdots\cdot,y\_{N}) = \mathbf{D}(\mathbf{x}\_{1}+y\_{1\prime}\cdots\cdot,\mathbf{x}\_{N}+y\_{N})\dots$$

Now, we prove the statement. Suppose that (17) is a solution of (16). Substituting (17) into (16) and then dividing both sides by ∏*<sup>N</sup> <sup>i</sup>*=<sup>1</sup> *<sup>ξ</sup>xi σ*(*i*) , then

$$\begin{split} \mathbf{ED}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N}) \mathbf{A} &= \mathbf{r}\_{1} \mathbf{D}(\mathbf{x}\_{1} - 1, \mathbf{x}\_{2}, \cdots, \mathbf{x}\_{N}) \mathbf{A} \mathbf{x}\_{\sigma(1)}^{-1} + \cdots \\ &\quad + \mathbf{r}\_{N} \mathbf{D}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N-1}, \mathbf{x}\_{N} - 1) \mathbf{A} \mathbf{x}\_{\sigma(N)}^{-1} - (\mathbf{r}\_{1} + \cdots + \mathbf{r}\_{N}) \mathbf{D}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N}) \mathbf{A} \\ &= \mathbf{D}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N}) \mathbf{A} \mathbf{x}\_{\sigma(1)}^{-1} + \cdots + \mathbf{D}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N}) \mathbf{A} \mathbf{x}\_{\sigma(N)}^{-1} \\ &\quad - (\mathbf{r}\_{1} + \cdots + \mathbf{r}\_{N}) \mathbf{D}(\mathbf{x}\_{1}, \cdots, \mathbf{x}\_{N}) \mathbf{A}. \end{split}$$

Multiplying by **<sup>A</sup>**−1**D**−1(*x*1, ··· , *xN*) on both sides, we obtain

$$\mathbf{E} = \left(\frac{1}{\xi\_{\nu(1)}} + \dots + \frac{1}{\xi\_{\nu(N)}}\right) \mathbf{I}\_{N^N} - (\mathbf{r}\_1 + \dots + \mathbf{r}\_N)\_{\nu}$$

and thus, the (*π*, *π*)-element of **E** is given by (18). The second part of the proof can be done via the reverse way of the first part of the proof.

The previous lemma implies that the general solution of (16) is given by

$$\mathbf{U}(\mathbf{x}\_{1\prime}\cdot\cdots\prime,\mathbf{x}\_{N}) = \sum\_{\sigma\in\mathcal{S}\_{N}} \mathbf{D}(\mathbf{x}\_{1\prime}\cdot\cdots\prime,\mathbf{x}\_{N}) \mathbf{A}\_{\sigma} \prod\_{i=1}^{N} \mathfrak{I}^{x\_{i}}\_{\sigma(i)}.\tag{19}$$

#### *2.3. Boundary Conditions*

Now, (19) should satisfy the spatial part of the *boundary condition* (12), that is,

$$\begin{aligned} \mathbf{r}\_{i+1} \mathbf{U}(\mathbf{x}\_{1\prime}, \cdots, \mathbf{x}\_{i-1\prime}, \mathbf{x}\_i, \mathbf{x}\_i, \mathbf{x}\_{i+2\prime}, \cdots, \mathbf{x}\_N) &= \\ (\mathbf{B}\_i + \mathbf{r}\_i - \mathbf{C}\_i) \mathbf{U}(\mathbf{x}\_{1\prime}, \cdots, \mathbf{x}\_{i-1\prime}, \mathbf{x}\_i, \mathbf{x}\_i + 1, \mathbf{x}\_{i+2\prime}, \cdots, \mathbf{x}\_N) \end{aligned} \tag{20}$$

for *<sup>i</sup>* <sup>=</sup> 1, ··· , *<sup>N</sup>* <sup>−</sup> 1. Extending the technique used in [4], we will find the formulas of **<sup>A</sup>***<sup>σ</sup>* in (19), which satisfy (20). Define an *N* × *N* diagonal matrix,

$$\mathbf{r}^{\chi\_i} := \begin{bmatrix} b\_1^{\chi\_i} \\ & \ddots \\ & & b\_N^{\chi\_i} \end{bmatrix}$$

and recall the definition of **R***βα* in (1). Then, we observe that

$$\mathbf{R}\_{\mathsf{f}\mathbf{x}} = -\left[\left(\mathbf{I}\_{N^2} - \mathbf{j}\_{\mathsf{f}}^{\mathsf{x}} (\mathbf{B} + \mathbf{r} \otimes \mathbf{I}\_{N} - \mathbf{C})\right) \left(\mathbf{r}^{\mathbf{x}\_{i}} \otimes \mathbf{r}^{\mathbf{x}\_{i} + 1}\right)\right]^{-1} \left[\left(\mathbf{I}\_{N^2} - \mathbf{j}\_{\mathsf{f}}^{\mathsf{x}} (\mathbf{B} + \mathbf{r} \otimes \mathbf{I}\_{N} - \mathbf{C})\right) \left(\mathbf{r}^{\mathbf{x}\_{i}} \otimes \mathbf{r}^{\mathbf{x}\_{i} + 1}\right)\right].$$

**Lemma 3.** *If*

$$\mathbf{A}\_{T\_i \sigma} = \left(\mathbf{I}\_{N^{i-1}} \circledast \mathbf{R}\_{\sigma(i+1)\sigma(i)} \circledast \mathbf{I}\_{N^{N-i-1}}\right) \mathbf{A}\_{\sigma} \tag{21}$$

*for all even permutations <sup>σ</sup> and i* <sup>=</sup> 1, ··· , *<sup>N</sup>* <sup>−</sup> <sup>1</sup>*, then (20) is satisfied.*

**Proof.** First, we note that

$$\begin{split} \mathbf{r}\_{i+1} &= \underbrace{\mathbf{I}\_{N} \otimes \cdots \otimes \mathbf{I}\_{N}}\_{i-1} \otimes \left( \mathbf{I}\_{N} \otimes \mathbf{r} \right) \otimes \underbrace{\mathbf{I}\_{N} \otimes \cdots \otimes \mathbf{I}\_{N}}\_{N-i-1} \\ \mathbf{B}\_{i} + \mathbf{r}\_{i} - \mathbf{C}\_{i} &= \underbrace{\mathbf{I}\_{N} \otimes \cdots \otimes \mathbf{I}\_{N}}\_{i-1} \otimes \left( \mathbf{B} + \mathbf{r} \otimes \mathbf{I}\_{N} - \mathbf{C} \right) \otimes \underbrace{\mathbf{I}\_{N} \otimes \cdots \otimes \mathbf{I}\_{N}}\_{N-i-1} \end{split}$$

and

$$\mathbf{D}(\mathbf{x}\_1, \cdots, \mathbf{x}\_N) = \mathbf{r}^{\mathbf{x}\_1} \otimes \cdots \otimes \mathbf{r}^{\mathbf{x}\_N}.$$

Substituting (19) into (20), we obtain

$$\sum\_{\mathbf{r}\in\mathbb{S}\_{N}} \left[ \left( \mathbf{r}^{\mathbf{r}\_{1}} \otimes \dots \otimes \mathbf{r}^{\mathbf{r}\_{i-1}} \otimes \left( \mathbf{r}^{\mathbf{r}\_{i}} \otimes \mathbf{r}^{\mathbf{r}\_{i}+1} \right) \otimes \mathbf{r}^{\mathbf{r}\_{i+2}} \otimes \dots \otimes \mathbf{r}^{\mathbf{r}\_{N}} \right) - \right. $$

$$\left( \mathbf{r}^{\mathbf{r}\_{1}} \otimes \dots \otimes \mathbf{r}^{\mathbf{r}\_{i-1}} \otimes \left( \mathbf{B} + \mathbf{r} \otimes \mathbf{I}\_{N} - \mathbf{C} \right) \left( \mathbf{r}^{\mathbf{r}\_{i}} \otimes \mathbf{r}^{\mathbf{r}\_{i}+1} \right) \otimes \left. \mathbf{r}^{\mathbf{r}\_{i+2}} \otimes \dots \otimes \mathbf{r}^{\mathbf{r}\_{N}} \right) \xi\_{\sigma(i+1)} \right] \tag{22}$$

$$\times \mathbf{A}\_{\sigma} \xi\_{\sigma(i)}^{\mathbf{x}\_{i}} \xi\_{\sigma(i+1)}^{\mathbf{x}\_{i}} \prod\_{j \neq i, i+1} \xi\_{\sigma(j)}^{\mathbf{x}\_{j}} = \mathbf{0}.$$

If we express (22) as a sum over the alternating group A*<sup>N</sup>*

∑ *σ*∈A*<sup>N</sup>* **<sup>r</sup>***x*<sup>1</sup> ⊗ ··· ⊗ **<sup>r</sup>***xi*−<sup>1</sup> <sup>⊗</sup> **<sup>r</sup>***xi* <sup>⊗</sup> **<sup>r</sup>***xi*+1 <sup>⊗</sup> **<sup>r</sup>***xi*+<sup>2</sup> ⊗···⊗ **<sup>r</sup>***xN* − **<sup>r</sup>***x*<sup>1</sup> ⊗ ··· ⊗ **<sup>r</sup>***xi*−<sup>1</sup> <sup>⊗</sup> **<sup>B</sup>** <sup>+</sup> **<sup>r</sup>** <sup>⊗</sup> **<sup>I</sup>***<sup>N</sup>* <sup>−</sup> **<sup>C</sup> r***xi* <sup>⊗</sup> **<sup>r</sup>***xi*+1 <sup>⊗</sup> **<sup>r</sup>***xi*+<sup>2</sup> ⊗···⊗ **<sup>r</sup>***xN ξσ*(*i*+1) **A***σ* + **<sup>r</sup>***x*<sup>1</sup> ⊗ ··· ⊗ **<sup>r</sup>***xi*−<sup>1</sup> <sup>⊗</sup> **<sup>r</sup>***xi* <sup>⊗</sup> **<sup>r</sup>***xi*+1 <sup>⊗</sup> **<sup>r</sup>***xi*+<sup>2</sup> ⊗···⊗ **<sup>r</sup>***xN* − **<sup>r</sup>***x*<sup>1</sup> ⊗ ··· ⊗ **<sup>r</sup>***xi*−<sup>1</sup> <sup>⊗</sup> **<sup>B</sup>** <sup>+</sup> **<sup>r</sup>** <sup>⊗</sup> **<sup>I</sup>***<sup>N</sup>* <sup>−</sup> **<sup>C</sup> r***xi* <sup>⊗</sup> **<sup>r</sup>***xi*+1 <sup>⊗</sup> **<sup>r</sup>***xi*+<sup>2</sup> ⊗···⊗ **<sup>r</sup>***xN ξσ*(*i*) **A***Ti<sup>σ</sup>* <sup>×</sup> *<sup>ξ</sup>xi σ*(*i*) *ξxi <sup>σ</sup>*(*i*+1) ∏ *<sup>j</sup>* <sup>=</sup>*i*, *<sup>i</sup>*+<sup>1</sup> *ξ xj <sup>σ</sup>*(*j*) <sup>=</sup> **<sup>0</sup>**. (23)

However, (23) is satisfied if

 **<sup>r</sup>***x*<sup>1</sup> ⊗ ··· ⊗ **<sup>r</sup>***xi*−<sup>1</sup> <sup>⊗</sup> **<sup>r</sup>***xi* <sup>⊗</sup> **<sup>r</sup>***xi*+1 <sup>⊗</sup> **<sup>r</sup>***xi*+<sup>2</sup> ⊗···⊗ **<sup>r</sup>***xN* − **<sup>r</sup>***x*<sup>1</sup> ⊗ ··· ⊗ **<sup>r</sup>***xi*−<sup>1</sup> <sup>⊗</sup> **<sup>B</sup>** <sup>+</sup> **<sup>r</sup>** <sup>⊗</sup> **<sup>I</sup>***<sup>N</sup>* <sup>−</sup> **<sup>C</sup> r***xi* <sup>⊗</sup> **<sup>r</sup>***xi*+1 <sup>⊗</sup> **<sup>r</sup>***xi*+<sup>2</sup> ⊗···⊗ **<sup>r</sup>***xN ξσ*(*i*+1) **A***σ* <sup>=</sup> <sup>−</sup> **<sup>r</sup>***x*<sup>1</sup> ⊗ ··· ⊗ **<sup>r</sup>***xi*−<sup>1</sup> <sup>⊗</sup> **<sup>r</sup>***xi* <sup>⊗</sup> **<sup>r</sup>***xi*+1 <sup>⊗</sup> **<sup>r</sup>***xi*+<sup>2</sup> ⊗···⊗ **<sup>r</sup>***xN* − **<sup>r</sup>***x*<sup>1</sup> ⊗ ··· ⊗ **<sup>r</sup>***xi*−<sup>1</sup> <sup>⊗</sup> **<sup>B</sup>** <sup>+</sup> **<sup>r</sup>** <sup>⊗</sup> **<sup>I</sup>***<sup>N</sup>* <sup>−</sup> **<sup>C</sup> r***xi* <sup>⊗</sup> **<sup>r</sup>***xi*+1 <sup>⊗</sup> **<sup>r</sup>***xi*+<sup>2</sup> ⊗···⊗ **<sup>r</sup>***xN ξσ*(*i*) **A***Ti<sup>σ</sup>*

for each even permutation *σ*, which is equivalent to (21).

In fact, (3) and (5) implies that the assumption of Lemma 3 is satisfied, hence (19) with **A***σ* in (5) satisfies (20).

#### *2.4. Consistency Relations—Proof of Lemma 1*

Lemma 1 confirms that the multi-species TASEP is *integrable* even when the rates are species-dependent.

2.4.1. Proof of Lemma 1 (a)

It suffices to show that

$$(\mathbf{R}\_{\mathbf{a}\notin} \otimes \mathbf{I}\_N \otimes \mathbf{I}\_N)(\mathbf{I}\_N \otimes \mathbf{I}\_N \otimes \mathbf{R}\_{\gamma\delta}) = (\mathbf{I}\_N \otimes \mathbf{I}\_N \otimes \mathbf{R}\_{\gamma\delta})(\mathbf{R}\_{\mathbf{a}\notin} \otimes \mathbf{I}\_N \otimes \mathbf{I}\_N).$$

This equality clearly holds because both sides are equal to **R***αβ* ⊗ **R***γδ*.

2.4.2. Proof of Lemma 1 (b)—Yang-Baxter Equation

It suffices to show

$$(\mathbf{R}\_{\gamma\beta}\otimes\mathbf{I}\_{N})(\mathbf{I}\_{N}\otimes\mathbf{R}\_{\gamma\mathfrak{a}})(\mathbf{R}\_{\text{fix}}\otimes\mathbf{I}\_{N})=(\mathbf{I}\_{N}\otimes\mathbf{R}\_{\text{fix}})(\mathbf{R}\_{\gamma\mathfrak{a}}\otimes\mathbf{I}\_{N})(\mathbf{I}\_{N}\otimes\mathbf{R}\_{\gamma\mathfrak{f}}).\tag{24}$$

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

If we re-arrange the columns and the rows of the *<sup>N</sup>*<sup>3</sup> <sup>×</sup> *<sup>N</sup>*<sup>3</sup> matrices in (24) so that all their labels from the same multi-set [*i*, *j*, *k*] are grouped together, then the matrices in (24) become block-diagonal (See Figure 1).

**Figure 1.** Form of the *<sup>N</sup>*<sup>3</sup> <sup>×</sup> *<sup>N</sup>*<sup>3</sup> matrices in (24) after re-ordering the rows and columns.

Let (**<sup>R</sup>** <sup>⊗</sup> **<sup>I</sup>***N*)[*i*,*j*,*k*] be the sub-matrix of (**<sup>R</sup>** <sup>⊗</sup> **<sup>I</sup>***N*) whose rows and columns are labelled by the permutations of the multi-set [*i*, *<sup>j</sup>*, *<sup>k</sup>*], and similarly, we define (**I***<sup>N</sup>* <sup>⊗</sup> **<sup>R</sup>**)[*i*,*j*,*k*]. Then, in order to show (24), it suffices to show

$$\begin{aligned} \left(\mathbf{R}\_{\gamma\notin}\circledast\mathbf{I}\_{N}\right)\_{[i,j,k]}(\mathbf{I}\_{N}\circledast\mathbf{R}\_{\gamma a})\_{[i,j,k]}(\mathbf{R}\_{\not\models a}\circledast\mathbf{I}\_{N})\_{[i,j,k]} &= \\ (\mathbf{I}\_{N}\circledast\mathbf{R}\_{\not\models a})\_{[i,j,k]}(\mathbf{R}\_{\gamma a}\circledast\mathbf{I}\_{N})\_{[i,j,k]}(\mathbf{I}\_{N}\circledast\mathbf{R}\_{\gamma\not\models \not k})\_{[i,j,k]} \end{aligned} \tag{25}$$

for each multi-set [*i*, *<sup>j</sup>*, *<sup>k</sup>*] whose elements are from {1, 2, 3} because all matrices are blockdiagonal matrices in the same form. If *i* = *j* = *k*, (25) is equivalent to

$$S\_{\gamma\beta}(i)S\_{\gamma\alpha}(i)S\_{\beta\alpha}(i) = S\_{\beta\alpha}(i)S\_{\gamma\alpha}(i)S\_{\gamma\beta}(i),$$

which is trivially true. If *i* = *j* > *k*, then (25) is equivalent to

$$
\begin{split}
&\left[\begin{array}{ccc}
S\_{\gamma\beta}(k) & T\_{\gamma\beta}(k) & 0 \\
0 & -1 & 0 \\
0 & 0 & S\_{\gamma\beta}(i)
\end{array}\right] \left[\begin{array}{ccc}
S\_{\gamma a}(i) & 0 & 0 \\
0 & S\_{\gamma a}(k) & T\_{\gamma a}(k) \\
0 & 0 & -1
\end{array}\right] \left[\begin{array}{ccc}
S\_{\beta a}(k) & T\_{\beta a}(k) & 0 \\
0 & -1 & 0 \\
0 & 0 & S\_{\beta a}(i)
\end{array}\right] \\
&= \left[\begin{array}{ccc}
S\_{\beta a}(i) & 0 & 0 \\
0 & S\_{\beta a}(k) & T\_{\beta a}(k) \\
0 & 0 & -1
\end{array}\right] \left[\begin{array}{ccc}
S\_{\gamma a}(k) & T\_{\gamma a}(k) & 0 \\
0 & -1 & 0 \\
0 & 0 & S\_{\gamma a}(i)
\end{array}\right] \left[\begin{array}{ccc}
S\_{\gamma\beta}(i) & 0 & 0 \\
0 & S\_{\gamma\beta}(k) & T\_{\gamma\beta}(k) \\
0 & 0 & -1
\end{array}\right],
\end{split}
$$

which can be easily verified by direct computation. Similarly, the other two cases of (25) for *i* = *j* < *k* and for the case that all *i*, *j*, *k* are distinct can be verified by direct computation.

#### 2.4.3. Proof of Lemma 1 (c)

It suffices to show that **<sup>R</sup>***βα***R***αβ* <sup>=</sup> **<sup>I</sup>***N*<sup>2</sup> . Let us re-arrange the rows and the columns in the same way as in the proof of Lemma 1 (b) to make **R***βα* and **R***αβ* block-diagonal. Then, each block on the diagonal of **R***βα* is either a 1 × 1 matrix or a 2 × 2 matrix. The 1 × 1 sub-matrix of **R***βα* consisting of the row *ii* and the column *ii* is *Sβα*(*i*) . The 2 × 2 sub-matrix of **R***βα* consisting of the rows *ij*, *ji* and the columns *ij*, *ji* with *i* < *j* is

$$
\begin{bmatrix}
\mathcal{S}\_{\beta\alpha}(i) & T\_{\beta\alpha}(i) \\
0 & -1
\end{bmatrix}.
$$

Similarly, the 1 × 1 sub-matrix of **R***αβ* consisting of the row *ii* and the column *ii* is *Sαβ*(*i*) , and the 2 × 2 sub-matrix of **R***αβ* consisting of the rows *ij*, *ji* and the columns *ij*, *ji* with *i* < *j* is

$$
\begin{bmatrix}
\mathcal{S}\_{\alpha\beta}(\vec{r}) & T\_{\alpha\beta}(\vec{r}) \\
0 & -1
\end{bmatrix}
$$

.

It is trivial that *Sβα*(*i*)*Sαβ*(*i*) = 1, and this can be verified,

$$
\begin{bmatrix} S\_{\beta\alpha}(i) & T\_{\beta\alpha}(i) \\ 0 & -1 \end{bmatrix} \begin{bmatrix} S\_{a\beta}(i) & T\_{a\beta}(i) \\ 0 & -1 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \prime
$$

by the direct computation.

#### *2.5. Initial Condition*

The contour integral of (19) multiplied by **J**(*t*)from the left and by **D** (*y*1, ··· , *yN*) <sup>∏</sup>*<sup>i</sup> <sup>ξ</sup>* −*yi*−1 *i* from the right, that is, the right-hand side of (6) still satisfies (11) and (12). (The contour is the one introduced in Theorem 1). Hence, it remains to show that all transition probabilities *<sup>P</sup>*(*Y*,*ν*)(*X*, *<sup>π</sup>*; *<sup>t</sup>*) satisfy the initial condition

$$P\_{\{Y,\nu\}}(X,\pi;0) \quad = \sum\_{\sigma \in \mathcal{S}\_N} \int\_{\mathcal{C}} \cdots \int\_{\mathcal{C}} A\_{\sigma}^{\pi \nu} \prod\_{i=1}^N b\_{\pi(i)}^{x\_i} b\_{\nu(i)}^{-y\_i} \prod\_{i=1}^N \left(\xi\_{\sigma(i)}^{x\_i - y\_{\sigma(i)} - 1}\right) d\xi\_1 \cdots d\xi\_N$$
 
$$= \begin{cases} 1 & \text{if } (Y,\nu) = (X,\pi); \\ 0 & \text{otherwise} \end{cases} \tag{26}$$

where *Aπν <sup>σ</sup>* is the (*π*, *<sup>ν</sup>*)-element of **<sup>A</sup>***σ*. We will show that the integral with the identity permutation in the sum satisfies (26), and other integral terms with non-identity permutations are zero.

**Proof.** If *<sup>σ</sup>* is the identity permutation, then **<sup>A</sup>***<sup>σ</sup>* is the identity matrix. Hence, if *<sup>π</sup>* <sup>=</sup> *<sup>ν</sup>*, then the integral is zero. It is easy to see that if *<sup>π</sup>* = *<sup>ν</sup>* and *xi* = *yi* for all *<sup>i</sup>*, then the integral is 1. If *<sup>π</sup>* = *<sup>ν</sup>* and *xi* > *yi* for some *<sup>i</sup>* (recall that our model is totally asymmetric), then the integral becomes zero when integrating with respect to *ξi*. Now, suppose that *σ* is not the identity permutation. Note that the factors in *Aπν <sup>σ</sup>* are from (1), all poles arising from *Aπν <sup>σ</sup>* , if any, are outside the contours. There exists an *i* such that *xi* − *yσ*(*i*) − 1 ≥ 0 because each *xi* ≥ *yi* and *σ* is not the identity permutation. Hence, integrating with respect to *i*, the integral is 0.

#### **3. Conclusions**

In this paper, we have shown that the Bethe ansatz method is still applicable to the multi-species TASEP with species-dependent rates. Theorem 1 provides the transition probabilities for all possible compositions of species, which is expected to be used to study further objects, such as the current distribution for certain special initial configurations. The methods used in this paper have limitations in extending to the ASEP (0 < *p* < 1) for now, but it would be interesting to see if the methods can be used to study the speciesinhomogeneity of other multi-species models.

**Funding:** This research was funded by Nazarbayev University (Faculty Development Competitive Grant grant number 090118FD5341 and 021220FD4251.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The author is grateful to anonymous referees for providing valuable comments.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


## *Article* **Shock Structure and Relaxation in the Multi-Component Mixture of Euler Fluids**

**Damir Madjarevi´c 1,†, Milana Pavi´c-Coli´ ˇ c2,† and Srboljub Simi´c 2,\*,†**


**Abstract:** The shock structure problem is studied for a multi-component mixture of Euler fluids described by the hyperbolic system of balance laws. The model is developed in the framework of extended thermodynamics. Thanks to the equivalence with the kinetic theory approach, phenomenological coefficients are computed from the linearized weak form of the collision operator. Shock structure is analyzed for a three-component mixture of polyatomic gases, and for various combinations of parameters of the model (Mach number, equilibrium concentrations and molecular mass ratios). The analysis revealed that three-component mixtures possess distinguishing features different from the binary ones, and that certain behavior may be attributed to polyatomic structure of the constituents. The multi-temperature model is compared with a single-temperature one, and the difference between the mean temperatures of the mixture are computed. Mechanical and thermal relaxation times are computed along the shock profiles, and revealed that the thermal ones are smaller in the case discussed in this study.

**Keywords:** shock structure; extended thermodynamics; kinetic theory of gases

#### **1. Introduction**

A shock wave is a paradigm of non-equilibrium process with irreversible character. It drives the system out of equilibrium and, at the same time, causes irreversible changes observed through the increase of entropy [1–3]. In idealized models, shock waves are described as singular surfaces on which jumps of field variables occur. However, when dissipation is taken into account, shock waves are smeared out into continuous profiles, thus creating the so-called shock structure. [4].

Analysis of the shock structure in gaseous media is strongly influenced by the model taken to describe dissipation. For the present study, two approaches are significant continuum (macroscopic) and kinetic (mesoscopic) one. Continuum approach uses macroscopic field variables and their governing equations to describe the system. Within it, classical continuum models assume non-local constitutive relations for non-convective fluxes like stress tensor, heat flux, or diffusion flux. Typical examples are Navier–Stokes model, Fourier law and Fick law [5]. Such models usually lead to parabolic partial differential equations, or their generalizations. In the same (continuum) setting, dissipation may also be described by means of relaxation, which assumes that non-equilibrium variables converge (relax) towards their equilibrium values during the process. They yield models which usually consist of the hyperbolic systems of balance laws. It is remarkable that in the classical (thermodynamic) limit, i.e., in the small relaxation time approximation, one may recover the classical parabolic models starting from the hyperbolic ones [6–9]. In more physical terms, classical models can be regarded as approximations of the hyperbolic ones when processes occur in the neighborhood of the local equilibrium state [10].

**Citation:** Madjarevi´c, D.; Pavi´c-Coli´ ˇ c, M.; Simi´c, S. Shock Structure and Relaxation in the Multi-Component Mixture of Euler Fluids. *Symmetry* **2021**, *13*, 955. https://doi.org/ 10.3390/sym13060955

Academic Editors: Rahmat Ellahi and Iver H. Brevik

Received: 14 April 2021 Accepted: 19 May 2021 Published: 27 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The kinetic approach is focused on the smaller (meso) scale. It takes statistical description—single particle velocity distribution function—and mutual interaction between the particles which causes its change. It is built around the Boltzmann equation and inherits dissipation by construction, through the collision operator [11]. A remarkable feature of the kinetic theory is that it can recover classical continuum models in the so-called hydrodynamic limit, when the Knudsen number is small, using the Chapman-Enskog asymptotic expansion. On the other hand, by means of the moment method kinetic theory yields hyperbolic systems of balance laws for macroscopic observables, taken as moments of the velocity distribution function. Under appropriate assumptions, relaxation continuum models of rarefied gases are completely equivalent to moment equations of the kinetic theory of gases [12,13].

Either approach described above has its advantages and shortcomings. For example, kinetic theory is obviously more refined than the continuum one, but evolution of the system is eventually monitored through macroscopic observables. In this study we shall try to take the benefits from both approaches and analyze the shock structure in the hyperbolic model of the multi-component mixture of Euler fluids. By an Euler fluid we consider the fluid in which viscous and heat conducting effects are negligible. By a phrase multicomponent mixture, although it seems to be a pleonasm, we reveal our intention to analyze the mixture of more than two constituents (as opposed to binary, i.e., two-component mixture). In fact, we shall analyze the shock structure in a three-component mixture as a paradigmatic example of the multi-component one. It is known that a multi-component mixture exhibits peculiar phenomena with respect to the binary one, as for instance the uphill diffusion [14–16]. Before we embark on the shock structure analysis, a short review of modelling issues and important results in this particular problem relevant to our study will be given.

Our analysis will be based upon the continuum model of mixture developed within the framework of rational extended thermodynamics (RET) [13,17] (for a review of the model one may consult [18]). This model is rather peculiar from macroscopic point of view since it takes velocities and temperatures of the constituents as field variables. RET formalism, i.e., Galilean invariance of the governing equations and their compatibility with the entropy inequality, makes the model thermodynamically consistent. The model is hyperbolic, and dissipation is taken into account through mutual interaction between the constituents, which turns out to be of the relaxation type. However, this approach is limited in the sense that phenomenological coefficients, which appear in the closure process, cannot be explicitly determined within its framework. To resolve this problem, we turn to the kinetic approach based upon the polyatomic gas mixture model, i.e., the system of Boltzmann equations, with the continuous internal energy variable [19–21], which is more tractable than the semi-classical model with discrete energy levels [22–25], yet shown to reproduce physical requirements [26–29] and to be well posed mathematically in the case of space homogeneity [30]. Moreover, the multi-temperature (MT) model for mixture of Euler fluids is fully compatible with moment equations derived from this system at an Euler level [31,32]. In particular, we find the complete closure of the MT model by using the phenomenological coefficients from the weak form of the continuous collision operator at an Euler level [33].

In this study, the shock wave will be treated as a continuous traveling wave solution of the governing equations which converges to equilibrium states at infinity. Such an assumption makes feasible the use of dynamical systems approach. This kind of theoretical shock structure analysis in continuum models of single-component fluids has its roots in [34], where it was applied to Navier–Stokes–Fourier model. Further application of this method to hyperbolic continuum models was studied in [35] (see also Chapter 12 in [13] for a comprehensive review). The main feature of this approach is that the shock structure is described by a heteroclinic orbit in the state space which connects stationary points related to upstream and downstream equilibrium. More precisely, this orbit is an intersection of unstable and stable manifolds of these stationary points [36], and can be related to

the bifurcation of the upstream equilibrium point [37]. This could also be useful for the numerical computation of the shock profiles, albeit in the limited way. Namely, it was shown [38] that in hyperbolic systems of balance laws, endowed with a convex entropy density, continuous shock structure does not exist if the shock speed is greater than the highest characteristic speed of the system; instead, there appears a continuous solution with the so-called sub-shock.

Apart from single-component models of gases, shock waves were also studied in the context of mixtures. In the early stage, they were analyzed by means of classical models [39] and DSMC [40]. Efficient computing devices moved the focus towards development of numerical methods for the solution of Boltzmann equations for mixtures [41–43]. On the other hand, shock structure was also analyzed for the moment equations of mixtures, derived from the kinetic theory for various types of interaction [44,45]. At the same time, detailed and systematic analysis of the shock structure in hyperbolic MT mixture of Euler fluids was given in [46,47], which revealed certain new features of the temperature overshoot. In this context, existence of continuous profiles and appearance of sub-shocks turned out to be more delicate than in the one-component case, as discussed in [48–50]. Shock structure in the mixture of noble gases was revisited by DSMC in [51]. Recently, attention was given to the entropy growth and the entropy production analysis in the mixture of Euler fluids [52], and the mixture with viscous and thermal dissipation [53]. It was found out that in former case diffusion makes the major contribution to dissipation, while in the latter one it is the thermal conduction.

It is remarkable that most of the studies mentioned above were concerned with binary mixtures. They certainly reveal the properties which distinguish them from singlecomponent fluids, like temperature overshoot or influence of diffusion on the shock thickness. However, there are phenomena which can be observed only when the multicomponent mixture is studied. This gives us strong motivation to analyze the shock structure problem in three-component mixture of Euler fluids. So far, only few references were found which analyze the shock waves in multi-component mixtures within the frameworks indicated above. In [54], the shock structure was analyzed in a mixture of gases with disparate molecular masses, and viscous and thermal dissipation; ternary mixture is studied as a final example, with a special attention on the presence of two zones within the shock profile. In [55], the analysis was based upon numerical solution of the Boltzmann equations; although the main focus was on binary mixture, three-species mixture was studied for the computation of parallel temperatures which exhibit an overshoot. Finally, in [56], three- and four-component mixtures are studied by numerical solution of the Boltzmann equations; profiles of macroscopic variables are found to have good agreement with experiments and computations by other authors.

Our study will be focused on the shock waves in mixtures of Euler fluids, using the MT model developed in [17], together with phenomenological coefficients from [33]. It will be assumed that the mixture is homogeneous in the sense that in each infinitesimal volume element all the constituents are present. It will also be assumed that all the constituents are ideal gases, described by usual thermal and caloric equations of state, and that chemical reactions do not occur between them. Once equipped with the closed, thermodynamically consistent system of governing equations, we shall analyze numerically computed shock structures in three-component mixtures. Our aim is to reveal and bring to light possible qualitative novelties which appear when the mixture contains more than two components.

In the rest of the paper, model of the mixture of Euler fluids is presented in Section 2. It brings together continuum modeling and kinetic theory results to get the closed system of equations. In Section 3, the shock structure problem for a three-component mixture is formulated. It contains all the information needed for derivation of equations, their transformation into dimensionless form, and building up an appropriate numerical procedure. Section 4 contains numerical solutions of the shock structure equations and their analysis. Shock profiles are analyzed for different values of parameters of the model (Mach number, equilibrium concentrations and molecular mass ratios), which provides an insight of their

importance in the shock structure analysis. Additionally, multi-temperature shock profiles are compared with single-temperature ones, and the norm of their difference is computed. Finally, Section 5 brings the analysis of relaxation processes captured by the model. It is of utmost importance since they manifest dissipation. Therefore, relative magnitudes of relaxation times for mechanical and thermal relaxations are computed throughout the profile. Paper ends up with proper conclusions and an outlook of future research directions.

#### **2. Multi-Component Mixture of Euler Fluids**

A mixture is a medium consisted of more than one identifiable constituent. One of the most comprehensive mathematical models for multi-component mixtures of Euler fluids is developed within the framework of rational extended thermodynamics (RET) [13]. It inherits the postulates, so-called metaphysical principles, which are the basis for rational thermodynamics of mixtures [57], but also exploits the principles of RET—local dependence of constitutive functions on field variables, Galilean invariance of governing equations and entropy principle with the convex entropy density. As a result, governing equations have the structure of hyperbolic system of balance laws. For a complete resolution of the closure problem, it is important to note that these models are fully compatible with transfer equations for moments, obtained from Boltzmann equation(s) by means of the moment method.

#### *2.1. General Structure of the Model*

In the RET approach to the mixture modelling it is assumed that the state of each constituent *α* = 1, ... , *n* is described by its own set of field variables—mass density *ρα*, velocity **v***α* and temperature *Tα*. There lies the main difference with respect to other, more traditional approaches—it is (primarily) multi-velocity, and multi-temperature model. Further assumptions constitute the set of postulates known as metaphysical principles [57]: (i) all properties of the mixture are mathematical consequences of properties of the constituents; (ii) behavior of each constituent is governed by the same equations as if it were isolated from the mixture, but mutual interactions with other constituents must be taken into account; (iii) behavior of the whole mixture is governed by the same equations as is a single body.

Application of postulate (ii) yields governing equations for the constituents:

$$\begin{aligned} \frac{\partial \rho\_{a}}{\partial t} + \text{div}(\rho\_{a} \mathbf{v}\_{a}) &= \mathbf{r}\_{a}, \\ \frac{\partial}{\partial t}(\rho\_{a} \mathbf{v}\_{a}) + \text{div}(\rho\_{a} \mathbf{v}\_{a} \otimes \mathbf{v}\_{a} - \mathbf{t}\_{a}) &= \mathbf{m}\_{a}, \\ \frac{\partial}{\partial t} \left(\frac{1}{2} \rho\_{a} \mathbf{v}\_{a}^{2} + \rho\_{a} \boldsymbol{\varepsilon}\_{a}\right) + \text{div}\left\{ \left(\frac{1}{2} \rho\_{a} \mathbf{v}\_{a}^{2} + \rho\_{a} \boldsymbol{\varepsilon}\_{a}\right) \mathbf{v}\_{a} - \mathbf{t}\_{a} \mathbf{v}\_{a} + \mathbf{q}\_{a} \right\} &= \mathbf{c}\_{a}, \end{aligned} \tag{1}$$

where *εα* are the specific internal energies, **t***<sup>α</sup>* are the stress tensors, **q***<sup>α</sup>* are the heat fluxes and *τα*, **m***<sup>α</sup>* and *e<sup>α</sup>* are the source terms describing interaction between the constituents. To apply postulate (iii) we must introduce restriction (axiom) regarding the source terms:

$$\sum\_{a=1}^{n} \pi\_{a} = 0, \quad \sum\_{a=1}^{n} \mathbf{m}\_{a} = \mathbf{0}, \quad \sum\_{a=1}^{n} e\_{a} = 0. \tag{2}$$

Then, summation of Equation (1) yields the conservation laws for the mixture:

$$\begin{cases} \frac{\partial \rho}{\partial t} + \text{div}(\rho \mathbf{v}) = 0, \\ \frac{\partial}{\partial t}(\rho \mathbf{v}) + \text{div}(\rho \mathbf{v} \otimes \mathbf{v} - \mathbf{t}) = \mathbf{0}, \\ \frac{\partial}{\partial t} \left( \frac{1}{2} \rho v^2 + \rho \varepsilon \right) + \text{div} \left\{ \left( \frac{1}{2} \rho v^2 + \rho \varepsilon \right) \mathbf{v} - \mathbf{t} \mathbf{v} + \mathbf{q} \right\} = 0, \end{cases} \tag{3}$$

provided we define the mixture variables as follows:

$$\begin{aligned} \boldsymbol{\rho} &= \sum\_{a=1}^{n} \rho\_{a\prime} & \mathbf{v} &= \frac{1}{\rho} \sum\_{a=1}^{n} \rho\_{a} \mathbf{v}\_{a\prime} & \mathbf{u}\_{a} &= \mathbf{v}\_{a} - \mathbf{v}\_{\prime} \\ \boldsymbol{\varepsilon}\_{I} &= \frac{1}{\rho} \sum\_{a=1}^{n} \rho\_{a} \boldsymbol{\varepsilon}\_{a\prime} & \boldsymbol{\varepsilon} &= \boldsymbol{\varepsilon}\_{I} + \frac{1}{2\rho} \sum\_{a=1}^{n} \rho\_{a} \mathbf{u}\_{a\prime}^{2} \\ \mathbf{t} &= \sum\_{a=1}^{n} (\mathbf{t}\_{a} - \rho\_{a} \mathbf{u}\_{a} \otimes \mathbf{u}\_{a}) & \mathbf{q} &= \sum\_{a=1}^{n} \left\{ \mathbf{q}\_{a} + \rho\_{a} \left( \boldsymbol{\varepsilon}\_{a} + \frac{1}{2} \boldsymbol{u}\_{a}^{2} \right) \mathbf{u}\_{a} - \mathbf{t}\_{a} \mathbf{u}\_{a} \right\}, \end{aligned} \tag{4}$$

where *ρ* is the mixture mass density, **v** is the mean velocity, **u***<sup>α</sup>* is the diffusion velocity, *ε<sup>I</sup>* is the intrinsic specific internal energy, *ε* is the specific internal energy, **t** is the mixture stress tensor, and **q** is the mixture flux of internal energy. Definitions (4) comprise the postulate (i).

In the sequel we shall introduce further restrictions into our model. First, we shall assume that all the constituents are Euler fluids, i.e., their stress tensors are symmetric and heat fluxes vanish:

$$\mathbf{t}\_{\mathfrak{k}} = -p\_a \mathbf{I}\_{\prime} \quad \mathbf{q}\_{\mathfrak{a}} = \mathbf{0},\tag{5}$$

where *pα* are partial pressures, and **I** is the identity tensor. Second, we shall assume that each constituent obey thermal and caloric equation of state of ideal gas:

$$p\_a = \rho\_a \frac{k\_\text{B}}{m\_a} T\_{a\prime} \quad \varepsilon\_a = \frac{k\_\text{B}}{m\_a(\gamma\_a - 1)} T\_a = \varepsilon\_{Va} T\_{a\prime} \tag{6}$$

where *k*<sup>B</sup> is the Boltzmann constant, and *m<sup>α</sup>* are the molecular masses of the constituents. Ratios of specific heats *γα* and constant volume specific heats *cV<sup>α</sup>* of the constituents are assumed to be constant. Finally, we shall assume that there are no mutual chemical interactions between the constituents, which amounts to:

$$
\pi\_a = 0, \quad \pi = 1, \dots, n. \tag{7}
$$

To describe behavior of the mixture, 5*n* scalar field variables are used. Thus, 5*n* scalar governing equations are needed, and this choice may not be unique. Here, we follow the standard procedure and choose conservation laws (3) for the mixture, and *n* − 1 balance laws (1) (we drop the balance laws for the constituent *<sup>α</sup>* <sup>=</sup> *<sup>n</sup>*, and use index *<sup>b</sup>* <sup>=</sup> 1, ... , *<sup>n</sup>* <sup>−</sup> <sup>1</sup> instead of *α*). Such a model is usually called the multi-temperature (MT) model of mixture.

The RET approach is peculiar in the sense that it requires Galilean invariance of governing equations from the outset, and the compatibility with the entropy inequality with convex entropy density. This resolves the closure problem to an extent standard for macroscopic theories. To summarize the consequences of these principles, primarily reflected on the source terms, we emphasize that Galilean invariance restricts their velocity dependence:

$$\mathbf{m}\_b = \mathbf{\hat{n}}\_b; \quad \varepsilon\_b = \mathbb{\ell}\_b + \mathbf{\hat{n}}\_b \cdot \mathbf{v}, \quad b = 1, \dots, n - 1,\tag{8}$$

while the entropy principle yields more specific structure of velocity independent terms **m**ˆ *<sup>b</sup>* and *e*ˆ*b*:

$$\mathbf{m}\_{\mathbb{b}} = -\sum\_{c=1}^{n-1} \psi\_{\mathbb{b}c}(\mathbf{w}) \left( \frac{\mathbf{u}\_{\mathbb{c}}}{T\_{\mathbb{c}}} - \frac{\mathbf{u}\_{\mathbb{n}}}{T\_{\mathbb{n}}} \right); \quad \mathbb{\mathcal{E}}\_{\mathbb{b}} = -\sum\_{c=1}^{n-1} \theta\_{\mathbb{b}c}(\mathbf{w}) \left( -\frac{1}{T\_{\mathbb{c}}} + \frac{1}{T\_{\mathbb{n}}} \right), \tag{9}$$

where *<sup>ψ</sup>bc*(**w**) and *<sup>θ</sup>bc*(**w**) are positive semi-definite matrix functions of objective quantities **w** [17,18].

If the conditions of thermodynamic process are such that large discrepancies between species temperatures are not expected, one may use the single-temperature (ST), but still multi-velocity model. All the constituents have the same temperature *T*, and governing equations consist of the conservation laws (3) and the balance laws (1)1,2, i.e., energy balance laws for the constituents are discarded. This system is a principal subsystem in the hierarchy of hyperbolic systems, as described in [17,58].

Local equilibrium conditions **<sup>m</sup>***<sup>b</sup>* = **<sup>0</sup>**, *eb* = 0 imply the following conditions:

$$\begin{aligned} T\_1 = T\_2 = \dots = T\_n = T; \\ \mathbf{u}\_1 = \mathbf{u}\_2 = \dots = \mathbf{u}\_{\mathbb{R}} (= \mathbf{0}) \quad \Leftrightarrow \quad \mathbf{v}\_1 = \mathbf{v}\_2 = \dots = \mathbf{v}\_{\mathbb{R}} = \mathbf{v}\_{\mathbb{A}} \end{aligned} \tag{10}$$

i.e., all the constituents have common temperature, and common velocity. Conditions (10) determine the so-called equilibrium subsystem [58].

A brief comment on source terms may shed a light on their role in the model. They are introduced in accordance with metaphysical principles [57]—principle (ii) in particular and they describe mutual interaction between the constituents. Their form (9) comes out as a consequence of the entropy inequality, and it secures the proper sign of the entropy production rate. In that sense, although the constituents are ideal gases, dissipation is brought into system through the source terms. Finally, their structure reveals also the physical basis of dissipation—it is a consequence of mutual exchange of momentum and energy between the constituents.

#### *2.2. Equations of State and Average Quantities for the Mixture*

The use of conservation laws (3) as governing equations implicitly introduces assumption that total pressure *p* of the mixture and intrinsic specific internal energy *ε<sup>I</sup>* mimic the equations of state of constituents. Namely:

$$p = \rho \frac{k\_{\rm B}}{m} T, \quad \varepsilon\_{I} = \frac{k\_{\rm B}}{(\gamma - 1)m} T. \tag{11}$$

where *T* is the average temperature of the mixture, *m* is the average atomic mass, and *γ* is the average ratio of specific heats. Total pressure and intrinsic specific internal energy must obey Dalton's law and defining relation (4):

$$p = \sum\_{a=1}^{n} p\_{a\prime} \quad \rho \varepsilon\_I = \sum\_{a=1}^{n} \rho\_a \varepsilon\_a. \tag{12}$$

In the sequel, we shall introduce the mass concentration variables *cα* along with restriction which they obey due to (4):

$$c\_{\mathfrak{a}} = \frac{\rho\_{\mathfrak{a}}}{\rho}, \quad \sum\_{\mathfrak{a}=1}^{n} c\_{\mathfrak{a}} = 1. \tag{13}$$

Since relations (11) and (12) must hold in local equilibrium, *<sup>T</sup>*<sup>1</sup> = ... = *Tn* = *<sup>T</sup>*, the following implicit definitions emerge:

$$\frac{1}{m} = \sum\_{a=1}^{n} \frac{c\_a}{m\_a}, \quad \frac{1}{(\gamma - 1)m} = \sum\_{a=1}^{n} \frac{c\_a}{(\gamma\_a - 1)m\_a}, \tag{14}$$

where, obviously, *m* = *m*(*cα*) and *γ* = *γ*(*cα*). On the other hand, in general case, from (11), (12) and (14) one obtains definition of the average temperature *T*:

$$T = (\gamma - 1)m \sum\_{a=1}^{n} \frac{\mathfrak{c}\_a}{(\gamma\_a - 1)m\_a} T\_a. \tag{15}$$

Note that although mixture equations of state (11) resemble the form of the same equations in the single-component gas, they depend on concentrations of the constituents due to (14). Furthermore, the mixture cannot be treated as ideal gas since the stress tensor **t**, apart from total pressure *p*, has non-vanishing diagonal and off-diagonal terms, and the flux of internal energy **q** does not vanish (see Equation (4)).

#### *2.3. Diffusion Fluxes and Diffusion Temperatures*

It is common in the mixture theory to introduce quantities which describe the difference between the state of constituents and certain properly defined state of the mixture. One such example is the diffusion velocity **<sup>u</sup>***<sup>α</sup>* <sup>=</sup> **<sup>v</sup>***<sup>α</sup>* <sup>−</sup> **<sup>v</sup>**. In this study we shall choose the so-called diffusion quantities which frequently appear in the literature. These are the diffusion fluxes **J***<sup>α</sup>* and the diffusion temperatures Θ*α*:

$$\mathbf{J}\_{\mathfrak{a}} = \rho\_{\mathfrak{a}} \mathbf{u}\_{\mathfrak{a}\mathfrak{a}} \quad \ominus \mathfrak{a}\_{\mathfrak{a}} = T\_{\mathfrak{a}} - T. \tag{16}$$

Due to (4), (14) and (15), they are not independent and must satisfy the following constraints: *<sup>n</sup>*

$$\sum\_{a=1}^{n} \mathbf{J}\_{a} = \mathbf{0}, \quad \sum\_{a=1}^{n} \frac{\mathbf{c}\_{a}}{(\gamma\_{a} - 1)m\_{a}} \boldsymbol{\Theta}\_{a} = \mathbf{0}. \tag{17}$$

Therefore, in the sequel whenever appears summation over *α*, or quantities related to *n*th constituent appear, the following relations will be assumed:

$$\mathbf{c}\_{\text{ll}} = 1 - \sum\_{b=1}^{n-1} \mathbf{c}\_{b\text{\prime}}; \quad \mathbf{J}\_{\text{ll}} = -\sum\_{b=1}^{n-1} \mathbf{J}\_{b\text{\prime}}; \quad \Theta\_{\text{ll}} = -\sum\_{b=1}^{n-1} \mathbf{Y}\_{b\text{\prime}} \Theta\_{b\text{\prime}} \tag{18}$$

where:

$$\mathcal{Y}\_{b\mathbb{N}} = \frac{c\_b}{c\_n} \frac{(\gamma\_n - 1)m\_n}{(\gamma\_b - 1)m\_b}. \tag{19}$$

Additionally, in definitions of all the mixture variables (4) partial mass densities *ρα* will be replaced with mass concentrations, i.e., *ρα* = *ρcα*, and diffusion velocities **u***<sup>α</sup>* will be replaced with diffusion fluxes, i.e., **u***<sup>α</sup>* = **J***α*/*ρα*.

#### *2.4. Phenomenological Coefficients*

As a continuum theory, RET is limited in the sense that phenomenological coefficients, which appear in the source terms (9), cannot be determined explicitly within its framework. This could be overcome by matching with an experimental evidence, or by means of some more refined approach. In this study, the latter path is taken, so that the kinetic theory of gases and moment equations appear as a natural choice. This procedure is completely developed in [33], albeit with slightly different notation. In the sequel, we shall briefly outline these results.

Basis for computation of the phenomenological coefficients, proposed in [33], is the multi-velocity and multi-temperature model of mixture of gases described by means of Boltzmann equations. The model is analyzed at an Euler level, i.e., using Maxwellian velocity distribution function for each constituent, with its own velocity and its own temperature. Moment equations consisted of mass, momentum and energy balance laws. However, source terms obtained as moments of collision operator had a structure that cannot be correlated to (9). Therefore, both source terms were linearized in the neighborhood of the local equilibrium state, and the explicit form of phenomenological coefficients *ψbc* and *θbc* was determined.

The source terms (9) can be linearized near the local equilibrium state determined with conditions (10). Firstly, the objective quantity **<sup>w</sup>** = **<sup>w</sup>**(**u1**, ... , **un**, *<sup>T</sup>*1, ... , *Tn*) evaluated at the local equilibrium state (**0**,..., **0**, *T*,..., *T*) is denoted with:

$$\mathbf{w}^0 = \mathbf{w}(\mathbf{0}, \dots, \mathbf{0}, T, \dots, T). \tag{20}$$

Then the source terms (9) can be approximated by:

$$\mathbf{f}\mathbf{\hat{n}}\_b \approx -\sum\_{c=1}^{n-1} \psi\_{bc}(\mathbf{w}^0) \left(\frac{\mathbf{u}\_c - \mathbf{u}\_n}{T}\right); \quad \mathbb{A}\_b \approx -\sum\_{c=1}^{n-1} \theta\_{bc}(\mathbf{w}^0) \left(\frac{T\_c - T\_n}{T^2}\right),$$

Such an approximation enables comparison with the source terms computed from the Boltzmann collision operator in a suitable asymptotics. Namely, in [33] the phenomenological coefficients are explicitly computed, depending on the parameter *sαβ* of the cross-section, and satisfying *<sup>s</sup>αβ* <sup>=</sup> *<sup>s</sup>βα* and *<sup>s</sup>αβ* <sup>&</sup>gt; <sup>−</sup>3/4 for each *<sup>α</sup>*, *<sup>β</sup>* <sup>=</sup> 1, ... , *<sup>n</sup>*. In our notation, they read:

$$\psi\_{\rm bc}(\mathbf{w}^{0}) = \begin{cases} -\frac{2}{3}T \,\mu\_{\rm bc} K\_{\rm bc'}^{\psi} & b \neq c, \\ \frac{2}{3}T \sum\_{\substack{\ell=1\\\ell \neq b}}^{n} \mu\_{\rm b\ell} K\_{\rm b\ell'}^{\psi} & b = c, \end{cases} \qquad \theta\_{\rm bc}(\mathbf{w}^{0}) = \begin{cases} -k\_{\rm B}T^{2} K\_{\rm bc'}^{\theta} & b \neq c, \\ k\_{\rm B}T^{2} \sum\_{\substack{\ell=1\\\ell \neq b}}^{n} K\_{\rm b\ell'}^{\theta} & b = c, \end{cases} \tag{21}$$

for any *<sup>b</sup>*, *<sup>c</sup>* <sup>=</sup> 1, ... , *<sup>n</sup>* <sup>−</sup> 1. Coefficients *<sup>K</sup><sup>ψ</sup> αβ* and *<sup>K</sup><sup>θ</sup> αβ* are defined for any *<sup>α</sup>*, *<sup>β</sup>* <sup>=</sup> 1, ... , *<sup>n</sup>*, as long as *<sup>α</sup>* <sup>=</sup> *<sup>β</sup>*, and read as follows:

$$K\_{a\beta}^{\Psi} = \mathbb{K} \, \kappa\_{a\beta}^{\Psi} \frac{\rho\_a}{m\_a} \frac{\rho\_\beta}{m\_\beta} \mu\_{a\beta}^{-\frac{s\_{a\beta}}{2}} (k\_\text{B}T)^{-\left(d\_a + d\_\beta\right) + \frac{s\_{a\beta}}{2}}, \quad K\_{a\beta}^{\theta} = \kappa\_{a\beta}^{\theta} K\_{a\beta}^{\psi} \tag{22}$$

where the following constants are introduced:

$$\begin{aligned} \kappa\_{a\beta}^{\Psi} &= \frac{2^{\frac{s\_{a\beta}}{2} + 4} \Gamma\left[\frac{s\_{a\beta} + 3}{2}\right]}{s\_{a\beta} + 5} \frac{\sqrt{\pi t}}{\Gamma[d\_{\alpha} + 1] \Gamma\left[d\_{\beta} + 1\right]},\\ \kappa\_{a\beta}^{\theta} &= \frac{2\left(s\_{a\beta} + 5\right)}{s\_{a\beta} + 7} \left(\frac{1}{s\_{a\beta} + 3} + \frac{\mu\_{a\beta}}{m\_{a} + m\_{\beta}}\right). \end{aligned} \tag{2.3}$$

In (22) and (23), *μαβ* and *d<sup>α</sup>* are the reduced mass and the parameter related to the ratio of specific heats, respectively:

$$
\mu\_{\alpha\beta} = \frac{m\_{\alpha}m\_{\beta}}{m\_{\alpha} + m\_{\beta}}, \quad d\_{\alpha} = \frac{-5\gamma\_{\alpha} + 7}{2(\gamma\_{\alpha} - 1)}.
$$

Dimensional constant K in (22) secures the proper dimension of the coefficients, and can be represented as a product of constituent-related part <sup>K</sup>*αβ* and common term <sup>K</sup><sup>∗</sup> in the following way:

$$\mathbb{K} = \mathbb{K}\_{\bullet} \, \mathbb{K}\_{a\beta\prime} \quad \mathbb{K}\_{a\beta} = \frac{m\_0^2}{\rho\_0^2} (k\_\mathrm{B} T\_0)^{d\_a + d\_\beta} \left(\frac{k\_\mathrm{B}}{\mu\_{a\beta}} T\_0\right)^{-\frac{s\_{a\beta}}{2}} \,\mathrm{}\,\tag{24}$$

where subscript 0 indicates the values of average quantities (mass, density and temperature) in a reference state.

#### *2.5. Relaxation Times*

Source terms (9) describe mutual interaction of the constituents and model dissipation in the system through relaxation towards (local) equilibrium. This process is usually expressed in terms of relaxation times—quantities which provide estimates of time-rate of convergence towards equilibrium. They can be related to phenomenological coefficients, and in the case of *n*-component mixture they read:

$$\tau\_{\rm Dbc} = (-1)^{1-\delta\_{\rm lc}} \frac{1}{\psi\_{\rm lc}(\mathbf{w}^0)} \frac{\rho\_{\rm c} \rho\_{\rm n}}{\sum\_{a=1}^n \rho\_a} T, \quad \tau\_{\rm Tbc} = (-1)^{1-\delta\_{\rm bc}} \frac{1}{\theta\_{\rm bc}(\mathbf{w}^0)} \frac{\rho\_{\rm c} \upsilon\_{\rm c} \rho\_{\rm n} \upsilon\_{\rm Vn}}{\sum\_{a=1}^n \rho\_a \upsilon\_{\rm Va}} T^2,\tag{25}$$

where *τDbc* and *τTbc* denote mechanical (diffusion) and thermal relaxation time, respectively, and factor (−1)1−*δbc* , with *<sup>δ</sup>bc* being Kronecker delta, is used to secure their positivity. Using relaxation times (25), source terms may be formally expressed in the following form:

$$\begin{split} \dot{\mathbf{m}}\_{b} &= -\sum\_{c=1}^{n-1} (-1)^{1-\delta\_{\rm lc}} \frac{1}{\tau\_{\rm Dlc}(\mathbf{w}^{0})} \frac{\rho\_{\rm c} \rho\_{\rm u}}{\sum\_{a=1}^{n} \rho\_{a}} T \left( \frac{\mathbf{u}\_{\rm c}}{T\_{\rm c}} - \frac{\mathbf{u}\_{\rm n}}{T\_{n}} \right); \\ \dot{\boldsymbol{\varepsilon}}\_{b} &= -\sum\_{c=1}^{n-1} (-1)^{1-\delta\_{\rm lc}} \frac{1}{\tau\_{\rm Dlc}(\mathbf{w}^{0})} \frac{\rho\_{\rm c} \mathbf{c}\_{\rm Vc} \rho\_{\rm n} \mathbf{c}\_{\rm Vn}}{\sum\_{a=1}^{n} \rho\_{a} \mathbf{c}\_{\rm Va}} T^{2} \left( -\frac{1}{T\_{\rm c}} + \frac{1}{T\_{\rm n}} \right) . \end{split} \tag{26}$$

It may be observed that relaxation times (25) present proper generalization of the relaxation times given in [47].

#### **3. Shock Structure Problem for Three-Component Mixture of Euler Fluids**

Our main concern in this study is the shock structure problem for three-component mixture of Euler fluids. In idealized situation (non-dissipative models) shock waves are (singular) surfaces, which move through space, on which jumps of field variables occur. If *s* denotes the shock speed—propagation velocity component perpendicular to shock surface—then the values of field variables in front and behind the shock are related to *s* by celebrated Rankine–Hugoniot conditions. In our analysis we shall analyze single shock whose surface is plane, and the mixture flows in direction perpendicular to the shock wave. If the dissipation is present in the model, shock wave is smeared out into the shock structure (profile). This is smooth solution, in the form of traveling wave, which has steep gradients in the small neighborhood of the shock wave and asymptotically connects equilibrium states in front (upstream) and behind (downstream) the shock. Equilibrium states are related to the shock speed *s* through the Rankine–Hugoniot equations.

With these preliminaries at hand, we may sketch the steps which will lead us to a formulation of the shock structure problem for three-component mixture of Euler fluids (*n* = 3). Since the analysis is restricted to plane waves, we shall first expose governing equations for one-dimensional flow. Then, taking the traveling wave solution ansatz we shall derive equations which determine the shock structure. Finally, the shock structure problem will be put into dimensionless form which will be subsequently numerically integrated and analyzed.

#### *3.1. Equations for One-Dimensional Flow*

In one dimensional settings, certain simplifications have to be done. If **<sup>e</sup>***i*, *<sup>i</sup>* = 1, 2, 3, represents standard basis in R3, then we assume motion in **e**<sup>1</sup> direction and respective coordinate denote by *<sup>x</sup>*. Velocities will be **<sup>v</sup>** = *<sup>v</sup>***e**<sup>1</sup> and **<sup>u</sup>***<sup>α</sup>* = *<sup>u</sup>α***e**1, and diffusion fluxes **<sup>J</sup>***<sup>α</sup>* <sup>=</sup> *<sup>J</sup>α***e**1. Mixture stress tensor will have the form **<sup>t</sup>** <sup>=</sup> *<sup>t</sup>*11**e**<sup>1</sup> <sup>⊗</sup> **<sup>e</sup>**1, and flux of internal energy will be **<sup>q</sup>** = *<sup>q</sup>***e**1, where:

$$t\_{11} = -p - \sum\_{a=1}^{3} \frac{1}{\rho\_a} l\_{a\prime}^2 \quad q = \sum\_{a=1}^{3} \left( \varepsilon\_a + \frac{p\_a}{\rho\_a} + \frac{1}{2\rho\_a^2} l\_a^2 \right) l\_{a\prime} \tag{27}$$

whose more explicit form will be given in the sequel. Source terms then read **<sup>m</sup>**<sup>ˆ</sup> *<sup>b</sup>* = *<sup>m</sup>*<sup>ˆ</sup> *<sup>b</sup>***e**1.

Simplifications introduced above lead us to the following one dimensional governing equations, conservation laws for the mixture:

$$\begin{aligned} \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x} (\rho v) &= 0, \\ \frac{\partial}{\partial t} (\rho v) + \frac{\partial}{\partial x} \left(\rho v^2 - t\_{11}\right) &= 0, \\ \frac{\partial}{\partial t} \left(\frac{1}{2} \rho v^2 + \rho \varepsilon\right) + \frac{\partial}{\partial x} \left[\left(\frac{1}{2} \rho v^2 + \rho \varepsilon\right) v - t\_{11} v + q\right] &= 0, \end{aligned} \tag{28}$$

and balance laws for the constituents:

$$\begin{split} \frac{\partial}{\partial t}(\rho c\_b) + \frac{\partial}{\partial x}(\rho c\_b v + f\_b) &= 0, \\ \frac{\partial}{\partial t}(\rho c\_b v + f\_b) + \frac{\partial}{\partial x} \left(\rho c\_b v^2 + \frac{f\_b^2}{\rho c\_b} + 2v f\_b + p\_b \right) &= \dot{m}\_{b\prime} \\ \frac{\partial}{\partial t} \left[ \frac{1}{2} \rho c\_b \left(v + \frac{f\_b}{\rho c\_b} \right)^2 + \rho c\_b \varepsilon\_b \right] \\ &+ \frac{\partial}{\partial x} \left\{ \left[ \frac{1}{2} \rho c\_b \left(v + \frac{f\_b}{\rho c\_b} \right)^2 + \rho c\_b \varepsilon\_b + p\_b \right] \left(v + \frac{f\_b}{\rho c\_b} \right) \right\} = \dot{\varepsilon}\_b + \dot{m}\_b v, \end{split} \tag{29}$$

where *b* = 1, 2.

#### *3.2. Shock Structure Equations*

To derive equations which determine the shock structure we assume the traveling wave profile of the solution. In that case generic field variable *U*(*t*, *x*) depends on a single composite variable *<sup>ξ</sup>* <sup>=</sup> *<sup>x</sup>* <sup>−</sup> *st*, *<sup>U</sup>* <sup>=</sup> *<sup>U</sup>*(*ξ*), where *<sup>s</sup>* is the shock speed. Partial derivatives are replaced with ordinary derivatives in the following manner, *<sup>∂</sup>*/*∂<sup>t</sup>* <sup>=</sup> <sup>−</sup>*s*(d/d*ξ*), *<sup>∂</sup>*/*∂<sup>x</sup>* <sup>=</sup> d/d*ξ*, and relative velocity with respect to the shock front is defined as *<sup>u</sup>* <sup>=</sup> *<sup>v</sup>* <sup>−</sup> *<sup>s</sup>*. After straightforward computations our system of governing equations is transformed into ODE system of shock structure equations, i.e., conservation laws for the mixture:

$$\begin{aligned} \frac{\mathbf{d}}{\mathbf{d}\_{\sharp}^{\mathcal{L}}}(\rho u) &= 0, \\ \frac{\mathbf{d}}{\mathbf{d}\_{\sharp}^{\mathcal{L}}} \Big(\rho u^2 - t\_{11}\Big) &= 0, \\ \frac{\mathbf{d}}{\mathbf{d}\_{\sharp}^{\mathcal{L}}} \Big[ \Big(\frac{1}{2}\rho u^2 + \rho \varepsilon\Big)u - t\_{11}u + q\Big] &= 0, \end{aligned} \tag{30}$$

and balance laws for the constituents:

$$\begin{aligned} \frac{\mathbf{d}}{\mathbf{d}\_b^x} (\rho c\_b u + l\_b) &= 0, \\ \frac{\mathbf{d}}{\mathbf{d}\_b^x} \left(\rho c\_b u^2 + \frac{l\_b^2}{\rho c\_b} + 2uI\_b + p\_b\right) &= \hat{m}\_{b\prime} \\ \frac{\mathbf{d}}{\mathbf{d}\_b^x} \left\{ \left[\frac{1}{2}\rho c\_b \left(u + \frac{l\_b}{\rho c\_b}\right)^2 + \rho c\_b \varepsilon\_b + p\_b \right] \left(u + \frac{l\_b}{\rho c\_b}\right) \right\} &= \pounds\_b + \hat{m}\_b u, \end{aligned} \tag{31}$$

where *b* = 1, 2. It may be noted that equations (30) express the conservation of fluxes of mass, momentum and energy of the mixture across the shock wave, while (31)1 expresses the conservation of mass fluxes of the constituents. On the other hand, equations (31)2,3 imply that fluxes of momenta and energies of the constituents are not conserved, but rather balanced by the mutual interactions with other constituents (source terms).

#### *3.3. Dimensionless Shock Structure Equations*

To obtain results as general as possible, we transform our problem into dimensionless form. Field variables will be scaled with respect reference values, which are taken to be ones in upstream equilibrium indicated by the subscript 0. To that end we shall introduce reference (equilibrium) mass density *<sup>ρ</sup>*0, temperature *<sup>T</sup>*0, and velocity *<sup>a</sup>*<sup>0</sup> = (*γ*0(*k*B/*m*0)*T*0)1/2, which is the speed of sound in upstream equilibrium. Here, *m*<sup>0</sup> and *γ*<sup>0</sup> are average mass and average ratio of specific heats in upstream equilibrium. Reference length *l*<sup>0</sup> is expressed in terms of *a*<sup>0</sup> and diffusion relaxation time *τ*<sup>0</sup> *<sup>D</sup>*<sup>12</sup> evaluated in upstream equilibrium, *<sup>l</sup>*<sup>0</sup> = *<sup>a</sup>*0*τ*<sup>0</sup> *<sup>D</sup>*12. Dimensionless variables now read:

$$\tilde{\rho} = \frac{\rho}{\rho\_0}, \quad \tilde{u} = \frac{u}{a\_0}, \quad \tilde{T} = \frac{T}{T\_0}, \quad \tilde{l}\_b = \frac{f\_b}{\rho\_0 a\_0}, \quad \tilde{\Theta}\_b = \frac{\Theta\_b}{T\_0}, \quad \tilde{\xi} = \frac{\tilde{\xi}}{a\_0 \tau\_{D12}^0}.\tag{32}$$

Without ambiguity we may drop tilde in subsequent equations. Mixture conservation laws for the shock structure read:

$$\begin{aligned} \frac{\mathbf{d}}{\mathbf{d}\xi}(\rho u) &= 0, \\ \frac{\mathbf{d}}{\mathbf{d}\xi} \left(\rho u^2 - t\_{11}\right) &= 0, \\ \frac{\mathbf{d}}{\mathbf{d}\xi} \left[\left(\frac{1}{2}\rho u^2 + \rho \varepsilon\right)u - t\_{11}u + q\right] &= 0, \end{aligned} \tag{33}$$

where we the following relations hold:

$$\begin{split} t\_{11} &= -\frac{1}{\gamma\_{0}} \frac{m\_{0}}{m} \rho T - \frac{1}{\rho c\_{1}} l\_{1}^{2} - \frac{1}{\rho c\_{2}} l\_{2}^{2} - \frac{1}{\rho c\_{3}} l\_{3}^{2}, \\ \rho \varepsilon &= \frac{m\_{0}}{\gamma\_{0}} \frac{1}{(\gamma - 1)m} \rho T + \frac{1}{2\rho c\_{1}} l\_{1}^{2} + \frac{1}{2\rho c\_{2}} l\_{2}^{2} + \frac{1}{2\rho c\_{3}} l\_{3}^{2}, \\ q &= \left[ \frac{1}{2\rho^{2} c\_{1}^{2}} l\_{1}^{2} + \frac{1}{\gamma\_{0}} \frac{m\_{0}}{m\_{1}} \frac{\gamma\_{1}}{\gamma\_{1} - 1} (T + \Theta\_{1}) \right] l\_{1} \\ &+ \left[ \frac{1}{2\rho^{2} c\_{2}^{2}} l\_{2}^{2} + \frac{1}{\gamma\_{0}} \frac{m\_{0}}{m\_{2}} \frac{\gamma\_{2}}{\gamma\_{2} - 1} (T + \Theta\_{2}) \right] l\_{2} \\ &+ \left[ \frac{1}{2\rho^{2} c\_{3}^{2}} l\_{3}^{2} + \frac{1}{\gamma\_{0}} \frac{m\_{0}}{m\_{3}} \frac{\gamma\_{3}}{\gamma\_{3} - 1} (T + \Theta\_{3}) \right] l\_{3}. \end{split} \tag{34}$$

In Equation (34) we have to take into account constraints (17):

$$c\_3 = 1 - c\_1 - c\_2, \quad l\_3 = -l\_1 - l\_2, \quad \Theta\_3 = -\Upsilon\_{13}\Theta\_1 - \Upsilon\_{23}\Theta\_2. \tag{35}$$

where *Υ*<sup>13</sup> and *Υ*<sup>23</sup> are determined by (19). Balance laws for the constituents read:

$$\begin{cases} \frac{\mathbf{d}}{\mathbf{d}\_b^T} (\rho c\_b \mu + f\_b) = 0, \\ \frac{\mathbf{d}}{\mathbf{d}\_b^T} \left( \rho c\_b \mu^2 + \frac{f\_b^2}{\rho c\_b} + 2u f\_b + \frac{1}{\gamma\_0} \frac{m\_0}{m\_b} \rho c\_b (T + \Theta\_b) \right) = \hat{m}\_b, \\ \frac{\mathbf{d}}{\mathbf{d}\_b^T} \left\{ \left[ \frac{1}{2} \rho c\_b \left( \mu + \frac{f\_b}{\rho c\_b} \right)^2 + \frac{m\_0}{m\_b} \frac{\gamma\_b}{\gamma\_0} \frac{1}{\gamma\_b - 1} \rho c\_b (T + \Theta\_b) \right] \left( \mu + \frac{f\_b}{\rho c\_b} \right) \right\} = \hat{\varepsilon}\_b + \hat{m}\_b \mu \end{cases} (36)$$

where *<sup>b</sup>* = 1, 2. Dimensionless source terms *<sup>m</sup>*<sup>ˆ</sup> *<sup>b</sup>* read:

$$\begin{split} \hat{m}\_{1} &= -c\_{20}c\_{30} \left[ \frac{\psi\_{11}}{\psi\_{12}^{0}} \left( \frac{J\_{1}}{\rho c\_{1}(T + \Theta\_{1})} - \frac{f\_{3}}{\rho c\_{3}(T + \Theta\_{3})} \right) \right. \\ &\left. + \frac{\psi\_{12}}{\psi\_{12}^{0}} \left( \frac{f\_{2}}{\rho c\_{2}(T + \Theta\_{2})} - \frac{f\_{3}}{\rho c\_{3}(T + \Theta\_{3})} \right) \right], \\ \hat{m}\_{2} &= -c\_{20}c\_{30} \left[ \frac{\psi\_{21}}{\psi\_{12}^{0}} \left( \frac{J\_{1}}{\rho c\_{1}(T + \Theta\_{1})} - \frac{f\_{3}}{\rho c\_{3}(T + \Theta\_{3})} \right) \right. \\ &\left. + \frac{\psi\_{22}}{\psi\_{12}^{0}} \left( \frac{J\_{2}}{\rho c\_{2}(T + \Theta\_{2})} - \frac{f\_{3}}{\rho c\_{3}(T + \Theta\_{3})} \right) \right]. \end{split} (37)$$

Dimensionless source terms *e*ˆ*<sup>b</sup>* read:

$$\begin{split} \mathfrak{E}\_{1} &= -\mathfrak{c}\_{20}\mathfrak{c}\_{30} \left[ \frac{\theta\_{11}}{a\_{0}^{2}\mathfrak{p}\_{12}^{0}} \frac{\Theta\_{1} - \Theta\_{3}}{(T + \Theta\_{1})(T + \Theta\_{3})} + \frac{\theta\_{12}}{a\_{0}^{2}\mathfrak{p}\_{12}^{0}} \frac{\Theta\_{2} - \Theta\_{3}}{(T + \Theta\_{2})(T + \Theta\_{3})} \right], \\ \mathfrak{E}\_{2} &= -\mathfrak{c}\_{20}\mathfrak{c}\_{30} \left[ \frac{\theta\_{21}}{a\_{0}^{2}\mathfrak{p}\_{12}^{0}} \frac{\Theta\_{1} - \Theta\_{3}}{(T + \Theta\_{1})(T + \Theta\_{3})} + \frac{\theta\_{22}}{a\_{0}^{2}\mathfrak{p}\_{12}^{0}} \frac{\Theta\_{2} - \Theta\_{3}}{(T + \Theta\_{2})(T + \Theta\_{3})} \right]. \end{split} \tag{38}$$

In (37) and (38), *c*<sup>20</sup> and *c*<sup>30</sup> denote mass concentrations in upstream equilibrium. Ratios of phenomenological coefficients are given in the Appendix A.

Governing Equations (33) and (36), along with (34), (37) and (38), constitute the system of (dimensionless) shock structure equations for a multi-temperature three-component mixture of Euler fluids. If the shock structure is analyzed in a single-temperature model, then equations (36)3 should be dropped from the system, and <sup>Θ</sup><sup>1</sup> <sup>=</sup> <sup>Θ</sup><sup>2</sup> <sup>=</sup> <sup>Θ</sup><sup>3</sup> <sup>≡</sup> 0 set in the remaining equations.

#### *3.4. Parameters of the Model*

The dimensionless form of the model substantially reduces the number of parameters. In previous study of binary mixtures [47] there were three parameters upon which the solution depend—Mach number and mass concentration of one constituent in upstream equilibrium, and mass ratio of the constituents. Increase of the number of constituents increases the number parameters. Thus, in this study, in which we deal with *n* = 3, constituents the following parameters will be used:

$$\begin{aligned} M\_0 &= \frac{u\_0}{a\_0}, & \text{Mach number in upstream equilibrium;}\\ c\_{10} &= \left(\frac{\rho\_1}{\rho}\right)\_0, c\_{20} = \left(\frac{\rho\_2}{\rho}\right)\_0, & \text{mass concentrations in upstream equilibrium;}\\ \mu\_1 &= \frac{m\_1}{m\_3}, \; \mu\_2 = \frac{m\_2}{m\_3}, & \text{mass ratios of the constituents.} \end{aligned}$$

To avoid ambiguity, we shall assume that the following inequality holds, *m*<sup>1</sup> < *m*<sup>2</sup> < *m*3, which implies 0 < *μ*<sup>1</sup> < *μ*<sup>2</sup> < 1. Additionally, the ratios of specific heats of the constituents *γα*, *α* = 1, 2, 3, have to be fixed.

#### *3.5. Numerical Procedure*

Governing equations of the shock structure consist of conservation laws (33) for the mixture, and balance laws for the constituents (36), where explicit form of the fluxes and source terms are given by (34), (37) and (38), respectively. Governing equations could be formally written in the form:

$$\frac{\mathbf{d}}{\mathbf{d}\tilde{\xi}}\mathbf{F}(\mathbf{U}(\xi)) = \mathbf{f}(\mathbf{U}(\xi)),\tag{39}$$

where **<sup>U</sup>** = (*ρ*, *<sup>u</sup>*, *<sup>T</sup>*, *<sup>c</sup>*1, *<sup>c</sup>*2, *<sup>J</sup>*1, *<sup>J</sup>*2, <sup>Θ</sup>1, <sup>Θ</sup>2)*<sup>T</sup>* denotes the vector of field variables, and **<sup>F</sup>**(**U**) and **f**(**U**) have obvious meaning.

Shock structure is continuous solution of the system (39) which asymptotically connects upstream equilibrium **U**<sup>−</sup> with downstream equilibrium **U**+, i.e., a heteroclinic orbit. Equilibrium states are not independent, but related by Rankine–Hugoniot equations. As a consequence, the system of first-order ordinary differential equations (39) is adjoined with the following boundary conditions:

$$\begin{aligned} \mathbf{U}(-\infty) &= \lim\_{\frac{\mathbf{r}}{\xi} \to -\infty} \mathbf{U}(\xi) = \mathbf{U}\_{-\prime} \quad \mathbf{U}(+\infty) = \lim\_{\frac{\mathbf{r}}{\xi} \to +\infty} \mathbf{U}(\xi) = \mathbf{U}\_{+\prime} \\ &\mathbf{U}^{\prime}(\pm\infty) = \lim\_{\frac{\mathbf{r}}{\xi} \to \pm\infty} \frac{\mathbf{d} \mathbf{U}(\frac{\mathbf{r}}{\xi})}{\mathbf{d}\_{\natural}^{\mathbb{Z}}} = \mathbf{0}, \end{aligned} \tag{40}$$

where the boundary data are as follows:

$$\mathbf{U}\_{-}=\begin{bmatrix}\rho\_{-}\\u\_{-}\\T\_{-}\\c\_{1}&\\c\_{2}&\\&c\_{2}\\
\end{bmatrix}=\begin{bmatrix}1\\M\_{0}\\1\\c\_{10}\\c\_{20}\\0\\0\\0\\0\\0\end{bmatrix},\quad\mathbf{U}\_{+}=\begin{bmatrix}\rho\_{+}\\u\_{+}\\T\_{+}\\c\_{1+}\\c\_{2+}\\0\\0\\0\end{bmatrix}=\begin{bmatrix}\frac{M\_{0}^{2}}{1-\rho\_{0}^{2}+\mu\_{0}^{2}M\_{0}^{2}}\\\-\frac{1}{\rho\_{0}^{2}+\mu\_{0}^{2}M\_{0}^{2}}\\\-\frac{1}{\rho\_{0}}+\frac{\mu\_{0}^{2}+\mu\_{0}^{2}M\_{0}^{2}}{M\_{0}}\\\-\frac{1}{\rho\_{0}}+\frac{\mu\_{0}^{2}+\mu\_{0}^{2}M\_{0}^{2}}{M\_{0}}\\\-\frac{1}{\rho\_{0}}\\0\\\-\frac{\mu\_{0}}{\rho\_{2}}\\0\\0\end{bmatrix},\tag{41}$$

where *μ*<sup>2</sup> <sup>0</sup> = (*γ*<sup>0</sup> <sup>−</sup> <sup>1</sup>)/(*γ*<sup>0</sup> <sup>+</sup> <sup>1</sup>). In the case of a ST mixture, field variables <sup>Θ</sup><sup>1</sup> and <sup>Θ</sup><sup>2</sup> should be dropped from the vector **U**, as well as from boundary data (41).

Note that the role of Rankine–Hugoniot equations in this problem is not the same as in the models without dissipation, like Euler equations for single-component gases [1]. In the Euler equations, they relate the jump of field variables at the shock wave (singular surface) to the shock speed. Here, instead, they relate the equilibrium states, upstream and downstream, which are connected by a continuous shock structure. However, they coincide with Rankine–Hugoniot equations for Euler equations because equilibrium states lie on equilibrium manifold, and the equilibrium subsystem for our model coincides with the Euler equations. They are derived by means of integration of conservation laws (33), taking into account equilibrium values of non-convective fluxes *<sup>t</sup>*<sup>11</sup> and *<sup>q</sup>* from (34) and *<sup>J</sup><sup>α</sup>* = 0.

Numerical computation of the shock structure in hyperbolic systems of balance laws is delicate because continuous solution does not exist for all values of the shock speed. Namely, it was shown that continuous solution breaks down, and there appears the profile with a sub-shock, if shock speed exceeds the greatest characteristic speed of the system in upstream equilibrium [38]. However, there may also appear a regular singularity, within the continuous profile, if the shock speed reaches the critical value in downstream equilibrium. These cases are carefully analyzed in [35], and existence of continuous shock profiles in mixtures was the subject of several recent studies [48,49].

To find numerical solution of the problem (39)–(41), the infinite domain of definition of heteroclinic orbit, −∞ < *ξ* < ∞, has to be truncated; however, the new finite domain, *<sup>ξ</sup>* <sup>∈</sup> [*ξ*0, *<sup>ξ</sup>*1], has to be large enough to secure that the whole profile is properly captured. Further computation may be based upon two different strategies. First one takes into account the type of stationary/equilibrium points and finds the shock structure as solution of the initial value problem. This approach is described in details in [47]. It is simple, straightforward, with low computational cost. Nevertheless, its applicability is limited since regular singularity, if it appears, cannot be overcome. Second strategy uses the finite-difference scheme to compute the shock structure as solution of the boundary value problem. This method may require several steps of computation. Adaptive step size has to be used, with smaller steps in the sub-domain in which steep gradients of field variables

occur. However, position of this sub-domain is not known beforehand. Additionally, an initial guess has to be provided for numerical computation. On the other hand, boundary value approach may produce continuous solution even when regular singularity appear in the domain. Details about this approach may be found in [18,53].

Our analysis will try to take advantages of both strategies described above. In cases in which initial value problem can be solved, we shall first use this method, and then exploit this solution as initial guess for the boundary value problem. In such a way we shall get the estimate of domain width and the location of sub-domain with steep gradients in first place. After that, finite-difference scheme will be used to further improve the solution. On the other hand, when initial value strategy is not applicable, we shall use as an initial guess the continuous solution obtained for some other (not too distinct) values of parameters.

#### **4. Analysis of the Shock Structure**

Analysis of the shock structure in three-component mixture of Euler fluids has several goals. First, we want to validate our model, especially phenomenological coefficients, by computing typical profiles that are comparable with the ones obtained by other methods. Second, we want to recognize the features typical for polyatomic gases. Third, we want get a glimpse of the influence of Mach number on the shock structure. Fourth, we want set up the data which mimic the binary mixture and find out what are the differences that three-component model brings. Finally, we shall compare the multi-temperature shock structure with the single-temperature one, and try to estimate the effect obtained by the MT assumption.

Computation are performed in the way explained in Section 3. However, field variables are scaled to change in the [0, 1] interval, where it is appropriate.

#### *4.1. Shock Structure in the MT Model*

Our model contains several parameters, whose values have to be given prior to numerical computation. First, phenomenological coefficients (21) depend on the type of particle interaction through the parameter *sαβ* of the cross-section, (22) and (23). In the sequel, we shall assume *sαβ* = 1, which corresponds to the hard-sphere model of particle interaction. Furthermore, in all the computations we shall assume the following values of ratios of specific heats: *<sup>γ</sup>*<sup>1</sup> = 7/5, *<sup>γ</sup>*<sup>2</sup> = *<sup>γ</sup>*<sup>3</sup> = 9/7. In most of the computations mass ratios will be *<sup>μ</sup>*<sup>1</sup> = 0.01 and *<sup>μ</sup>*<sup>2</sup> = 0.1, unless stated otherwise (see Case 4). This assumption tends to cover the most interesting case from the MT point of view—the one with a mixture whose constituents have disparate masses.

Case 1. In this case it is analyzed the situation with small equilibrium concentration *c*<sup>30</sup> of the heaviest constituent, and large equilibrium concentration *c*20. Shock structure is computed for *<sup>c</sup>*<sup>10</sup> = 0.10, *<sup>c</sup>*<sup>20</sup> = 0.75 and *<sup>c</sup>*<sup>30</sup> = 0.15, and *<sup>M</sup>*<sup>0</sup> = 1.30, and shown in Figure 1. Average mixture variables *ρ*, *u* and *T* have proper monotonic profiles, as well as velocities of the constituents. However, it may be observed that temperature *T*<sup>3</sup> of the heaviest constituent has an overshoot—the zone in which the variable has the value greater than the downstream equilibrium value. This result corresponds to similar observations in binary mixtures, where an overshoot appears when heavier constituents have small concentration [41,42,47]. It may also be noted that the difference between the velocities of the constituents is significant. Therefore, this case corresponds to typical profiles which can be found in [55,56]. Implicitly, it validates the structure of phenomenological coefficients (21) obtained by matching the continuum and kinetic approach.

Case 2. In the second case, two profiles were computed for the same value of Mach number, *<sup>M</sup>*<sup>0</sup> = 1.20. The difference appears in equilibrium concentrations: for *<sup>c</sup>*<sup>10</sup> = 0.60, *<sup>c</sup>*<sup>20</sup> = 0.30, and *<sup>c</sup>*<sup>30</sup> = 0.10, profiles are presented in Figure 2, graphs (a) and (b); for *<sup>c</sup>*<sup>10</sup> = 0.60, *<sup>c</sup>*<sup>20</sup> = 0.20, and *<sup>c</sup>*<sup>30</sup> = 0.20, profiles are presented in Figure 2, graphs (c) and (d). The profiles consist of two zones: upstream one with steeper gradients of field variables, and downstream one with moderate gradients. This appears to be more prominent in the profiles of mixture variables, especially on graph (a). This kind of profile, named Type-C, was firstly observed in [59] for the shock structure in ET14 model of a polyatomic gas (see also [18] for a contextualized account). The result was further supported by numerical solution of the Boltzmann equation for polyatomic gases [60], as well as two-temperature Navier–Stokes equations for polyatomic gases [61]. The main cause of such a profile is attributed to the presence of dynamic pressure and its relaxation (bulk viscosity), which is longer than the relaxation of viscous stress and heat flux. Although our model does not take into account these processes, there exists certain resemblance to this phenomenon in our profiles as well. To the best of our knowledge, such profiles do not appear in monatomic gases and their mixtures, and we are convinced that polyatomic molecular structure is the cause for this result. Note that such profiles appear in the case with an overshoot and without it alike.

**Figure 1.** Shock structure for *<sup>c</sup>*<sup>10</sup> = 0.10, *<sup>c</sup>*<sup>20</sup> = 0.75, *<sup>c</sup>*<sup>30</sup> = 0.15, and *<sup>M</sup>*<sup>0</sup> = 1.30: (**a**) mixture field variables; (**b**) velocities; (**c**) temperatures; (**d**) diffusion fluxes; (**e**) concentrations.

**Figure 2.** Shock structure for *<sup>M</sup>*<sup>0</sup> = 1.20, *<sup>c</sup>*<sup>10</sup> = 0.60, *<sup>c</sup>*<sup>20</sup> = 0.30, and *<sup>c</sup>*<sup>30</sup> = 0.10: (**a**) mixture field variables; (**b**) temperatures. Shock structure for *<sup>M</sup>*<sup>0</sup> = 1.20, *<sup>c</sup>*<sup>10</sup> = 0.60, *<sup>c</sup>*<sup>20</sup> = 0.20, and *<sup>c</sup>*<sup>30</sup> = 0.20: (**c**) mixture field variables; (**d**) temperatures.

Case 3. In this case we examined the profiles with the highest equilibrium concentration *<sup>c</sup>*30. In particular, we computed the profiles for *<sup>c</sup>*<sup>10</sup> = 0.20, *<sup>c</sup>*<sup>20</sup> = 0.20 and *<sup>c</sup>*<sup>30</sup> = 0.60, and for two different values of Mach number, *<sup>M</sup>*<sup>0</sup> = 1.20 and *<sup>M</sup>*<sup>0</sup> = 1.40, presented in Figure 3a,b, respectively. For these profiles it is remarkable that temperatures *Tα* of the constituents have negligibly small difference, which is almost indistinguishable on the scale of graph (and thus not presented here). Therefore, we present the two cases which are in accordance with the results for binary mixture [47]—shock thickness decreases with the increase of Mach number. This conclusion is valid for small values of *M*0.

**Figure 3.** Shock structure for *<sup>c</sup>*<sup>10</sup> = 0.20, *<sup>c</sup>*<sup>20</sup> = 0.20 and *<sup>c</sup>*<sup>30</sup> = 0.60: (**a**) mixture field variables for *<sup>M</sup>*<sup>0</sup> = 1.20; (**b**) mixture field variables for *<sup>M</sup>*<sup>0</sup> = 1.40.

Case 4. Here we computed the profiles for *<sup>c</sup>*<sup>10</sup> = 0.40, *<sup>c</sup>*<sup>20</sup> = 0.50 and *<sup>c</sup>*<sup>30</sup> = 0.10. However, we analyzed two different cases: (a) *<sup>M</sup>*<sup>0</sup> = 1.20 and *<sup>μ</sup>*<sup>2</sup> = 0.10, and (b) *<sup>M</sup>*<sup>0</sup> = 1.40 and *<sup>μ</sup>*<sup>2</sup> = 0.80, presented in Figure 4. These results give an insight into the influence of mass ratio on the shock structure. In case (a), in the mixture with disparate masses of constituents, one may observe the temperature overshoot in *T*3. However, in case (b) where constituents 2 and 3 have comparable masses, temperatures of the constituents are almost indistinguishable.

**Figure 4.** Shock structure for *<sup>c</sup>*<sup>10</sup> = 0.40, *<sup>c</sup>*<sup>20</sup> = 0.50 and *<sup>c</sup>*<sup>30</sup> = 0.10: (**a**) temperatures for *<sup>M</sup>*<sup>0</sup> = 1.20 and *<sup>μ</sup>*<sup>2</sup> = 0.10; (**b**) temperatures for *<sup>M</sup>*<sup>0</sup> = 1.40 and *<sup>μ</sup>*<sup>2</sup> = 0.80.

#### *4.2. Comparison with a Single-Temperature Model*

There always appears a question of the relevance of the MT assumption in the mixtures. Certainly, there are non-equilibrium processes in which this assumption cannot be dropped (e.g., in plasma). On the other hand, it also proved to be relevant in mixtures whose constituents have disparate atomic/molecular masses, like noble gases. Within the framework of RET, even if the multi-temperature assumption is dropped, one ends up with a multi-velocity model which is still hyperbolic and does not produce the paradox of infinite speed of pulse propagation.

Our intention is to quantify the discrepancy in mixture field variables, primarily in the average temperature of the mixture, which appears within the shock profile computed for MT model, and the one computed for ST model. To that end, we compute the sup-norm of the difference of numerically computed solutions, i.e., *<sup>T</sup>* <sup>−</sup> *<sup>T</sup>*ST <sup>=</sup> sup*ξ*∈[*ξ*0,*ξ*1] <sup>|</sup>*T*(*ξ*) <sup>−</sup> *<sup>T</sup>*ST(*ξ*)|, where *<sup>T</sup>*ST(*ξ*) denotes the temperature field along the shock profile in ST model.

For the shock structure with a temperature overshoot, presented in Case 1. of this study, sup-norm of the difference of normalized temperature profiles (the ones whose equilibrium values are set to 0.0 and 1.0) is *T* − *T*ST ∼ 0.01 (see Figure 5). This implies that processes of energy exchange between the constituents, which appear due to the temperature difference, do not have a significant influence on the average value of temperature in MT model. However, this does not mean that they may be neglected. It was shown that appearance of temperature overshoot is directly related to insufficient energy exchange/transfer [43,47]. Additionally, preliminary computations, not reported here, show that there are cases with extremely large temperature overshoot, (*T*max <sup>3</sup> <sup>−</sup> *<sup>T</sup>*+)/(*T*<sup>+</sup> <sup>−</sup> *<sup>T</sup>*−) <sup>≈</sup> 0.70, in which ST is not capable of proper capturing the behavior of the mixture temperature. This will be the subject of our prospective study.

**Figure 5.** Difference of the mean temperatures in MT and ST model, computed for the Case 1, *<sup>c</sup>*<sup>10</sup> = 0.10, *<sup>c</sup>*<sup>20</sup> = 0.75, *<sup>c</sup>*<sup>30</sup> = 0.15, and *<sup>M</sup>*<sup>0</sup> = 1.30.

#### **5. Relaxation in the Shock Structure**

One of the main features of hyperbolic systems of balance laws is that dissipation is achieved through the relaxation process, i.e., convergence of non-equilibrium variables to a certain (local) equilibrium state. Quantitative measure of the rate of convergence towards equilibrium is relaxation time. In simpler mathematical models, relaxation time may be unique for the whole system. However, if the model describes complex physical processes, there may be a need for more than one relaxation time. This is the case in our model of multi-component mixtures.

Mechanical (diffusion) and thermal relaxation times were introduced in (25), while their influence on the source terms was described in (26). It was shown recently [53] that the rate of convergence may not be the same for different processes. Therefore, our aim is to compare different relaxation times and to get an insight to the rate of convergence towards equilibrium for different processes. We shall compare the relaxation times in a relative sense, i.e., by computing the ratio of relaxation times *τDbc* and *τTbc* to the reference relaxation time *τ*0 *<sup>D</sup>*<sup>12</sup> along the shock profiles. After straightforward calculation, one may find the explicit dimensionless form of ratios of the diffusion relaxation times to the reference one:

$$\begin{split} \frac{\tau\_{D11}}{\tau\_{D12}^{0}} &= \frac{\psi\_{12}^{0}}{\psi\_{11}} \frac{c\_{1}c\_{3}}{c\_{20}c\_{30}} \rho T = \widetilde{\tau}\_{D11}; \qquad \frac{\tau\_{D12}}{\tau\_{D12}^{0}} = -\frac{\psi\_{12}^{0}}{\psi\_{12}} \frac{c\_{2}c\_{3}}{c\_{20}c\_{30}} \rho T = \widetilde{\tau}\_{D12};\\ \frac{\tau\_{D21}}{\tau\_{D12}^{0}} &= -\frac{\psi\_{12}^{0}}{\psi\_{21}} \frac{c\_{1}c\_{3}}{c\_{20}c\_{30}} \rho T = \widetilde{\tau}\_{D21}; \quad \frac{\tau\_{D22}}{\tau\_{D12}^{0}} = \frac{\psi\_{12}^{0}}{\psi\_{22}} \frac{c\_{2}c\_{3}}{c\_{20}c\_{30}} \rho T = \widetilde{\tau}\_{D22}. \end{split}$$

Ratios of the thermal relaxation times to the reference one in dimensionless form read:

$$\begin{split} & \frac{\tau\_{T11}}{\tau\_{D12}^{0}} = \frac{a\_0^2 \psi\_{12}^0}{\theta\_{11}} \frac{\gamma - 1}{\gamma\_0 (\gamma\_1 - 1)(\gamma\_3 - 1)} \frac{m\_0 m}{m\_1 m\_3} \frac{c\_1 c\_3}{c\_{20} c\_{30}} \rho T^2 = \tilde{\tau}\_{T11}; \\ & \frac{\tau\_{T12}}{\tau\_{D12}^0} = -\frac{a\_0^2 \psi\_{12}^0}{\theta\_{12}} \frac{\gamma - 1}{\gamma\_0 (\gamma\_2 - 1)(\gamma\_3 - 1)} \frac{m\_0 m}{m\_2 m\_3} \frac{c\_2 c\_3}{c\_{20} c\_{30}} \rho T^2 = \tilde{\tau}\_{T12}; \\ & \frac{\tau\_{T21}}{\tau\_{D12}^0} = -\frac{a\_0^2 \psi\_{12}^0}{\theta\_{21}} \frac{\gamma - 1}{\gamma\_0 (\gamma\_1 - 1)(\gamma\_3 - 1)} \frac{m\_0 m}{m\_1 m\_3} \frac{c\_1 c\_3}{c\_{20} c\_{30}} \rho T^2 = \tilde{\tau}\_{T21}; \\ & \frac{\tau\_{T22}}{\tau\_{D12}^0} = \frac{a\_0^2 \psi\_{12}^0}{\theta\_{22}} \frac{\gamma - 1}{\gamma\_0 (\gamma\_2 - 1)(\gamma\_3 - 1)} \frac{m\_0 m}{m\_2 m\_3} \frac{c\_2 c\_3}{c\_{20} c\_{30}} \rho T^2 = \tilde{\tau}\_{T22}. \end{split}$$

Ratios of the phenomenological coefficients needed to complete these expressions are given in the Appendix A.

Relaxation times, presented in Figure 6, are computed for the profiles analyzed in Case 1. They share one common property—a decrease along the shock profile, although it may not be monotonic. However, the most remarkable feature of the relaxation times in the multi-component mixture, at least in this case, is that the thermal relaxation times *τ*˜*Tbc* are considerably smaller than the diffusion relaxation times *τ*˜*Dbc*. These results oppose the ones obtained in the case of binary mixture [40,47], even when the viscous and thermal dissipation are taken into account [53]. This certainly calls for deeper analysis of relaxation processes in the mixture. Additionally, it must be kept in mind that relaxation times may be strictly related to the rate of convergence towards equilibrium only in the simplest possible space-homogeneous cases. Therefore, relaxation times (25) in our model cannot be simply identified as the time-rates of convergence of non-equilibrium variables.

**Figure 6.** Relaxation times in the shock structure for *<sup>c</sup>*<sup>10</sup> = 0.10, *<sup>c</sup>*<sup>20</sup> = 0.75, *<sup>c</sup>*<sup>30</sup> = 0.15, and *<sup>M</sup>*<sup>0</sup> = 1.30: (**a**) mechanical (diffusion) relaxation times; (**b**) thermal relaxation times.

#### **6. Conclusions and Outlook**

This paper was devoted to the study of shock structure in multi-component mixture of Euler fluids. Analysis was based upon the hyperbolic MT model for mixtures developed within the framework of RET. Explicit form of phenomenological coefficients, which appear in the source terms, was determined by matching with a linearized weak form of the collision operator in the system of Boltzmann equations for mixtures of polyatomic gases. The analysis was restricted to the normal shock waves, smeared out into a shock structure due to dissipation. Therefore, the shock structure was assumed in the form of a traveling wave. The original governing equations are transformed into a system of ODE's that determine the shock structure as a heteroclinic orbit which asymptotically connects stationary points. Shock profiles were found as numerical solutions to the system of shock structure equations in dimensionless form, and analyzed for the several different sets of values of the parameters.

Our attention was focused on cases which possess certain distinguishing features, and thus can be taken as representatives of a broader class of qualitatively similar shock profiles. In particular, the following representative cases were recognized:


We also analyzed the ST model and compared its profiles with the MT ones. It was recognized that MT assumption does not produce significant difference between the profiles of the average temperatures in MT and ST case. However, certain profiles exist with larger magnitude of temperature overshoot, and systematic analysis of MT assumption in these cases is a matter of ongoing research.

Finally, to perceive the rate at which dissipation occurs, we computed the relaxation times which appear in the source terms. It was found out that relaxation times decreased along the shock profile (but not monotonically). Analysis also shed new light on their magnitudes by giving an estimate that thermal relaxation times could be smaller than the diffusion ones.

The present study opens a new chapter in analysis of the shock structure in gaseous mixtures. It is important to note that the model took advantage of the continuum approach, which secured thermodynamic consistency and peculiar structure of the source terms, and the kinetic theory, which equipped us with explicit form of the phenomenological coefficients. In such a way we obtain reliable and tractable mathematical model at the same time. Results obtained in this way are in good agreement with the results of other studies, and therefore confirm validity of the model. Further analysis, which is a work in progress, will be concerned with a systematic study of the shock structure in the regions of parameter space in which continuous solutions exist. This will give a deeper insight into the influence of certain parameters on quantitative characteristics of the shock profiles, like the temperature overshoot or the shock thickness. Additionally, in continuation, it is our intention to extend the analysis of the entropy growth and the entropy production rate, given in [52] for the binary mixture, to three-component mixtures. Our aim is also to take into account viscous and thermal dissipation and evaluate their influence on the shock structure in multi-component mixtures, like it was done in [53] for binary mixtures. Finally, it is our intent to analyze the cases of specific mixtures. To that end we plan to use novel models of polyatomic gases, derived in the kinetic theory [21], to get better estimates of phenomenological coefficients and to relate them to measurable physical quantities like diffusivity. This will also lead to possible explicit estimate of the shock thickness in terms of the average mean free path of the molecules.

**Author Contributions:** Conceptualization, M.P.-C. and S.S.; methodology, D.M., M.P.- ˇ C. and S.S.; ˇ software, D.M.; formal analysis, M.P.-C. and S.S.; investigation, D.M. and M.P.- ˇ C.; writing—original ˇ draft preparation, M.P.-C. and S.S.; writing—review and editing, D.M., M.P.- ˇ C. and S.S.; visualization, ˇ D.M.; funding acquisition, M.P.-C. and S.S. All authors have read and agreed to the published version ˇ of the manuscript.

**Funding:** This research was funded (D. Madjarevi´c and M. Pavi´c-Coli´ ˇ c) by the support of the Program for Excellent Projects of Young Researchers (PROMIS) of the Science Fund of the Republic of Serbia as members of the project MaKiPol #6066089, Mathematical methods in the kinetic theory of polyatomic gas mixtures: modelling, analysis and computation, whose support enabled publication of this article in an open access manner. The research was also supported (S. Simi´c) by the Ministry of Education, Science and Technological Development of the Republic of Serbia (Grant No. 451-03- 9/2021-14/200125).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **Appendix A**

Here we list the ratios of phenomenological coefficients which appear in the dimensionless form of the source terms (37) and (38). We use the abbreviated notation, *<sup>ψ</sup>bc* = *<sup>ψ</sup>bc*(**w**0), *<sup>θ</sup>bc* = *<sup>θ</sup>bc*(**w**0), and *<sup>ψ</sup>*<sup>0</sup> <sup>12</sup> denotes the phenomenological coefficient evaluated at upstream equilibrium. We also take into account that coefficients are symmetric, i.e., *<sup>ψ</sup>*<sup>21</sup> = *<sup>ψ</sup>*<sup>12</sup> and *<sup>θ</sup>*<sup>21</sup> = *<sup>θ</sup>*12.

Ratios of coefficients in (37) read:

$$\begin{aligned} \frac{\psi\_{11}}{\psi\_{12}^{0}} &= T \left[ \frac{K\_{12}^{\psi}}{K\_{12}^{\psi 0}} + \frac{\mu\_{13}}{\mu\_{12}} \frac{K\_{13}^{\psi}}{K\_{12}^{\psi 0}} \right], & \quad \frac{\psi\_{12}}{\psi\_{12}^{0}} &= -T \frac{K\_{12}^{\psi}}{K\_{12}^{\psi 0}},\\ \frac{\psi\_{21}}{\psi\_{12}^{0}} &= \frac{\psi\_{12}}{\psi\_{12}^{0}} = -T \frac{K\_{12}^{\psi}}{K\_{12}^{\psi 0}}, & \quad \frac{\psi\_{22}}{\psi\_{12}^{0}} &= T \left[ \frac{K\_{12}^{\psi}}{K\_{12}^{\psi 0}} + \frac{\mu\_{23}}{\mu\_{12}} \frac{K\_{23}^{\psi}}{K\_{12}^{\psi 0}} \right]. \end{aligned}$$

Ratios of coefficients in (38) read:

$$\begin{split} \frac{\theta\_{11}}{a\_{0}^{2}\psi\_{12}^{0}} &= \frac{3}{2} \frac{1}{\gamma\_{0}} \frac{m\_{0}}{\mu\_{12}} T^{2} \left[ \kappa\_{12}^{\theta} \frac{K\_{12}^{\theta}}{K\_{12}^{\theta 0}} + \kappa\_{13}^{\theta} \frac{K\_{13}^{\theta}}{K\_{12}^{\theta 0}} \right], & \quad \frac{\theta\_{12}}{a\_{0}^{2}\psi\_{12}^{0}} &= -\frac{3}{2} \frac{1}{\gamma\_{0}} \frac{m\_{0}}{\mu\_{12}} T^{2} \kappa\_{12}^{\theta} \frac{K\_{12}^{\theta}}{K\_{12}^{\theta 0}},\\ \frac{\theta\_{21}}{a\_{0}^{2}\psi\_{12}^{0}} &= \frac{\theta\_{12}}{a\_{0}^{2}\psi\_{12}^{0}} = -\frac{3}{2} \frac{1}{\gamma\_{0}} \frac{m\_{0}}{\mu\_{12}} T^{2} \kappa\_{12}^{\theta} \frac{K\_{12}^{\theta}}{K\_{12}^{\theta 0}}, & \quad \frac{\theta\_{22}}{a\_{0}^{2}\psi\_{12}^{0}} &= \frac{3}{2} \frac{1}{\gamma\_{0}} \frac{m\_{0}}{\mu\_{12}} T^{2} \left[ \kappa\_{12}^{\theta} \frac{K\_{12}^{\theta}}{K\_{12}^{\theta 0}} + \kappa\_{23}^{\theta} \frac{K\_{23}^{\theta}}{K\_{12}^{\theta 0}} \right]. \end{split}$$

Ratios of coefficients given above are expressed in terms of *K<sup>ψ</sup> bc* given in (22). Taking into account (22)–(24), one may find their explicit form:

$$\begin{split} \frac{K\_{12}^{\Psi}}{K\_{12}^{\Psi0}} &= \rho^2 \frac{c\_1 c\_2}{c\_{10} c\_{20}} T^{-(d\_1 + d\_2) + \frac{\nu\_{12}}{2}},\\ \frac{K\_{13}^{\Psi}}{K\_{12}^{\Psi0}} &= \rho^2 \frac{c\_1 c\_3}{c\_{10} c\_{20}} \frac{\kappa\_{13}^{\Psi}}{\kappa\_{12}^{\Psi}} \mu\_2 T^{-(d\_1 + d\_3) + \frac{\nu\_{13}}{2}},\\ \frac{K\_{23}^{\Psi}}{K\_{12}^{\Psi0}} &= \rho^2 \frac{c\_2 c\_3}{c\_{10} c\_{20}} \frac{\kappa\_{23}^{\Psi}}{\kappa\_{12}^{\Psi}} \mu\_1 T^{-(d\_2 + d\_3) + \frac{\nu\_{23}}{2}}. \end{split}$$

To facilitate numerical computation, we list below some useful expressions which appear in the model:

,

$$\begin{aligned} \frac{m\_0}{m} &= \left(\frac{c\_1}{\mu\_1} + \frac{c\_2}{\mu\_2} + c\_3\right) \left(\frac{c\_{10}}{\mu\_1} + \frac{c\_{20}}{\mu\_2} + c\_{30}\right)^{-1}, & \quad \frac{m\_0}{m\_1} = \frac{1}{\mu\_1} \left(\frac{c\_{10}}{\mu\_1} + \frac{c\_{20}}{\mu\_2} + c\_{30}\right)^{-1}, \\\frac{m\_0}{m\_2} &= \frac{1}{\mu\_2} \left(\frac{c\_{10}}{\mu\_1} + \frac{c\_{20}}{\mu\_2} + c\_{30}\right)^{-1}, & \quad \frac{m\_0}{m\_3} = \left(\frac{c\_{10}}{\mu\_1} + \frac{c\_{20}}{\mu\_2} + c\_{30}\right)^{-1}, \\\frac{m\_0}{\mu\_{12}} &= \frac{\mu\_1 + \mu\_2}{\mu\_1 \mu\_2} \left(\frac{c\_{10}}{\mu\_1} + \frac{c\_{20}}{\mu\_2} + c\_{30}\right)^{-1}, \\\frac{\mu\_{13}}{\mu\_{12}} &= \frac{\mu\_1}{\mu\_1 + 1} \frac{\mu\_1 + \mu\_2}{\mu\_1 \mu\_2}, & \quad \frac{\mu\_{23}}{\mu\_{12}} = \frac{\mu\_2}{\mu\_2 + 1} \frac{\mu\_1 + \mu\_2}{\mu\_1 \mu\_2}. \end{aligned}$$

#### **References**


## *Article* **Cyclic Control Optimization Algorithm for Stirling Engines**

**Raphael Paul \* and Karl Heinz Hoffmann**

Institute of Physics, Chemnitz University of Technology, 09107 Chemnitz, Germany; hoffmann@physik.tu-chemnitz.de

**\*** Correspondence: raphael.paul@physik.tu-chemnitz.de

**Abstract:** The ideal Stirling cycle describes a specific way to operate an equilibrium Stirling engine. This cycle consists of two isothermal and two isochoric strokes. For non-equilibrium Stirling engines, which may feature various irreversibilities and whose dynamics is characterized by a set of coupled ordinary differential equations, a control strategy that is based on the ideal cycle will not necessarily yield the best performance—for example, it will not generally lead to maximum power. In this paper, we present a method to optimize the engine's piston paths for different objectives; in particular, power and efficiency. Here, the focus is on an indirect iterative gradient algorithm that we use to solve the cyclic optimal control problem. The cyclic optimal control problem leads to a Hamiltonian system that features a symmetry between its state and costate subproblems. The symmetry manifests itself in the existence of mutually related attractive and repulsive limit cycles. Our algorithm exploits these limit cycles to solve the state and costate problems with periodic boundary conditions. A description of the algorithm is provided and it is explained how the control can be embedded in the system dynamics. Moreover, the optimization results obtained for an exemplary Stirling engine model are discussed. For this Stirling engine model, a comparison of the optimized piston paths against harmonic piston paths shows significant gains in both power and efficiency. At the maximum power point, the relative power gain due to the power-optimal control is ca. 28%, whereas the relative efficiency gain due to the efficiency-optimal control at the maximum efficiency point is ca. 10%.

**Keywords:** cyclic optimal control; finite time thermodynamics; endoreversible thermodynamics; Stirling; optimization

#### **1. Introduction**

Stirling engines are closed-cycle regenerative heat engines that harness a temperature difference between two external heat baths. At the cost of taking entropy from the hot heat bath and disposing it to the cold heat bath, they generate mechanical work. This means that Stirling engines are very flexible regarding the possible technical representations of those two heat baths. Therefore, they can be applied in various scenarios. Some current example applications are power generation from burnable industrial waste gases [1], domestic combined heat and power generation [2], electro-thermal energy storage systems for renewable energies as an alternative to lithium-ion batteries [3], and low-maintenance power generation in remote regions with harsh climatic conditions [4], among others. If properly designed, Stirling engines can achieve high efficiency, durability, low-maintenance, and low-noise operation. Hence, they are considered to be one good candidate technology for enhancing the sustainability of power systems.

One certain ideal representation of the thermodynamic cycle performed in Stirling engines is commonly referred to as the ideal Stirling cycle. It consists of four distinct strokes (processes) that an enclosed working gas cyclically performs:


**Citation:** Paul, R.; Hoffmann, K.H. Cyclic Control Optimization Algorithm for Stirling Engines. *Symmetry* **2021**, *13*, 873. https:// doi.org/10.3390/sym13050873

Academic Editors: Robert Kovacs, Patrizia Rogolino and Francesco Oliveri

Received: 15 April 2021 Accepted: 9 May 2021 Published: 13 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

where *T*ExH and *T*ExL are the temperatures of the hot and cold heat bath, respectively. If those four ideal processes are repeated in the prescribed order, the cycle represents an ideal heat engine, operating with Carnot efficiency *<sup>η</sup>*<sup>C</sup> <sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>T</sup>*ExL/*T*ExH. In the two regenerative processes 4 and 2, heat needs to be stored and provided at a range of temperatures. The combined process is called regeneration, and is accomplished with help of a multi-temperature heat storage, referred to as a regenerator.

In real Stirling engines, the working gas is typically contained in two working spaces. One of them, the hot working space, is in thermal contact with the hot heat bath. The other one, the cold working space, is in thermal contact with the cold heat bath. The two working spaces' volumes can be changed by moving pistons. As a result of those pistons' movements, the working spaces exchange working gas through the regenerator, which is often implemented as a porous metal structure.

In contrast to the ideal Stirling cycle, in real Stirling engines, the four above-described strokes can typically not be clearly distinguished. This is due to several reasons. First of all, the engine's piston movements (piston paths) usually do not exactly reproduce the two isochoric strokes. Moreover, the regenerator, as well as the hot and cold working spaces, have non-zero dead volumes. The most important reason, however, is connected to the fact that real engines are supposed to produce a finite amount of work in finite time, which is in contrast to the assumptions of equilibrium thermodynamics, and thus to the ideal Stirling cycle. Consequently, in real Stirling engines, various inevitable loss phenomena occur, for example, finite heat transfer between the working gas and the external heat baths or the regenerator matrix, pressure drop across the regenerator, and heat leaks through the engine's internal components.

Correspondingly, real Stirling engines are non-equilibrium devices. Hence, piston paths aiming to emulate the ideal Stirling cycle's four strokes would not necessarily yield best performance. The question of what piston paths yield best performance, for example, optimal net power, constitutes a cyclic optimal control problem.

In this paper we present an algorithm that can be used to solve such cyclic optimal control problems. It is an indirect iterative gradient algorithm and based on Optimal Control Theory. The algorithm starts with prescribing initial control functions determining the piston paths. These control functions are then gradually modified over the course of the iterations so as to gradually improve an objective and approach the optimal control. To determine the gradual shifts in every iteration, not only does the system of ordinary differential equations, describing the thermodynamics of the Stirling engine, need to be solved, but also a conjugate system of differential equations. For the efficient practical feasibility of this task, low numerical effort of the utilized thermodynamic model is crucial. In detailed models of Stirling engines, significant numerical effort is usually connected to the description of the regenerator. Therefore, here we will apply a Stirling engine model with a reduced-order endoreversible regenerator model that provides a proper tradeoff between accuracy and numerical effort for optimal control problems.

As indicated above, the core discrepancy between the ideal Stirling cycle and a real Stirling engine is that, in the former, all processes are quasi-static, and hence reversible, whereas the latter is required to produce a finite amount of work in finite time. Therefore, real Stirling engines are non-equilibrium devices that could by no means achieve the performance measures of the ideal cycle, such as Carnot efficiency. A non-equilibrium thermodynamics field that studies how constraints on time and rate influence the performance of such devices is Finite-Time Thermodynamics (FTT) [5,6]. FTT typically approaches this and related questions with variational principles and global—rather than local—descriptions of the irreversible systems considered [5].

Endoreversible Thermodynamics [7–11] can be regarded as a subfield of FTT that intends to provide a toolbox of model building blocks from which an FTT model can be constructed in a comprehensible way. The emphasis is placed on including the main loss phenomena, while keeping the model structure clear and minimizing mathematical complexity and numerical effort. The basic approach is to decompose the considered thermodynamic system into a network of reversible subsystems that are connected with each other through (ir)reversible interactions. That is, in its endoreversible representation, the original system's irreversibilities are captured by the interactions, whereas for the reversible subsystems the known relations from equilibrium thermodynamics hold.

In the past, various endoreversible systems have been subject to investigations involving control optimizations. In particular, for heat engines, the piston trajectories can be optimized for objectives, such as power or efficiency. For example, such studies have been carried out for engines with Diesel [11–14] and Otto [15–17] cycles, light-driven engines [18–20], and Stirling engines [21–23].

The optimization approach used in the current study to optimize the Stirling engine differs from [21–23] in that a cyclic version of Optimal Control Theory is applied. Beyond Endoreversible Thermodynamics, this has been done by Craun and Bamieh [24] for an ideal beta-Stirling engine model with an actively controlled displacer. In [24], the heat transfer between the heat baths and the working gas was infinitely fast and regeneration was perfect. The pressure drop across the regenerator was thermodynamically neglected, but considered in terms of mechanical losses. As opposed to this, in the current study, an alpha-Stirling engine model with finite heat transfer between the heat baths and the working gas, an external heat leak, gas leakages and friction of the piston rings, and an irreversible regenerator is considered. Compared to the endoreversible regenerator models used in [22,23], the *EEn*-regenerator model [25] used in this study features a significantly higher degree of detail.

#### **2. Stirling Engine Model**

In this paper, our focus is on a cyclic control optimization algorithm that can be applied to various Stirling engine models, or models of other cyclically operating systems. Nevertheless, we here briefly introduce a Stirling engine model [25] that we will use to demonstrate the algorithm's usage. In particular, we will present optimization results for this model in Section 5.

We consider an alpha-type Stirling engine as schematized in Figure 1, where the heat transfer between the external heat baths and the working gas is realized through the cylinders.

**Figure 1.** Schematics of the alpha-type Stirling engine considered. The corresponding endoreversible model [25] has ten degrees of freedom (state variables), being the working space's volumes *Vi*, temperatures *Ti*, and particle numbers *ni*, the energies *E*R.h and *E*R.l of the hot and cold halves of the regenerator matrix as well as the particle number *n*R.d and temperature *T*R.d of the gas inside the regenerator dead space.

For the sake of simplicity, we assume a constant finite heat transfer coefficient. Further irreversibilities result from an external heat leak between the heat baths, gas leakage of the piston rings, friction of the pistons, pressure drop across the regenerator, and thermal mixing at the interfaces of the regenerator and the working spaces.

The regenerator is described with the reduced-order endoreversible *EEn*-regenerator model, developed in [25]. Essentially, this is based on the assumptions that the spatial temperature distribution inside the regenerator is linear and that the temperature difference between the gas and matrix is small. Aside from external irreversibilities that automatically occur due to thermal mixing in the regenerator's interactions with the working spaces, internal irreversibilities can be included via entropy source terms in the *EEn*-regenerator model. In the current study, the irreversibility due to the pressure drop across the regenerator is accounted for by such an entropy source term.

The Stirling engine performance measures that are to be optimized in this study are power and efficiency. The power, or net power output, is defined as

$$P = \frac{W}{\pi r},\tag{1}$$

with the fixed cycle time *τ* and the net cycle work

$$\mathcal{W} = \int\_0^\tau \left( p\_\text{H} \, \dot{V}\_\text{H} + p\_\text{L} \, \dot{V}\_\text{L} \right) - \gamma \left( \dot{V}\_\text{H}^2 + \dot{V}\_\text{L}^2 \right) \, \text{d}t,\tag{2}$$

where *Vi* and *pi* are the volume and pressure of the working spaces *i* ∈ {H, L} and *γ* is the friction coefficient of the pistons. Then, the efficiency is

$$\eta = \frac{\int\_0^\tau \left(p\_\text{H} \,\dot{V}\_\text{H} + p\_\text{L} \,\dot{V}\_\text{L}\right) - \gamma \left(\dot{V}\_\text{H}^2 + \dot{V}\_\text{L}^2\right) \,\text{d}t}{\int\_0^\tau I\_\text{H} \,\text{ExH} + K\_\text{Ex} \,\left(T\_\text{ExH} - T\_\text{ExL}\right) \,\text{d}t},\tag{3}$$

where *I*H,ExH is the instantaneous heat flux to the hot working space and *K*Ex (*T*ExH − *T*ExL) corresponds to an external heat leak. We here use the same transfer laws and parameter values as in [25], for example, the heat bath temperatures are defined as *<sup>T</sup>*ExH <sup>=</sup> 500 K and *<sup>T</sup>*ExL <sup>=</sup> 300 K. Solely for the heat conductance of the external heat leak, we chose a different parameter value: *<sup>K</sup>*Ex <sup>=</sup> 5 W/K.

This Stirling engine model has ten state variables, which are the working spaces' volumes *Vi*, temperatures *Ti*, and particle numbers *ni* with *i* ∈ {H,L}, the energies *E*R.h and *E*R.l of the hot and cold halves of the regenerator matrix as well as the particle number *n*R.d and temperature *T*R.d of the gas inside the regenerator dead space.

The state variables of this Stirling engine model can be arranged in a state vector *x* and the state dynamics can be expressed in terms of *x*˙ = *f*(*x*, *u*). Here, *f*(*x*, *u*) is a vectorvalued function that depends on the state vector *x* and a control vector *u*. In our Stirling engine model we define the control vector *u* in terms of an explicitly time-dependent, *<sup>τ</sup>*-periodic control function *<sup>u</sup>*(*t*), which has the entries *<sup>u</sup>*H(*t*) and *<sup>u</sup>*L(*t*). Each of those two sub-functions determines the dynamics of one working volume, as indicated in Figure 1. An essential requirement that we raise is that the state dynamics *x*˙ = *f*(*x*, *u*(*t*)) features a limit cycle with regard to all state variables, which will be necessary for the applicability of the optimization algorithm. Therefore, care must be taken regarding the embedding of the control in the system dynamics, which will be explained in Section 4.

#### **3. Cyclic Optimal Control Problem**

Optimal Control Theory provides a framework for the dynamic (or indirect) optimization of dynamical systems, such as the above-described Stirling engine model. The type of optimal control problem considered in this work can be formulated as: find an unconstrained, smooth, *τ*-periodic control function *u*(*t*) that maximizes the objective functional

$$J = \int\_0^\pi \zeta(\mathbf{x}, \mathbf{u}(t)) \, \mathrm{d}t \tag{4}$$

subject to the constraint

$$
\dot{\mathbf{x}} = f(\mathbf{x}, \boldsymbol{\mu}(t)) \tag{5}
$$

for a fixed period *τ*. Here, *x* is the state vector of the system. The dynamics of the system is determined by the ordinary differential Equation (5) and is influenced by the vector-valued, explicitly time-dependent, *τ*-periodic control function *u*(*t*).

For the formulation of necessary conditions of optimality, a Hamilton function *H* is defined as [26,27]

$$H(\mathbf{x}, \boldsymbol{\mu}, \boldsymbol{\lambda}) := \mathbb{\mathcal{J}}(\mathbf{x}, \boldsymbol{\mu}) + \boldsymbol{\lambda}^T f(\mathbf{x}, \boldsymbol{\mu}) \tag{6}$$

with an adjoint vector *λ*, which varies over time and will henceforth also be called a costate vector. As we will see later, there exists an interesting symmetry between the dynamics of the state vector *x* and the costate vector *λ*. The first order necessary conditions for the control *u* with the optimal value are [26,27]

$$
\dot{\mathbf{x}} = \nabla\_{\lambda} H(\mathbf{x}, \mathbf{u}, \lambda),
\tag{7}
$$

$$
\dot{\lambda} = -\nabla\_{\mathbf{x}} H(\mathbf{x}, \mu, \lambda),
\tag{8}
$$

$$\mathbf{0} = \nabla\_{\mathbf{u}} H(\mathbf{x}, \mathbf{u}, \lambda), \tag{9}$$

combined with the respective transversality conditions that specify the boundary conditions for the state and costate dynamics from Equations (7) and (8), respectively. In the case of the optimal cyclic regime, as valid for the stationary operational state of a cyclically operating engine, these boundary conditions are [27]

$$\left.\mathbf{x}\right|\_{\tau} = \left.\mathbf{x}\right|\_{0^{\prime}}\tag{10}$$

$$
\lambda|\_{\tau} = \lambda|\_{0}.\tag{11}
$$

The notation *x*|*<sup>t</sup>* refers to the value of *x* at time *t*, while it does not indicate explicit time dependence. These boundary conditions do not render the problem overdetermined, since the initial and final values are not quantified, but only required to be equal.

In this way, the global dynamic optimization problem from Equation (4) with the constraint Equation (5) is turned into the continuous set of local (static) optimization problems from Equation (9) with the constraints Equations (7) and (8). This means that *H*(*x*, *u*, *λ*) must be locally maximized for all *<sup>t</sup>* <sup>∈</sup> [0, *<sup>τ</sup>*) with respect to *<sup>u</sup>*. For the optimal solution *<sup>x</sup>*<sup>∗</sup> , *u*∗ , *λ*∗ of the considered kind of cyclic optimization problem, without explicit time dependence of *ζ*(*x*, *u*), the temporal value of the Hamilton function *H*<sup>∗</sup> will then be constant [26,27].

#### *3.1. Maximum Power*

In the case of the optimization being performed for the objective functional, which is the cycle-averaged net power output *P*, as defined in Equation (1), the definition of the path target function *ζ*(*x*, *u*) for Equations (4) and (6) is straightforward:

$$\mathcal{Z}\_{\rm W}(\mathbf{x}, \boldsymbol{\mu}) := \left( p\_{\rm H} \dot{\boldsymbol{V}}\_{\rm H} + p\_{\rm L} \dot{\boldsymbol{V}}\_{\rm L} \right) - \gamma \left( \dot{\boldsymbol{V}}\_{\rm H}^{2} + \dot{\boldsymbol{V}}\_{\rm L}^{2} \right), \tag{12}$$

which is simply the instantaneous power output. Since the optimization is performed for fixed cycle time *τ*, maximum power corresponds to maximum work. Hence, a prefactor of 1/*τ* can be disregarded here.

#### *3.2. Maximum Efficiency*

In the case of the optimization being performed for the objective of efficiency *η*, as defined in Equation (3), the definition of *ζ*(*x*, *u*) is slightly more involved. This is because the definition of the efficiency corresponds to a ratio of two functionals (integrals) which need to be separately evaluated before calculating their ratio. This can be seen in Equation (3) and is in contrast to the structure of Equation (4). In order to resolve this discrepancy, we investigate how variations in the cycle work *W* and the heat *<sup>Q</sup>*<sup>H</sup> :<sup>=</sup> # *<sup>τ</sup>* <sup>0</sup> *<sup>I</sup>*H,ExH <sup>+</sup> *<sup>K</sup>*Ex (*T*ExH <sup>−</sup> *<sup>T</sup>*ExL) <sup>d</sup>*<sup>t</sup>* translate to the first variation of the efficiency:

$$
\delta\eta = \delta \left(\frac{W}{Q\_{\rm H}}\right) = \frac{\partial\eta}{\partial W}\,\delta W + \frac{\partial\eta}{\partial Q\_{\rm H}}\,\delta Q\_{\rm H} = \frac{1}{Q\_{\rm H}}\,\delta W - \frac{W}{Q\_{\rm H}^2}\,\delta Q\_{\rm H}.\tag{13}
$$

This means that the problem of optimizing the efficiency can be morphed into an equivalent optimization problem of the weighted sum *<sup>λ</sup>*<sup>W</sup> *<sup>W</sup>* <sup>+</sup> *<sup>λ</sup>*QH *<sup>Q</sup>*<sup>H</sup> with the weights being defined according to the partial derivatives:

$$
\lambda\_{\mathcal{W}} = \partial \eta / \partial \mathcal{W} = 1 / Q\_{\mathcal{H}} \text{ and } \begin{aligned}
\lambda\_{\mathcal{Q}\mathcal{H}} = \partial \eta / \partial Q\_{\mathcal{H}} = -\mathcal{W} / Q\_{\mathcal{H}}^2. \end{aligned} \tag{14}$$

Within these weights, *W* and *Q*<sup>H</sup> are the results obtained for a specific prescribed *τ*-periodic control function *u*† (*t*) for which, at a given time *t*, the values of the control and state vectors are *<sup>u</sup>* <sup>≈</sup> *<sup>u</sup>*† (*t*) and *<sup>x</sup>* <sup>≈</sup> *<sup>x</sup>*† (*t*). Here, *x*† (*t*) is the periodic solution of the state dynamics for *u*† (*t*). Based on this, we can now define the path target function for maximizing the efficiency:

$$\mathcal{J}\_{\eta}(\mathbf{x},\boldsymbol{\mu}) := \frac{\boldsymbol{W}}{Q\_{\mathrm{H}}\,\boldsymbol{\tau}} + \frac{1}{Q\_{\mathrm{H}}} \mathcal{J}\_{\mathrm{W}}(\mathbf{x},\boldsymbol{\mu}) - \frac{\boldsymbol{W}}{Q\_{\mathrm{H}}^{2}} I\_{\mathrm{H},\mathrm{ExH}}(\mathbf{x},\boldsymbol{\mu}).\tag{15}$$

Here, *W* and *Q*<sup>H</sup> are again the results obtained for a *τ*-periodic control function *u*† (*t*) with the above-described closeness requirements for the instantaneous values of *x* and *<sup>u</sup>*. In a strict sense, the path target function *ζη*(*x*, *<sup>u</sup>*) is therefore also a functional of *u*† (*t*) and *x*† (*t*), that is: *ζη*(*x*, *<sup>u</sup>*) = *ζη*[ *<sup>x</sup>*† (·), *<sup>u</sup>*† (·)](*x*, *<sup>u</sup>*) because *<sup>W</sup>* <sup>=</sup> *<sup>W</sup>*[ *<sup>x</sup>*† (·), *<sup>u</sup>*† (·)] and *<sup>Q</sup>*<sup>H</sup> <sup>=</sup> *<sup>Q</sup>*H[ *<sup>x</sup>*† (·), *<sup>u</sup>*† (·)]. If the closeness requirement of *<sup>x</sup>* and *<sup>u</sup>*, regarding *<sup>u</sup>*† (*t*) is valid for all *t*, the second and third terms in Equation (15) will (approximately) cancel upon integration from 0 to *τ*. Therefore, the first term was added here to ensure that # *<sup>τ</sup>* <sup>0</sup> *ζη*(*x*, *<sup>u</sup>*) <sup>d</sup>*<sup>t</sup>* does actually represent *η*. This first term does, in fact, not influence the optimization problem from Equation (7) to Equation (11), since it does not depend on the instantaneous values of *x* and *u*.

#### *3.3. Penalty Function*

The cyclic optimal control problem introduced so far does not contain inequality constraints for the state *x* or the control *u*. However, we can introduce inequality constraints for *x* and *u* by means of a penalty function. In particular, we will use such a penalty function to introduce lower and upper bounds for the Stirling engine's working volumes:

$$\mathcal{Z}\_{\text{penalty}}(\mathbf{x},\boldsymbol{\mu}) := -v\_0 \left( e^{v\_1 \frac{V\_{\text{H,min}} - V\_{\text{H}}}{V\_{\text{H,sw}}}} + e^{v\_1 \frac{V\_{\text{H}} - V\_{\text{H,max}}}{V\_{\text{H,sw}}}} + e^{v\_1 \frac{V\_{\text{L,min}} - V\_{\text{L}}}{V\_{\text{L,sw}}}} + e^{v\_1 \frac{V\_{\text{L}} - V\_{\text{L,max}}}{V\_{\text{L,sw}}}} \right) \tag{16}$$

where the working volumes *Vi* are contained in *<sup>x</sup>*, *<sup>v</sup>*<sup>0</sup> <sup>=</sup> 1W and *<sup>v</sup>*<sup>1</sup> <sup>=</sup> 500 are penalty factors, *Vi*,sw <sup>=</sup> *Vi*,max <sup>−</sup> *Vi*,min refers to the admissible swept volumes and *Vi*,max and *Vi*,min are the maximum and minimum admissible working volumes with *i* ∈ {H, L}.

Accordingly, in the case of optimizing the Stirling engine's piston paths for maximum power, the overall path target function is defined as

$$\mathcal{L}(\mathbf{x}, \mathfrak{u}) := \mathcal{J}\_W(\mathfrak{x}, \mathfrak{u}) + \mathcal{J}\_{\text{penalty}}(\mathfrak{x}, \mathfrak{u})\_\prime \tag{17}$$

whereas in case of the optimization for maximum efficiency, it is defined as

$$\mathcal{Z}(\mathbf{x}, \mathfrak{u}) := \mathcal{Z}\_{\eta}(\mathbf{x}, \mathfrak{u}) + \mathcal{Z}\_{\text{peralty}}(\mathbf{x}, \mathfrak{u}). \tag{18}$$

#### **4. Optimization Algorithm**

The indirect iterative gradient method described in the following was inspired by a contribution of Craun and Bamieh [24,28], who solved an optimal control problem of a modified beta-Stirling engine using an ideal thermodynamic model, namely the Schmidt model. Similar methods have earlier been used, for example, by Horn and Lin [29], Kowler and Kadlec [30], as well as Noorden et al. [31] for other applications. The method presented in this paper differs from the aforementioned ones in the means by which the periodic solutions of the state and costate equations are obtained. A pseudocode representation of the cyclic control optimization algorithm is shown in Algorithm 1.

#### **Algorithm 1** Cyclic optimal control algorithm.

```
1: //############ Basic functions of the control problem ############
2: define f(x, u){ // RHS of state dynamics
3: return temporal rate of state vector;
4: };
5: define ζ(x, u, ¯a){ // Path cost function, ¯a may contain W and QH
6: return path cost value;
7: };
8: define H(x, u, λ, ¯a){ // Hamilton function
9: return ζ(x, u, ¯a) + λT f(x, u);
10: };
11: define -
          ∇app
            x H(x, u, λ, ¯a){ // Approximation of RHS of costate dynamics
12: declare ε = small number;
13: declare eˆm = unit vector of mth state component;
14: return ∑m -
                eˆm( H(x + ε eˆm, u, λ, ¯a) − H(x, u, λ, ¯a) )/ ε;
15: };
16: define ∇app
           u H(x, u, λ, ¯a){ // Approximation of control gradient function
17: declare ε = small number;
18: declare eˆm = unit vector of mth control component;
19: return ∑m eˆm( H(x, u + ε eˆm, λ, ¯a) − H(x, u − ε eˆm, λ, ¯a) )/(2 ε);
20: };
21: //########## Declare control, state, and costate functions ###########
22: declare τ = cycle time; // Define cycle time
23: declare u(t) with t ∈ [0, τ); // Control vector function
24: declare x(t) with t ∈ [0, τ); // State vector function
25: declare λ(t) with t ∈ [0, τ); // Costate vector function
26: initialize u(t) with smooth guess of optimal control vector function;
27: initialize x(0) with arbitrary admissible values;
28: initialize λ(0) = 0;
29: //############ Iterative optimization of the control #############
30: for(n = 0; n < large number; n = n + 1){
31: while(ε > small number){ // Solve state dynamics
32: for(t = 0; t < τ; t += Δt){ x(t + Δt) = x(t) + f(x(t), u(t)) Δt; }
33: ε = ||x(τ) − x(0)||;
34: x(0) = x(τ);
35: }
36: declare a(t) = some_function(x(t),u(t)) for t ∈ [0, τ); // Auxiliary vector function
37: declare ¯a = cycle average of a(t) for t ∈ [0, τ); // Auxiliary vector cycle average
38: while(ε > small number){ // Solve costate dynamics
39: λ(τ) = λ(0);
40: for(t = τ; t > 0; t −= Δt){ λ(t − Δt) = λ(t) − -
                                                   ∇app
                                                     x H(x(t), u(t), λ(t), ¯a) Δt; }
41: ε = ||λ(τ) − λ(0)||;
42: }
43: declare s(t) = ∇app
                    u H(x(t), u(t), λ(t)) for t ∈ [0, τ); // Define search direction
44: smooth s(t) in the periodic domain t ∈ [0, τ); // Smooth search direction
45: declare Λ = small number; // Step-size factor
46: update u(t) = u(t) + Λ s(t) for t ∈ [0, τ); // Shift control
47: smooth u(t) in the periodic domain t ∈ [0, τ); // Smooth control
48: }
```
In the first section, between lines 1 and 20, the basic function definitions are made. First, the right-hand side (RHS) of the state dynamics is defined as a vector function *f*(*x*, *u*). Afterwards, the path cost function *ζ*(*x*, *u*, *a*¯) that is to be maximized in terms of the objective functional from Equation (4) is defined. Here, the argument *a*¯ may, for example, contain the cycle work *W* and heat *Q*H, as required in Equation (15). Then, based on the definition of the Hamilton function *H*(*x*, *u*, *λ*, *a*¯) from lines 8 to 10, the right-hand side of the costate dynamics −∇app *<sup>x</sup> <sup>H</sup>*(*x*, *<sup>u</sup>*, *<sup>λ</sup>*, ¯*a*) is defined in terms of an approximation:

$$-\frac{\partial H(\mathbf{x}, \boldsymbol{\mu}, \boldsymbol{\lambda}, \bar{\mathbf{a}})}{\partial \mathbf{x}\_m} \approx -\frac{H(\mathbf{x} + \varepsilon \, \mathfrak{E}\_m, \boldsymbol{\mu}, \boldsymbol{\lambda}, \bar{\mathbf{a}}) - H(\mathbf{x}, \boldsymbol{\mu}, \boldsymbol{\lambda}, \bar{\mathbf{a}})}{\varepsilon} \tag{19}$$

where *xm* is the *m*th state component, *e***ˆ***<sup>m</sup>* is the corresponding unit vector and *ε* is a sufficiently small difference. For simple systems, it is reasonable to derive this right-hand side in terms of algebraic expressions. However, for involved systems it can be much more practical to do this by numerical partial differentiation as, for example, according to Equation (19). Note that for systems with discontinuous dynamics additional care must be taken here in order to make sure that the derivatives are properly approximated in the vicinity of the discontinuity. Similarly to −∇app *<sup>x</sup> <sup>H</sup>*(*x*, *<sup>u</sup>*, *<sup>λ</sup>*, *<sup>a</sup>*¯), the vector <sup>∇</sup>app *<sup>u</sup> <sup>H</sup>*(*x*, *<sup>u</sup>*, *<sup>λ</sup>*, *<sup>a</sup>*¯) of partial derivatives of the Hamilton function with respect to the control is defined from lines 16 to 20. Here, the following approximation is used:

$$\frac{\partial H(\mathbf{x}, \boldsymbol{\mu}, \boldsymbol{\lambda}, \bar{\mathbf{a}})}{\partial \boldsymbol{u}\_{m}} \approx \frac{H(\mathbf{x}, \boldsymbol{\mu} + \varepsilon \,\mathfrak{E}\_{m}, \boldsymbol{\lambda}, \bar{\mathbf{a}}) - H(\mathbf{x}, \boldsymbol{\mu} - \varepsilon \,\mathfrak{E}\_{m}, \boldsymbol{\lambda}, \bar{\mathbf{a}})}{2\,\varepsilon} \tag{20}$$

where *um* is the *m*th control component, *e***ˆ***<sup>m</sup>* is the corresponding unit vector and *ε* is, again, a sufficiently small difference.

In the second section between lines 21 and 28, the most relevant variables are declared and initialized. These are the control vector function *u*(*t*), the state vector function *x*(*t*), and the costate vector function *λ*(*t*). They are defined in the time domain [0, *τ*) and may be implemented as vectors of arrays, where the array lengths correspond to the number of time steps per period. That means that the values of *u*, *x*, and *λ* are saved at every time step in the periodic domain. The control vector function *u*(*t*) is initialized with a smooth *τ*-periodic guess of the optimal control. Moreover, *x*(0) is set to arbitrary but physically admissible values, whereas *λ*(0) is set to **0**.

The iterative control optimization itself is described in the third section between lines 29 and 48. This starts at iteration *n* = 0 with the initial guess of the smooth *τ*-periodic control function *<sup>u</sup>*(0) (*t*), and the previously defined state and costate initial values *<sup>x</sup>*(0) (0) and *<sup>λ</sup>*(0) (0), respectively. Here, the bracketed superscript identifies the iteration *n* of the forloop from lines 30 to 48. Since the control is predefined, the state and costate problem can be solved separately. The steps to obtain the subsequent iteration *<sup>u</sup>*(*n*+1) (*t*), according to Algorithm 1, can be summarized in simplified terms:

0. If *<sup>n</sup>* <sup>=</sup> 0, set *<sup>x</sup>*(0) (0) to arbitrary physically admissible values and *<sup>λ</sup>*(0) (0) := **0**. Otherwise, set *<sup>x</sup>*(*n*) (0) :<sup>=</sup> *<sup>x</sup>*(*n*−1) (*τ*) and *<sup>λ</sup>*(*n*) (0) :<sup>=</sup> *<sup>λ</sup>*(*n*−1) (*τ*).

1. For given *<sup>u</sup>*(*n*) (*t*): Solve

$$\dot{\mathbf{x}}^{(n)}(t) = f(\mathbf{x}^{(n)}(t), \boldsymbol{\mu}^{(n)}(t))$$

by temporal forward integration in the periodic domain until the cyclic equilibrium is reached. Save auxiliary quantities, such as the cycle work *W* and the heat *Q*<sup>H</sup> in *a*¯ (*n*) . (Algorithm 1, lines 31 to 37: example with Euler method.)

2. For given *<sup>u</sup>*(*n*) (*t*), *<sup>x</sup>*(*n*) (*t*), ¯*a*(*n*) : Solve

$$\dot{\boldsymbol{\lambda}}^{(n)}(t) = -\nabla\_{\mathbf{x}} H(\mathbf{x}^{(n)}(t), \boldsymbol{\mu}^{(n)}(t), \boldsymbol{\lambda}^{(n)}(t), \bar{\boldsymbol{\mu}}^{(n)}),$$

by temporal backward integration in the periodic domain until the cyclic equilibrium is reached. (Algorithm 1, lines 38 to 42: example with Euler method.)

3. Calculate the search direction (See Algorithm 1, line 43.)

$$\mathbf{s}^{(n)}(t) := \nabla\_{\
u} H(\mathbf{x}^{(n)}(t), \mathbf{u}^{(n)}(t), \mathbf{\lambda}^{(n)}(t)).$$

4. Calculate the next control vector function

$$\mathfrak{u}^{(n+1)}(t) := \mathfrak{u}^{(n)}(t) + \Lambda \mathfrak{s}^{(n)}(t)$$

with a sufficiently small positive step-size factor Λ. (Algorithm 1, line 46.)

These five steps are repeated for increasing *<sup>n</sup>* until *<sup>u</sup>*(*n*) (*t*) has converged to the optimal control. This can, for example, be checked with the below indicators which should both approach a numerical zero:

$$\mathbb{C}\_{\nabla}^{(n)} := \int\_0^{\tau} \left( \nabla\_{\mathfrak{u}} H(\mathbf{x}^{(n)}(t), \mathfrak{u}^{(n)}(t), \mathsf{A}^{(n)}(t)) \right)^2 \mathrm{d}t,\tag{21}$$

$$\mathbb{C}\_{\mathbb{H}}^{(n)} := \max\_{t} H(\mathbf{x}^{(n)}(t), \boldsymbol{\mu}^{(n)}(t), \boldsymbol{\lambda}^{(n)}(t)) - \min\_{t} H(\mathbf{x}^{(n)}(t), \boldsymbol{\mu}^{(n)}(t), \boldsymbol{\lambda}^{(n)}(t)). \tag{22}$$

#### *4.1. Exploiting Limit Cycles*

The dynamic systems considered here are required to posses a limit cycle. This means that, given a proper periodic control function, the systems feature a closed state trajectory *<sup>x</sup>*∞(*t*) in the state space, which is attractive. Here, the subscript <sup>∞</sup> refers to this limit cycle. If the state dynamics is integrated forward in time starting from *t* = 0 with a proper initial value *<sup>x</sup>*|<sup>0</sup> in the vicinity of *<sup>x</sup>*∞(0), then *<sup>x</sup>*|*<sup>t</sup>* <sup>→</sup> *<sup>x</sup>*∞(*t*) for *<sup>t</sup>* <sup>→</sup> <sup>∞</sup>. This is a property of many dissipative systems, and it is directly made use of in the first step (Algorithm 1, lines 31 to 35) for finding the periodic solution of the state vector.

Numerical studies of optimal control problems carried out in the context of this work show an interesting type of symmetry in that the costate dynamics then features a limit cycle, too, for prescribed *<sup>u</sup>*(*t*) and *<sup>x</sup>*|*<sup>t</sup>* :<sup>=</sup> *<sup>x</sup>*∞(*t*). However, the respective closed trajectory *<sup>λ</sup>*∞(*t*) is repulsive when the costate dynamics is considered forward in time, and becomes attractive in the reversed time direction: *<sup>λ</sup>*|*<sup>t</sup>* <sup>→</sup> *<sup>λ</sup>*∞(*t*) for *<sup>t</sup>* <sup>→</sup> -∞. Such a relation is not surprising, as the overall Hamiltonian system is conservative and the expansive costate dynamics compensates for the contractive state dynamics. Therefore, in the second step (Algorithm 1, line 38 to 42), the temporal direction of integration is reversed to obtain the periodic solution of the costate vector.

The correspondingly enhanced stability of the costate problem under temporal backward integration can also be exploited without directly making use of the limit-cycle property for obtaining the periodic solution, but for solving boundary value problems, for example, see [26,31].

#### *4.2. Proper Embedding of the Control*

Note that the way in which the control is embedded in the system's state dynamics, can influence the existence of these limit cycles in the state and costate problems. In the Stirling engine optimization problem, considered in this paper, our goal is to optimize the piston paths. Hence, we chose to structure the system dynamics in such a way that the control function determines the temporal evolution of the working space volumes *Vi* with *i* ∈ {H, L}. A rather obvious choice for embedding the control *ui*(*t*) would be:

$$
\dot{V}\_i = u\_i(t). \tag{23}
$$

For an explicitly time-dependent, *<sup>τ</sup>*-periodic control *ui*(*t*) that has an average value # *τ* <sup>0</sup> *ui*(*t*) <sup>d</sup>*<sup>t</sup>* <sup>=</sup> 0, the resulting volumes *Vi* will be *<sup>τ</sup>*-periodic. However, this dynamics does not feature a limit cycle: The solution of *Vi*|*<sup>t</sup>* will even for *t* being large integer multiples of *τ* remain equal to the initial value *Vi*|<sup>0</sup> from which integration was started. This will translate to the costate problem, so that the solution algorithm from above will, in general, not converge without further manipulation of the optimization problem.

A proper choice for embedding the control in the system's state dynamics that leads to the existence of a limit cycle for a continuous *<sup>τ</sup>*-periodic *ui*(*t*) is:

$$
\dot{V}\_i = \frac{\nu}{\tau} (\mu\_i(t) - V\_i) \tag{24}
$$

with a sufficiently large number *ν*. This dynamics is similar to that of an overdamped mass-spring-damper system, where the spring support moves according to *ui*(*t*). If *<sup>ν</sup>* <sup>→</sup> <sup>∞</sup> then *Vi*|*<sup>t</sup>* <sup>→</sup> *ui*(*t*) for *<sup>t</sup>* <sup>→</sup> <sup>∞</sup>, independent of the initial value *Vi*|0. That is, for *<sup>ν</sup>* <sup>→</sup> <sup>∞</sup>, the control vector function *<sup>u</sup>*(*t*) represents the limit cycle *<sup>V</sup>*∞(*t*) of the volumes.

#### **5. Results**

The previously described optimization algorithm was applied to the exemplary alpha-Stirling engine model presented in Section 2 to optimize the control (piston paths) regarding both power and efficiency. The optimizations were performed for varying cycle time *τ*, or correspondingly varying engine speed 1/*τ*. With help of the penalty function from Equation (16), the working volumes were restricted to an approximate range from 100 cm<sup>3</sup> to 1100 cm<sup>3</sup> . The initial control, from which the optimizations were started, corresponds to harmonic piston paths exploiting the complete admissible working volume range and featuring a phase shift of 90◦ . This harmonic control will, in the following, also be used as a benchmark.

In Figure 2, the Stirling engine's power-optimal working volume trajectories are shown for different engine speeds.

It can be seen that the power-optimal volume trajectories significantly deviate from harmonic shapes. At low engine speeds, it is favorable to let the pistons rest in their extreme positions for about half the cycle time. This feature is less pronounced at higher engine speeds. Above ca. 1100 rpm and ca. 1300 rpm the power-optimal working volume trajectories detach from the maximum volume bounds, as can be seen in Figure 2 at *t*/*τ* ≈ 1/4 and *t*/*τ* ≈ 1, respectively. At the highest considered engine speed of 2000 rpm, the swept volumes are considerably reduced. The power-optimal engine speed is ca. 920 rpm, which is highlighted in Figure 2 by thick black lines. It is interesting that the power-optimal piston motions found in [22] by direct optimization are similar to the trajectories from Figure 2 for low and medium engine speeds. This is the case even though, in the current study, a Stirling engine with different parameters, and a much more detailed regenerator model, are considered. Nevertheless, compared to the parametric AS class of piston motions used in [22] for direct optimization, the trajectories from Figure 2 feature more details.

In Figure 3, the Stirling engine's efficiency-optimal working volume trajectories are shown, again for different engine speeds.

**Figure 2.** Power-optimal trajectories of the working volumes *V*<sup>H</sup> (solid lines) and *V*<sup>L</sup> (dotted lines) of an exemplary alpha-Stirling engine for varying engine speed 1/*τ*. The power-optimal engine speed is ca. 920 rpm (thick black lines).

**Figure 3.** Efficiency-optimal trajectories of the working volumes *V*<sup>H</sup> (solid lines) and *V*<sup>L</sup> (dotted lines) of an exemplary alpha-Stirling engine for varying engine speed 1/*τ*. The efficiency-optimal engine speed is ca. 430 rpm (thick black lines). The power-optimal engine speed for the efficiency-optimal control is ca. 670 rpm (thick gray lines).

Obviously, the efficiency-optimal volume trajectories are different from the poweroptimal ones. The tendency to let the pistons rest in their extreme positions at low engine speeds can also be observed here. However, especially for the cold working volume, this feature is less distinctive than with the power-optimal trajectories. The efficiencyoptimal volume trajectories tend to involve smaller piston velocities and reduced swept volumes, when compared to the power-optimal trajectories for the same engine speed. The efficiency-optimal engine speed is ca. 430 rpm, which is highlighted in Figure 3 by thick black lines.

In Figure 4 the Stirling engine's power and efficiency are plotted against the engine speed for both the power-optimal control (solid, blue) and the efficiency-optimal control (dashed, green). Moreover, the values obtained with the harmonic control (dotted, grey) from which the optimizations were started is shown as a benchmark.

**Figure 4.** Net output power *P* and efficiency *η* for the harmonic (gray, dotted), power-optimal (blue, solid), and efficiency-optimal (green, dasehd) control plotted against the engine speed.

At the maximum power point, the relative power gain due to the power-optimal control is ca. 28%. The relative efficiency gain at the maximum efficiency point, which is achieved with the efficiency-optimal control, is ca. 10%.

Obviously, for a given engine speed, the power-optimal and efficiency-optimal controls each lead to maximum power and maximum efficiency, respectively. However, the optimization for one performance measure may come at the cost of the other. This can, in particular, be seen in the left-hand subfigure for the efficiency-optimal control: For medium engine speeds, the optimization for maximum efficiency leads to remarkable reductions in the output power. There, the power even drops below that obtained with the harmonic control. As opposed to that, the power-optimal control also leads to higher efficiency than the harmonic control over the complete considered range of engine speeds. Note, however, that this may be different for a different set of model parameters.

In the left-hand subfigure of Figure 4, it can be seen that the efficiency-optimal control (green, dashed) leads to a sharp bend at the maximum power point. This occurs at an engine speed of ca. 670 rpm, where the hot working volume trajectory (solid gray line in Figure 3) starts to detach from the maximum volume bound at *t*/*τ* ≈ 1.

Note that optimal controls obtained with the optimization algorithm introduced above generally constitute local optima. This becomes obvious in Figure 4 when, for example, comparing the powers resulting from the power-optimal controls at *τ* = 0.06 s (1000 rpm) and *τ* = 0.12 s (500 rpm). At *τ* = 0.06 s the power is close to the maximum value, whereas at *τ* = 0.12 s the power is much lower. However, it is clear that the power-optimal control obtained for *τ* = 0.06 s is an admissible periodic control for *τ* = 0.12 s as well, since two

oscillations could be performed in one cycle. Hence, the power obtained for *τ* = 0.06 s (1000 rpm) constitutes a lower bound for the global optimum at *τ* = 0.12 s (500 rpm). It can be concluded that the optimal control found for *τ* = 0.12 s (500 rpm) is a local optimum and not the global one. Similar arguments are applicable for the efficiency-optimal control, for example, with the efficiency values obtained at 400 rpm and 200 rpm.

#### **6. Summary**

In this paper, we presented an iterative gradient method for optimizing the piston paths of Stirling engines, which is based on Cyclic Optimal Control Theory. After a brief introduction to Stirling engines, we outlined an exemplary irreversible alpha-Stirling engine model. This model served for illustrative purposes and is structured in such a way that its state dynamics can be expressed as a set of coupled ordinary differential equations. The Stirling engine's state dynamics is influenced by an explicitly time-dependent, vectorvalued control function. In particular, this control function determines the temporal paths of the engine's two working pistons.

The question of which realization of those piston paths—or correspondingly what control function—leads to optimal engine performance regarding a specific objective for a fixed cycle time, constitutes a cyclic optimal control problem. We introduced the necessary conditions of optimality for this problem, involving the definition of a Hamilton function, a path target function, and a set of adjoint ordinary differential equations (costate dynamics). For the optimization of Stirling engines, we presented proper definitions of path target functions to maximize either power or efficiency. Maximum and minimum volume constraints were accounted for here via an additional penalty term in the respective path target function.

Then, we gave a detailed description of an iterative optimization algorithm that starts with a predefined periodic control function, which is gradually shifted over the course of the iterations, so as to gradually enlarge the objective and approach the optimal control. To determine those gradual shifts of the control in every iteration, not only the engine's state dynamics needs to be solved, but also the costate dynamics. The cyclic optimal control problem features a symmetry that manifests itself in attractive and repulsive limit cycles in the state and costate dynamics, respectively. The algorithm exploits these limit cycles to solve the state and costate dynamics for periodic boundary conditions.

This optimization algorithm was applied to both, power-optimize and efficiencyoptimize the control (piston paths) of the above-mentioned exemplary alpha-Stirling engine model for a range of cycle times. The optimization results were compared to piston paths that correspond to a harmonic control with 90◦ phase shift, which had been used as an initial control for the optimizations. At the maximum power point, the relative power gain due to the power-optimal control is ca. 28%, whereas the relative efficiency gain due to the efficiency-optimal control at the maximum efficiency point is ca. 10%.

The developed optimization algorithm can easily be applied to other Stirling engine models that involve additional irreversibilities or transfer laws that are adapted to describe a specific experimental engine setup. Design limitations such as maximum pressure or maximum piston acceleration can be included in the optimal control problem in terms of additional constraints. Moreover, other types of machines, such as beta-type or free-piston Stirling engines, Stirling cryocoolers, or Vuilleumier refrigerators can be optimized analogously.

#### **7. Conclusions**

For the cyclic optimal control problems of a Stirling engine investigated here, we observed that the existence of an attractive limit cycle in the state dynamics translates to a repulsive limit cycle in the costate dynamics. We conclude that a very stable and numerically efficient way to solve such problems is through using indirect iterative gradient algorithms that exploit these limit cycles to solve the state and costate dynamics by temporal forward and backward integration, respectively.

In the case of the considered exemplary Stirling engine, the conducted piston path optimizations lead to significant performance gains of ca. 28% in maximum power and ca. 10% in maximum efficiency. This generally confirms earlier results [22–24] and we conclude that it can be very worthwhile to perform such optimizations during the design process of new engines or to improve the control strategy of existing ones. Numerous types of energy conversion systems operate cyclically. We expect that cyclic optimal control theory and, in particular, the developed algorithm can be applied to a wide class of these systems, in order to identify potential performance improvements and enhance their sustainability.

**Author Contributions:** Investigation, R.P. and K.H.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the German Federal Ministry of Education and Research grant number FKZ01LY1706B.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Review* **Recent Advances on Boundary Conditions for Equations in Nonequilibrium Thermodynamics**

**Wen-An Yong 1,2,\* and Yizhou Zhou <sup>3</sup>**


**\*** Correspondence: wayong@tsinghua.edu.cn

**Abstract:** This paper is concerned with modeling nonequilibrium phenomena in spatial domains with boundaries. The resultant models consist of hyperbolic systems of first-order partial differential equations with boundary conditions (BCs). Taking a linearized moment closure system as an example, we show that the structural stability condition and the uniform Kreiss condition do not automatically guarantee the compatibility of the models with the corresponding classical models. This motivated the generalized Kreiss condition (GKC)—a strengthened version of the uniform Kreiss condition. Under the GKC and the structural stability condition, we show how to derive the reduced BCs for the equilibrium systems as the classical models. For linearized problems, the validity of the reduced BCs can be rigorously verified. Furthermore, we use a simple example to show how thus far developed theory can be used to construct proper BCs for equations modeling nonequilibrium phenomena in spatial domains with boundaries.

**Keywords:** hyperbolic relaxation system; structural stability condition; generalized Kreiss condition

#### **1. Introduction**

Irreversible thermodynamics is a theory in physics for the mathematical modeling of nonequilibrium processes. The resultant models are usually time-dependent partial differential equations (PDEs) [1]. So far, the theory is still developing, and there are no well-accepted rules in establishing the equations. Therefore a number of different theories exist [2–8], leading to various PDEs. It is a challenging issue to evaluate the reasonableness of the different equations. To do so, four fundamental requirements were proposed and expounded in [9]. They are the observability of physical phenomena, time irreversibility, long-time tendency to equilibrium, and compatibility with possibly existing classical theories.

Recall well-known theories extended irreversible thermodynamics (EIT) [4,5]), rational extended thermodynamics (RET) [6]), and conservation–dissipation formalism (CDF) [8]). Irreversible processes can be described with PDEs of the form:

$$
\mathcal{U}I\_t + \sum\_{j=1}^d F\_j(\mathcal{U}I)\_{x\_j} = \tilde{\mathcal{Q}}(\mathcal{U}).\tag{1}
$$

Here, *<sup>U</sup>* <sup>=</sup> *<sup>U</sup>*(*x*, *<sup>t</sup>*) <sup>∈</sup> *<sup>G</sup>* <sup>⊂</sup> <sup>R</sup>*<sup>n</sup>* is the unknown function of *<sup>x</sup>* = (*x*1, *<sup>x</sup>*2, , *xd*) with *d* (=1, 2 or 3) for the spatial dimension and *t* > 0, *G* is an open set called the state space, *Fj*(*U*) is the flux along the *xj*-direction, source term *<sup>Q</sup>*˜(*U*) represents the thermodynamical force, and subscripts *t*, *xj* denote partial derivatives with respect to the corresponding variables. Very often, the source term can be written as

$$
\bar{\mathbb{Q}}(U) = \begin{pmatrix} 0 \\ \\ q(U) \end{pmatrix}.
$$

**Citation:** Yong, W.-A.; Zhou, Y. Recent Advances on Boundary Conditions for Equations in Nonequilibrium Thermodynamics. *Symmetry* **2021**, *13*, 1710. https:// doi.org/10.3390/sym13091710

Academic Editors: Róbert Kovács, Patrizia Rogolino, Francesco Oliveri, Charles F. Dunklm, Marin Marin and Mariano Torrisi

Received: 10 May 2021 Accepted: 1 September 2021 Published: 16 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Accordingly, we partition

$$\mathcal{U} = \begin{pmatrix} u \\ v \end{pmatrix}, \quad F\_{\vec{\jmath}}(\mathcal{U}) = \begin{pmatrix} f\_{\vec{\jmath}}(u\_{\vec{\nu}}v) \\ g\_{\vec{\jmath}}(u\_{\vec{\nu}}v) \end{pmatrix}.$$

With this partition, *<sup>u</sup>* <sup>∈</sup> <sup>R</sup>*n*−*<sup>r</sup>* stands for the conserved variables, and *<sup>v</sup>* <sup>∈</sup> <sup>R</sup>*<sup>r</sup>* is referred to as nonequilibrium variables. To cover more equations from EIT [4,5] and general equation for nonequilibrium reversible–irreversible coupling (GENERIC) [7,10]), we also consider the following equation.

$$
\mathcal{U}I\_t + \sum\_{j=1}^d A\_j(\mathcal{U})\mathcal{U}\_{x\_j} = \mathcal{Q}(\mathcal{U}), \tag{2}
$$

which is more general than (1). Indeed, (1) can be written as (2) with *Aj*(*U*) = *<sup>∂</sup>Fj <sup>∂</sup><sup>U</sup>* (*U*).

For such first-order PDEs, the first fundamental requirement corresponds to the hyperbolicity of Equation (2) [11,12]. To see the other requirements, thermodynamic force *Q*(*U*) contains a relaxation time in general, and we explicitly write it as *Q*˜(*U*) = *Q*(*U*)/. The requirements are closely related to the limit as the relaxation time tends to zero the so-called zero relaxation limit. Without considering the boundary conditions, the four requirements are fulfilled by (2) if the PDEs satisfy the structural stability condition proposed in [13,14]. See [14,15].

When the nonequilibrium phenomena occur in a spatial domain with boundaries, the PDEs alone cannot completely describe the physical processes, and proper boundary conditions (BCs) are indispensable. However, it is a challenging task to have such BCs, since the physical meaning of the nonequilibrium variables is often unclear (e.g., the higher-order moments in moment-closure systems [6,16]). On the other hand, formulating proper BCs for hyperbolic systems of PDEs is mathematically a tough issue [17]. For example, the number of BCs for (2) should be equal to the number of positive eigenvalues of coefficient matrix *<sup>A</sup>*1(*U*) in case that the boundary is *<sup>x</sup>*<sup>1</sup> = 0; moreover, the BCs should satisfy the uniform Kreiss condition (UKC) [18] to ensure the well-posedness of the complete model (PDEs together with BCs)—the first fundamental requirement. Furthermore, the rest requirements also apply to the complete model. This suggests to study the zero relaxation limit of initial-boundary-value problems (IBVPs) for hyperbolic systems (1) or (2). The resultant results are expected to be useful in providing reasonable constraints for the possible BCs and in developing a systematical method to construct the required BCs.

The paper presents the first author's considerations in the past two decades and recent results on these issues. It is organized as follows. Section 2 contains some basic knowledges about relaxation system (3) and boundary conditions. In Section 3, we show that the structural stability condition and the UKC do not automatically guarantee compatibility, and present a strengthened version of the UKC—the generalized Kreiss condition (GKC). Section 4 is devoted to the derivation of the boundary condition for the zero relaxation limit under the GKC and the structural stability condition. In Section 5, we use an example to show how the theory presented above can be used to construct BCs for the PDEs (3). Some concluding remarks are given in Section 6.

#### **2. Preliminaries**

To address the long-term tendency and compatibility requirements, we explicitly introduce relaxation time and write Equation (3) as

$$
\hbar L\_l + \sum\_{j=1}^d A\_j(\mathcal{U}) \mathcal{U}\_{x\_j} = \frac{1}{\epsilon} \mathcal{Q}(\mathcal{U}).\tag{3}
$$

For such a small parameter problem, the most fundamental is the following structural stability condition proposed in [13,14]:

(i) There is an invertible *<sup>n</sup>* <sup>×</sup> *<sup>n</sup>*-matrix *<sup>P</sup>*(*U*) and an invertible *<sup>r</sup>* <sup>×</sup> *<sup>r</sup>*-matrix *<sup>S</sup>*(*U*) (<sup>0</sup> <sup>&</sup>lt; *<sup>r</sup>* <sup>≤</sup> *<sup>n</sup>*), defined on equilibrium manifold <sup>E</sup> <sup>=</sup> {*<sup>U</sup>* <sup>∈</sup> *<sup>G</sup>* : *<sup>Q</sup>*(*U*) = <sup>0</sup>}, such that

$$P(\mathcal{U})Q\_{\mathcal{U}}(\mathcal{U}) = \begin{pmatrix} 0 & 0\\ 0 & S(\mathcal{U}) \end{pmatrix} P(\mathcal{U}) \quad \text{for } \mathcal{U} \in \mathcal{E};$$

(ii) As a hyperbolic system, (1) is symmetrizable, that is, there is a positive definite symmetric matrix *<sup>A</sup>*0(*U*), such that

$$A\_0(\mathcal{U})A\_j(\mathcal{U}) = A\_j^\*(\mathcal{U})A\_0(\mathcal{U}) \quad \text{for } \mathcal{U} \in G;$$

(iii) The hyperbolic part and the source term are coupled in the following sense.

$$A\_0(\iota lI)Q\_{lI}(\iota lI) + Q\_{lI}^\*(\iota lI)A\_0(\iota lI) \le -P^\*(\iota lI)\begin{pmatrix} 0 & 0\\ 0 & I\_r \end{pmatrix} P(\iota lI) \quad \text{for } \iota lI \in \mathcal{E}\text{-}I$$

Here and below, *QU* denotes the Jacobian matrix of *<sup>Q</sup>* = *<sup>Q</sup>*(*U*), *Ik* denotes the unit matrix of order *k*, and superscript ∗ denotes the (conjugate) transpose of a matrix.

In the structural stability condition, (i) is just the usual assumption in the corresponding theory for ordinary differential equations [14], that is, only for the source term; (ii) is for the part of spatial derivatives and means that (3) is symmetrizable hyperbolic; and (iii) is a coupling condition that involves the both parts. As observed in [14,19], the structural stability condition is fulfilled by many classical models from mathematical physics. It is related to the Onsager reciprocal relation and can be viewed as a stability criterion for nonequilibrium thermodynamics. Regarding the Onsager relation, matrix *<sup>A</sup>*0(*U*)*QU*(*U*) is often symmetric. For simplicity, this symmetry is assumed throughout the paper.

Under the structural stability condition, it was proved in [13,14] that the solution to the initial-value problem of (3) converges to that of the corresponding equilibrium system as goes to zero. For the case where

$$Q(U) = \begin{pmatrix} 0 \\ \\ q(U) \end{pmatrix}$$

with *<sup>q</sup>* <sup>=</sup> *<sup>q</sup>*(*U*) <sup>∈</sup> <sup>R</sup>*<sup>r</sup>* , we write

$$\iota \mathcal{U} = \begin{pmatrix} u \\ v \end{pmatrix}, \qquad A\_{\dot{j}}(\mathcal{U}) = \begin{pmatrix} A\_{\dot{j}11} & A\_{\dot{j}12} \\ A\_{\dot{j}21} & A\_{\dot{j}22} \end{pmatrix} (\mathcal{U}).$$

In this case, it is reasonable to assume that *q*(*U*) = 0 if and only if *v* = *h*(*u*). Then, the equilibrium system can be expressed as

$$v = h(u),$$

$$u\_t + \sum\_{j=1}^d A\_{j11}(u, h(u))u\_{x\_j} + \sum\_{j=1}^d A\_{j12}(u, h(u))h(u)\_{x\_j} = 0.$$

Under the structural stability condition, it was also proved in [13,20] that the equilibrium system is symmetrizable hyperbolic.

When Equation (3) is given in a spatial domain with boundaries, proper boundary conditions (BCs) should be prescribed to solve the equations. By referring to [21], it is adequate to consider that the spatial domain is half-space

$$(\mathfrak{x}\_1, \mathfrak{x}\_2, \dots, \mathfrak{x}\_d) \equiv (\mathfrak{x}\_1, \mathfrak{x}) \in [0, +\infty) \times \mathbb{R}^{d-1}.$$

For this domain, assume that coefficient matrix *<sup>A</sup>*1(*U*) has *<sup>p</sup>* positive eigenvalues. According to the classical theory of hyperbolic equations [17,22], *p* independent BCs

$$B(\mathcal{U}(0, \mathfrak{x}, t); t) = 0 \tag{4}$$

should be prescribed at boundary *<sup>x</sup>*<sup>1</sup> = 0. Moreover, they should satisfy the uniform Kreiss condition (UKC) [18] to ensure the well-posedness, namely, there is a positive constant *cK*, such that

$$\left| \det \left\{ B\_{l\varPi} \boldsymbol{R}\_{\varPi}^{S} (\boldsymbol{\xi}, \boldsymbol{\omega}) \right\} \right| \geq c\_{K} \sqrt{\det \left\{ \boldsymbol{R}\_{\vartilde{M}}^{S\*} (\boldsymbol{\xi}, \boldsymbol{\omega}) \boldsymbol{R}\_{\vartilde{M}}^{S} (\boldsymbol{\xi}, \boldsymbol{\omega}) \right\} } $$

for all *<sup>ω</sup>* <sup>∈</sup> <sup>R</sup>*d*−<sup>1</sup> and all complex number *<sup>ξ</sup>* with *Re<sup>ξ</sup>* <sup>&</sup>gt; 0. To decode this UKC, we assume that *<sup>A</sup>*1(*U*) is invertible, and know from [18] that matrix

$$
\tilde{M}(\xi,\omega) := A\_1^{-1} \left( -\xi I + i \sum\_{j=2}^d \omega\_j A\_j \right),
$$

has *p* stable eigenvalues (eigenvalues with negative real parts) for *Reξ* > 0 and *<sup>ω</sup>* = (*ω*2,..., *<sup>ω</sup>d*) <sup>∈</sup> <sup>R</sup>*d*−1. Then, there is a full-rank *<sup>n</sup>* <sup>×</sup> *<sup>p</sup>*-matrix *<sup>R</sup><sup>S</sup> <sup>M</sup>*˜ (*ξ*, *<sup>ω</sup>*) (so-called rightstable matrix of *M*˜ [23]) satisfying

$$
\tilde{M}(\lnot\omega) \mathcal{R}^S\_{\mathcal{M}}(\lnot\omega) = \mathcal{R}^S\_{\mathcal{M}}(\lnot\omega) \tilde{\Lambda}(\lnot\omega),
$$

with <sup>Λ</sup>˜ (*ξ*, *<sup>ω</sup>*) <sup>a</sup> *<sup>p</sup>* <sup>×</sup> *<sup>p</sup>*-matrix having *<sup>p</sup>* stable eigenvalues.

When *<sup>A</sup>*1(*U*) is not invertible, boundary *<sup>x</sup>*<sup>1</sup> = 0 is referred to as characteristic. In this case, a modified UKC was formulated in [21] as a sufficient and essentially necessary condition for well-posedness. In addition, it was indicated in [21] that the boundary condition should not involve characteristic modes corresponding to the zero eigenvalue.

#### **3. Generalized Kreiss Condition**

This section shows that the structural stability condition together with the UKC do not automatically guarantee a well-behaved zero relaxation limit. To do this, we start with the linearized version of PDEs in (3):

$$
\delta \mathcal{U}\_t + \sum\_{j=1}^d A\_j \mathcal{U}\_{\mathcal{X}\_j} = \frac{1}{\epsilon} \mathcal{Q} \mathcal{U}\_t \tag{5}
$$

together with homogeneous boundary conditions

$$BLI(0, \pounds, t) = 0\tag{6}$$

for *<sup>x</sup>* = (*x*1, *<sup>x</sup>*ˆ)=(*x*1, *<sup>x</sup>*2, ... , *xd*) with *<sup>x</sup>*<sup>1</sup> <sup>&</sup>gt; 0. Here, *Aj*(*<sup>j</sup>* <sup>=</sup> 1, 2, ... , *<sup>d</sup>*) and *<sup>Q</sup>* are *<sup>n</sup>* <sup>×</sup> *<sup>n</sup>* constant matrices and *B* is a full-rank constant *p* × *n*-matrix with *p* the number of positive eigenvalues of *A*1.

Following [18,23,24], we consider the corresponding eigenvalue problem:

$$\begin{cases} \, \, ^\zeta \hat{\mathcal{U}} + A\_1 \hat{\mathcal{U}}\_{x\_1} - i \sum\_{j=2}^d \omega\_j A\_j \hat{\mathcal{U}} = Q \hat{\mathcal{U}}, \\\, \, \, ^B \hat{\mathcal{U}}(0) = 0 \end{cases} \tag{7}$$

for complex number *<sup>ξ</sup>* with *Re<sup>ξ</sup>* <sup>&</sup>gt; 0 and *<sup>ω</sup>* = (*ω*2, ... , *<sup>ω</sup>d*) <sup>∈</sup> <sup>R</sup>*d*−1. Let *<sup>U</sup>*<sup>ˆ</sup> <sup>=</sup> *<sup>U</sup>*<sup>ˆ</sup> (*x*1) be a bounded solution to the above problem. It is not difficult to verify that

$$\mathcal{U}^{\mathfrak{c}}(\mathfrak{x},t) := \exp\left(\frac{\mathfrak{x}t}{\mathfrak{c}} - \frac{i\omega \cdot \mathfrak{x}}{\mathfrak{c}}\right) \mathcal{U}\left(\frac{\mathfrak{x}\_1}{\mathfrak{c}}\right)$$

solves (5) and (6) with a bounded initial value. Since *Reξ* > 0, such a solution exponentially increases as goes to zero for any *t* > 0. Apparently, the zero relaxation limit is not well-behaved or does not exist in this situation.

In order to see when problem (7) has bounded solutions, we assume *A*<sup>1</sup> to be invertible and rewrite (7) as

$$\frac{d\hat{\mathcal{U}}}{d\mathbf{x}\_1} = M(\xi, \omega, \mathbf{1})\hat{\mathcal{U}}\tag{8}$$

with

$$M(\xi, \omega, \eta) = A\_1^{-1} \left(\eta Q - \xi I + i \sum\_{j=2}^d \omega\_j A\_j\right).$$

Here, *η* ≥ 0 was introduced to cover the case without source term *Q*, namely, *M*(*ξ*, *ω*, 0) = *M*˜ (*ξ*, *ω*) in Section 2.

According to Lemma 2.3 in [23], under the structural stability condition, matrix *<sup>M</sup>* <sup>=</sup> *<sup>M</sup>*(*ξ*, *<sup>ω</sup>*, *<sup>η</sup>*) has precisely *<sup>p</sup>* stable eigenvalues (with negative real parts) and (*<sup>n</sup>* <sup>−</sup> *<sup>p</sup>*) unstable eigenvalues. Then, there is a full-rank *<sup>n</sup>* <sup>×</sup> *<sup>p</sup>*-matrix *<sup>R</sup><sup>S</sup> <sup>M</sup>* <sup>=</sup> *<sup>R</sup><sup>S</sup> <sup>M</sup>*(*ξ*, *<sup>ω</sup>*, *<sup>η</sup>*) satisfying

$$MR\_M^S = R\_M^S \Lambda$$

with <sup>Λ</sup> <sup>a</sup> *<sup>p</sup>* <sup>×</sup> *<sup>p</sup>*-matrix having *<sup>p</sup>* stable eigenvalues. If the *<sup>p</sup>* <sup>×</sup> *<sup>p</sup>*-matrix *BR<sup>S</sup> <sup>M</sup>*(*ξ*, *<sup>ω</sup>*, 1) is not invertible, then there is a nonzero *p*-vector *ζ*, such that *BR<sup>S</sup> <sup>M</sup>*(*ξ*, *<sup>ω</sup>*, 1)*<sup>ζ</sup>* <sup>=</sup> 0. With *<sup>U</sup>*<sup>ˆ</sup> (0) = *RS <sup>M</sup>*(*ξ*, *<sup>ω</sup>*, 1)*<sup>ζ</sup>* as initial value, the ordinary differential Equation (8) has a bounded solution for *x*<sup>1</sup> > 0.

Consequently, we obtain

**Proposition 1.** *If there exists <sup>ξ</sup> with Re<sup>ξ</sup>* <sup>&</sup>gt; <sup>0</sup> *such that the <sup>p</sup>* <sup>×</sup> *p-matrix BR<sup>S</sup> <sup>M</sup>*(*ξ*, *<sup>ω</sup>*, 1) *is singular, then the problem in* (5) *and* (6)*, with a bounded initial value, admits an exponentially increasing solution for t* > 0 *as goes to zero.*

The above discussion can be illustrated with the one-dimensional linearized moment closure system [16,25]:

$$\begin{cases} \, \mathcal{U}l + A\_1 \mathcal{U}\_x = \mathcal{Q}lI / \mathfrak{e}\_\prime \\\\ \, \mathcal{B}lI(0,t) = 0 \end{cases} \tag{9}$$

with *<sup>Q</sup>* <sup>=</sup> diag(0, 0, 0, <sup>−</sup>1) and

$$A\_1 = \begin{pmatrix} 0 & \rho\_0 & 0 & 0\\ \theta\_0/\rho\_0 & 0 & 1 & 0\\ 0 & 2\theta\_0 & 0 & 6/\rho\_0\\ 0 & 0 & \rho\_0\theta\_0/2 & 0 \end{pmatrix}, \quad B = \begin{pmatrix} 0 & 1 & 0 & 0\\ 0 & 0 & \chi\rho\_0\sqrt{\theta\_0} & -1 \end{pmatrix}.$$

Here, *U* = (*ρ*, *u*, *θ*, *f*) is the unknown. *ρ*, *u*, *θ* are fluid density, velocity, and temperature, respectively, while *f* is related to the heat flux. Constants *ρ*<sup>0</sup> and *θ*<sup>0</sup> in coefficient matrix *A*<sup>1</sup> denote density and temperature at the equilibrium. Constant *χ* in the BC is a positive parameter. Obviously, the moment-closure system satisfies the structural stability condition with symmetrizer

$$A\_0 = \begin{pmatrix} 2\theta\_0^3 & 0 & 0 & 0\\ 0 & 2\rho\_0^2 \theta\_0^2 & 0 & 0\\ 0 & 0 & \rho\_0^2 \theta\_0 & 0\\ 0 & 0 & 0 & 12 \end{pmatrix}.$$

It was shown in [25] that the BC in (9) satisfies the UKC if and only if

$$\chi \neq \frac{\sqrt{3+\sqrt{6}} + \sqrt{3-\sqrt{6}}}{2(3+\sqrt{3})}.$$

Moreover, Proposition 3.3 in [25] says that if

$$\chi \in (0, \frac{\sqrt{3+\sqrt{6}} + \sqrt{3-\sqrt{6}}}{2(3+\sqrt{3})}),$$

there is a positive *ξ* such that *BR<sup>S</sup> <sup>M</sup>*(1, *<sup>ξ</sup>*) is singular. According to Proposition 1, the structural stability condition and the UKC do not automatically imply the existence of the zero relaxation limit.

Proposition 1 and the example above indicate the following necessary condition for the existence of zero relaxation limit:

$$\det\{BR\_M^S(\mathfrak{J},\omega,1)\} \neq 0,\quad \text{for any}\quad \text{Re}\mathfrak{J} > 0,\omega \in \mathbb{R}^{d-1}.$$

Motivated by this, the following condition was proposed in [23]: **Generalized Kreiss Condition.** There exists a constant *cK* > 0, such that

$$\left| \det \left\{ BR\_M^S(\mathfrak{J}, \omega, \eta) \right\} \right| \geq c\_K \sqrt{\det \{ R\_M^{S\*}(\mathfrak{J}, \omega, \eta) R\_M^S(\mathfrak{J}, \omega, \eta) \} }$$

for all *<sup>η</sup>* <sup>≥</sup> 0, *<sup>ω</sup>* <sup>∈</sup> <sup>R</sup>*d*−<sup>1</sup> and *<sup>ξ</sup>* with *Re<sup>ξ</sup>* <sup>&</sup>gt; 0.

Clearly, this GKC does not depend on the special choice of *R<sup>S</sup> <sup>M</sup>*(*ξ*, *<sup>ω</sup>*, *<sup>η</sup>*). It recovers the standard UKC when *η* = 0. Unlike the UKC, the GKC involves parameters even for one-dimensional problems. Moreover, it was shown in [23] that the GKC holds if the BCs are strictly dissipative in the following sense:

**Definition 1.** *The boundary condition is referred to as strictly dissipative with respect to symmetrizer A*<sup>0</sup> *and A*1*, if there is a positive constant cd, such that*

$$y^\* A\_0 A\_1 y \le -c\_d |y|^2 + c\_d^{-1} |By|^2$$

*for all y* <sup>∈</sup> <sup>R</sup>*n.*

We conclude this section with the following remark:

**Remark 1.** *In studying IBVPs for hyperbolic systems, it is important to distinguish whether the boundary is characteristic or not [21], that is, whether A*<sup>1</sup> *is invertible or not. In the above discussions, the invertibility of A*<sup>1</sup> *plays the role of the starting point. When A*<sup>1</sup> *is not invertible, a modified GKC was proposed in [26]. It covers the above GKC for noncharacteristic problems.*

#### **4. Reduced Boundary Conditions**

Besides the GKC, deriving reduced boundary conditions is another crucial issue in studying the zero relaxation limit for the initial-boundary-value problems. The relaxation limit satisfies the equilibrium system under the structural stability condition [14]. When the problem is given in the half-space, the equilibrium system alone can not determine the limit. It must be supplemented with proper boundary conditions. As part of the system determining the limit, such BCs should be completely derived from the relaxation system and its BCs. BCs thus derived are called reduced boundary conditions.

For linearized system (3) satisfying the structural stability condition, we assume without loss of generality that *Q* = diag(0, *S*) with *S* being a negative definite matrix. Corresponding to the partition of *Q*, we write coefficient matrix *Aj* as

$$A\_{\bar{\jmath}} = \begin{pmatrix} A\_{\bar{\jmath}11} & A\_{\bar{\jmath}12} \\ A\_{\bar{\jmath}21} & A\_{\bar{\jmath}22} \end{pmatrix}$$

Thus, the aforementioned equilibrium system reads as

$$\begin{cases} \upsilon = 0, \\\\ u\_t + \sum\_{j=1}^d A\_{j11} u\_{x\_j} = 0. \end{cases} \tag{10}$$

.

Recall that coefficient matrices *Aj*<sup>11</sup> are symmetrizable.

Let *p*<sup>1</sup> be the number of positive eigenvalues for matrix *A*111. According to the classical theory of IBVPs for hyperbolic systems [17,22], equilibrium system (10) needs *p*<sup>1</sup> BCs satisfying the UKC to be well-posed. As mentioned before, these BCs should be completely determined from the relaxation system and its BCs:

$$BLI(0, \pounds, t) = (B\_{\mu \nu} B\_{\upsilon}) \binom{\mu}{\upsilon} (0, \pounds, t) = b(t). \tag{11}$$

In this regard, it was established in [23,27] that

**Theorem 1.** *Under the structural stability condition and the GKC, there exists a p*<sup>1</sup> × *p-matrix Bp, unique up to an invertible p*<sup>1</sup> × *p*1*-matrix multiplying Bp from left, such that relation*

$$B\_p B\_u u(0, \pounds, t) = B\_p b(\pounds, t) \tag{12}$$

*satisfies the UKC as a BC for the equilibrium system. Moreover, if boundary <sup>x</sup>*<sup>1</sup> = <sup>0</sup> *is characteristic for the equilibrium system, BC* (12) *does not involve the characteristic modes corresponding to the zero eigenvalue.*

The equilibrium system with the reduced BC (12) consist of a well-posed IBVP since the reduced BC satisfies the UKC. The proof of this theorem involves the perturbation theory of linear operators [28] and subtle matrix analysis. The key is to analyze the limit of right-stable matrix *R<sup>S</sup> <sup>M</sup>*(*ξ*, *<sup>ω</sup>*, *<sup>η</sup>*) as *<sup>η</sup>* <sup>→</sup> <sup>∞</sup>. According to the GKC, *BR<sup>S</sup> <sup>M</sup>*(*ξ*, *<sup>ω</sup>*, <sup>+</sup>∞) is a *p* × *p* invertible matrix. Then, *Bp* can be chosen as the *p*<sup>1</sup> × *p*-matrix consisting of the first *<sup>p</sup>*<sup>1</sup> rows of [*BR<sup>S</sup> <sup>M</sup>*(*ξ*, *<sup>ω</sup>*, <sup>+</sup>∞)]−1. Such a matrix *Bp* is independent of *<sup>ξ</sup>* and *<sup>ω</sup>*. We omit the details here, and the interested reader is referred to to [23,27].

Regarding this theorem:

#### **Remark 2.**

• *For nonlinear problems, reduced BCs are quite complicated, but can be derived by considering boundary-layer equations*

$$A\_1(\mathcal{U})\mathcal{U}\_y = \mathcal{Q}(\mathcal{U}), \qquad y = \frac{x\_1}{\epsilon} \ge 0.$$

*For the details, the interested reader is referred to [23].*

• *In case the boundary is characteristic (for the relaxation system), a similar result can be found in [26].*

To see the validity of the reduced BCs, we consider asymptotic solutions of the form:

$$
\begin{pmatrix} \mu\_{\varepsilon} \\\\ \upsilon\_{\varepsilon} \end{pmatrix} (\mathbf{x}\_{1}, \mathbf{\hat{x}}, t; \varepsilon) = \begin{pmatrix} \vec{u} \\\\ \vec{v} \end{pmatrix} (\mathbf{x}\_{1}, \mathbf{\hat{x}}, t) + \begin{pmatrix} \mu\_{0} \\\\ \upsilon\_{0} \end{pmatrix} (\frac{\mathbf{x}\_{1}}{\varepsilon}, \mathbf{\hat{x}}, t)
$$

$$
+ \begin{pmatrix} \mu\_{1} \\\\ \upsilon\_{1} \end{pmatrix} (\frac{\mathbf{x}\_{1}}{\sqrt{\varepsilon}}, \mathbf{\hat{x}}, t) + \sqrt{\varepsilon} \begin{pmatrix} \mu\_{2} \\\\ \upsilon\_{2} \end{pmatrix} (\frac{\mathbf{x}\_{1}}{\sqrt{\varepsilon}}, \mathbf{\hat{x}}, t) \tag{13}
$$

for relaxation system (5) and its BC (11). Here, (*u*¯, *<sup>v</sup>*¯) is the outer solution, and (*μi*, *<sup>ν</sup>i*) (*i* = 0, 1, 2) are boundary-layer corrections satisfying matching conditions

$$(\mu\_{i\prime}\nu\_i)(\lhd\_\prime\pounds, t) = 0, \qquad i = 0, 1, 2.$$

By a standard procedure [27], the outer solution solves equilibrium system (10) with reduced BC (12), while the other terms in (13) can also be determined.

If the boundary is not characteristic for the equilibrium system, there is no boundarylayer with length <sup>√</sup>. In other words, the last two terms in (13) are void, and the three scales in the expansion degenerate to two scales.

The validity of the reduced BC (12) is verified with the following theorem.

**Theorem 2** ([27])**.** *Assuming that the structural stability condition and the GKC hold, the initial data are given at the equilibrium, and the BC is compatible with the initial condition, there exists a constant K* > 0*, such that the exact solution* (*u*, *v*) *fulfils the following estimate.*

$$\left\| \left( \mu^{\mathfrak{c}} - \mu\_{\mathfrak{c}\prime} \upsilon^{\mathfrak{c}} - \upsilon\_{\mathfrak{c}} \right) (\cdot, \cdot, t) \right\|\_{L^{2}(\mathbb{R}^{+} \times \mathbb{R}^{d-1})} \leq K \mathfrak{e}^{1/2} $$

*for all time t* <sup>∈</sup> [0, *<sup>T</sup>*]*.*

Next, we return to one-dimensional linearized moment closure system (9) to illustrate the reduced BCs. For this model, we have

$$A\_1 = \begin{pmatrix} A\_{11} & A\_{12} \\ A\_{21} & A\_{22} \end{pmatrix} = \begin{pmatrix} 0 & \rho\_0 & 0 \\ \theta\_0/\rho\_0 & 0 & 1 \\ 0 & 2\theta\_0 & 0 \\ \hline 0 & 0 & \rho\_0\theta\_0/2 & 0 \end{pmatrix}^\prime$$

and

$$B = \begin{pmatrix} B\_{\mu\_1} B\_{\upsilon} \end{pmatrix} = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & \chi \rho\_0 \sqrt{\theta\_0} \end{pmatrix} \begin{pmatrix} 0 \\ -1 \end{pmatrix}, \qquad b(t) = 0.$$

As shown in [25], coefficient matrix *A*<sup>1</sup> has two positive eigenvalues, 3 <sup>3</sup> <sup>±</sup> <sup>√</sup><sup>6</sup> <sup>√</sup>*θ*0, and *<sup>A</sup>*<sup>11</sup> has one positive eigenvalue <sup>√</sup>3*θ*0, that is, *<sup>p</sup>* <sup>=</sup> 2 and *<sup>p</sup>*<sup>1</sup> <sup>=</sup> 1. Thus equilibrium system

$$
\begin{pmatrix} \bar{\rho} \\ \bar{u} \\ \bar{\theta} \end{pmatrix}\_t + \begin{pmatrix} 0 & \rho\_0 & 0 \\ \theta\_0/\rho\_0 & 0 & 1 \\ 0 & 2\theta\_0 & 0 \end{pmatrix} \begin{pmatrix} \bar{\rho} \\ \bar{u} \\ \bar{\theta} \end{pmatrix}\_{x\_1} = 0
$$

needs one reduced boundary condition. On the other hand, in this case, we have

$$BR\_M^S(\xi, +\infty) = \begin{pmatrix} \sqrt{3\theta\_0}/\rho\_0 & 0\\ 2\chi\theta\_0\sqrt{\theta\_0} & -\chi\theta\_0\sqrt{\theta\_0} \end{pmatrix}.$$

This matrix is obviously invertible, and its inverse is a lower-triangular matrix. Therefore, *Bp* can be chosen as the 1 <sup>×</sup> 2-matrix (1, 0), and the reduced BC is

$$B\_{\overline{\rho}}B\bar{\mathcal{U}}(0,t) = (1,0)\begin{pmatrix} 0 & 1 & 0 & 0\\ 0 & 0 & \chi\rho\_0\sqrt{\theta\_0} & -1 \end{pmatrix} \bar{\mathcal{U}}(0,t) = B\_{\overline{\rho}}b(t) = 0$$

with *U*¯ = (*ρ*¯, *u*¯, ¯ *θ*, ¯ *f*). That is, just no-slip velocity BC

$$
\bar{u}(0, t) = 0.
$$

#### **5. Construction of BCs**

In this section, we show how the theory presented above can be used to construct BCs for PDEs (2). Usually, the relaxation system is derived by augmenting certain classical equations to take account into more refined physical processes. However, proper BCs for the relaxation systems are not always available, since the physical meaning of the nonequilibrium variables is often unclear (e.g., the higher-order moments in momentclosure systems [6,16]). Because of this, no physical means can be expected to obtain the BCs for certain relaxation systems.

For our construction of the BCs, we assume that both the relaxation systems and the well-posed BCs for the corresponding equilibrium systems are given. This assumption is reasonable because there already exist many well-known approaches to construct the relaxation systems, and much knowledge on equilibrium systems and their BCs exists.

We illustrate the construction with the following simple model.

$$\begin{aligned} \partial\_t u - \partial\_x v &= 0, \\ \partial\_t v - \partial\_x \sigma &= 0, \\ \partial\_t (\sigma - Eu) &= (g(u) - \sigma) / \epsilon. \end{aligned} \tag{14}$$

This model was proposed in [29] for the isothermal motions of a viscoelastic material, and can be viewed as an augmentation of classical equations of isothermal elastodynamics:

$$
\partial\_t \mu - \partial\_x \upsilon = 0,
$$

$$
\partial\_t \upsilon - \partial\_x \mathfrak{g}(\mu) = 0.\tag{15}
$$

Here, *u* and *σ* denote strain and stress, *v* is related to particle velocity, *E* is a positive constant called dynamic Young's modulus, is the relaxation time, and *g*(*u*) is a given function of *u*. For this model, the structural stability condition reads as

$$0 < \lg'(u) < E.$$

For simplicity, we assume that the model was defined in the domain where *x* ≥ −*αt* and *g*(*u*) = *Gu*, with *G* > 0 a constant and *α* a parameter. Under change in variables (*x*, *<sup>t</sup>*) <sup>→</sup> (*<sup>x</sup>* <sup>+</sup> *<sup>α</sup>t*, *<sup>t</sup>*), the system (14) becomes

$$
\begin{pmatrix} u \\ v \\ p \end{pmatrix}\_t + \begin{pmatrix} a & -1 & 0 \\ -G & a & -1 \\ 0 & G - E & a \end{pmatrix} \begin{pmatrix} u \\ v \\ p \end{pmatrix}\_x = \frac{1}{\epsilon} \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix} \begin{pmatrix} u \\ v \\ p \end{pmatrix} \tag{16}
$$

with *<sup>p</sup>* <sup>=</sup> *<sup>σ</sup>* <sup>−</sup> *Gu* and (15) reads as

$$
\left(\begin{array}{cc} u \\ v \end{array}\right)\_t + \left(\begin{array}{cc} a & -1 \\ -G & a \end{array}\right) \left(\begin{array}{cc} u \\ v \end{array}\right)\_x = 0. \tag{17}
$$

Moreover, the given BC for the equilibrium system is

$$
\hat{B}\begin{pmatrix} u \\ v \end{pmatrix}(0,t) = \hat{b}(t)
$$

with *B*ˆ = (*B*ˆ 1, *B*ˆ <sup>2</sup>) being a constant full-rank matrix. The number of rows of *<sup>B</sup>*<sup>ˆ</sup> is equal to the number of positive eigenvalues of coefficient matrix

$$A\_{11} \equiv \left( \begin{array}{cc} \alpha & -1 \\ -G & \alpha \end{array} \right).$$

This matrix has eigenvalues *<sup>α</sup>* <sup>±</sup> <sup>√</sup>*G*, while coefficient matrix

$$A\_1 \equiv \begin{pmatrix} \alpha & -1 & 0 \\ -G & \alpha & -1 \\ 0 & G - E & \alpha \end{pmatrix}$$

has eigenvalues *<sup>α</sup>* <sup>±</sup> <sup>√</sup>*<sup>E</sup>* and *<sup>α</sup>*.

To be precise and simple, we chose *<sup>α</sup>* <sup>∈</sup> (<sup>−</sup> <sup>√</sup>*E*, <sup>−</sup> <sup>√</sup>*G*), while other choices can be found in [30]. For such an *α*, boundary *x* = 0 is noncharacteristic for both (17) and (16), boundary matrix *B*ˆ is void, and we need to construct one BC of the form

$$B\left(\begin{array}{c}\mu(0,t)\\v(0,t)\\p(0,t)\end{array}\right) = b\_{\mathfrak{e}}(t)\tag{18}$$

for the relaxation system (16).

To find such a nonzero vector *<sup>B</sup>* <sup>∈</sup> <sup>R</sup>3, we consider the asymptotic solution of form

$$
\begin{pmatrix} u\_{\varepsilon} \\ v\_{\varepsilon} \\ p\_{\varepsilon} \end{pmatrix} (\mathfrak{x}, t; \mathfrak{e}) = \begin{pmatrix} \mathfrak{a} \\ \mathfrak{v} \\ \mathfrak{p} \end{pmatrix} (\mathfrak{x}, t) + \begin{pmatrix} \mu \\ \nu \\ \pi \end{pmatrix} (\frac{\mathfrak{x}}{\mathfrak{e}}, t),
$$

and require it to satisfy the BC (18):

$$B\begin{pmatrix}\vec{u}+\mu\\\vec{v}+\nu\\\vec{p}+\pi\end{pmatrix}(0,t)=b\_0(t).$$

In this way, we obtain some algebraic relations or restrictions on *B*, from which we deduce that

$$B^\* = \begin{pmatrix} 1 + \frac{\mathcal{C}\_1}{\mathcal{a}^2 - \mathcal{G}} \\\\ \sqrt{\mathcal{G}} + \frac{\mathcal{a}\mathcal{C}\_1}{\mathcal{a}^2 - \mathcal{G}} \\\\ \mathcal{C}\_1 \end{pmatrix} \times \begin{pmatrix} 1 + \frac{\mathcal{C}\_2}{\mathcal{a}^2 - \mathcal{G}} \\\\ -\sqrt{\mathcal{G}} + \frac{\mathcal{a}\mathcal{C}\_2}{\mathcal{a}^2 - \mathcal{G}} \\\\ \mathcal{C}\_2 \end{pmatrix}.$$

with *C*1, *C*<sup>2</sup> two free parameters. BCs thus constructed are not unique, which is expected.

It was shown in [30] that, when *C*1, *C*<sup>2</sup> are both sufficiently small, the constructed BCs are strictly dissipative. Thus, they satisfy the GKC. On this basis, the validity of the constructed BCs was rigorously proved in [30].

#### **6. Concluding Remarks**

In this paper, we presented a systematical review on boundary conditions (BCs) for partial differential equations (PDEs) from nonequilibrium thermodynamics. From a stabil-

ity point of view, such PDEs should satisfy the structural stability condition. In particular, they constitute hyperbolic systems, for which the uniform Kreiss condition (UKC) is a sufficient and essentially necessary condition for the well-posedness of the corresponding models (PDEs with BCs). Taking a linearized moment closure system as an example, we showed that the structural stability condition and UKC do not automatically guarantee the compatibility of the models with the corresponding classical models. This motivated the generalized Kreiss condition (GKC)—a strengthened version of the UKC. It is expected that the GKC and the structural stability condition are the very criteria to verify the four fundamental requirements when modeling nonequilibrium phenomena in spatial domains with boundaries.

The compatibility and long-time tendency requirements suggest to study the zero relaxation limit of the PDEs with BCs. It is trivial to obtain the equilibrium system for the limit. However, the equilibrium system alone cannot determine the limit, and proper BCs are needed. The needed BCs are referred to as reduced BCs, and they should be completely determined by the PDEs with the given BCs. The derivation of the reduced BCs is by no means trivial.

Under the GKC and the structural stability condition, we showed how to derive the reduced BCs for general relaxation models. For linearized problems, the validity of the reduced BCs was rigorously verified. For nonlinear problems, such a verification is open. Furthermore, we used a simple example to show how the thus far developed theory can be used to construct proper BCs accompanying PDEs modeling nonequilibrium phenomena in spatial domains with boundaries.

**Author Contributions:** Conceptualization, W.-A.Y.; investigation, W.-A.Y., Y.Z.; writing—original draft preparation, Y.Z.; writing—review and editing, W.-A.Y., Y.Z.; project administration, W.-A.Y.; funding acquisition, W.-A.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (grant number 12071246).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Symmetry* Editorial Office E-mail: symmetry@mdpi.com www.mdpi.com/journal/symmetry

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-7763-0