1. Introduction
The behaviour of aerosols (small droplets suspended in vapour/gases) covers a wide range of phenomena of interest to both scientists and engineers (Sazhin Reference Sazhin2014). In the context of airborne disease transmission (Tang et al. Reference Tang, Li, Eames, Chan and Ridgway2006; Xie et al. Reference Xie, Li, Chwang, Ho and Seto2007), large droplets (
$\gtrsim \mathrm {\mu }\textrm {m}$ in radius) fall rapidly to the ground due to gravity, and are therefore transmitted over short distances; smaller droplets (
$\lesssim \mathrm {\mu }\textrm {m}$ in radius), on the other hand, remain suspended in the air for a significant period of time and, hence, are transmitted over large distances. The immediate intuitive questions are ‘How far does a droplet travel?’ and ‘How long does a droplet survive?’ The former is answered by evaluating the momentum of the droplet while the latter by the heat and mass transfer properties of the droplet, with the two being inevitably coupled. At small scales, it is difficult to measure the relevant quantities (e.g. drag forces) on single and multiple droplets; therefore, an accurate modelling is key to understanding transport phenomena (the trajectory and thermal state) during evaporation.
When a spherical particle of radius 
$\hat {R}_0$ moves slowly with a steady translational velocity in a gas, with 
$\hat {R}_0$ being much larger than the mean free path (
$\hat {\lambda }$) of the gas, the flow fields around the particle can be obtained by solving the Stokes equation, and, therefore, one can obtain the drag force acting on the sphere (Lamb Reference Lamb1932). However, the Stokes equation as well as the other classical continuum theories of gas dynamics based on the Euler or Navier–Stokes–Fourier (NSF) equations, which are in general accurate for describing a flow in the so-called hydrodynamic regime (Knudsen number 
${Kn}\lesssim 0.001$), cannot explain the gas flow past a micro-/nano-sized drop (for which 
${Kn}=\hat {\lambda }/\hat {R}_0 \sim 1$) accurately (Cercignani Reference Cercignani2000; Sone Reference Sone2002; Struchtrup Reference Struchtrup2005). Such flows can only be predicted by employing tools of kinetic theory.
Within the framework of kinetic theory, the vapour/gas is described by the Boltzmann equation (Cercignani Reference Cercignani1975; Kremer Reference Kremer2010), which is the evolution equation for the (velocity) distribution function of the gas and accurately describes flows for all Knudsen numbers. Notwithstanding, the exact analytical solution of the Boltzmann equation is intractable for a general flow problem while numerical solutions of the Boltzmann equation demand a very high computational cost in general (particularly, for transition-regime flows). Although, very recently, some general synthetic iterative schemes have been proposed for solving the Boltzmann equation (Su, Zhang & Wu Reference Su, Zhang and Wu2021; Zhu et al. Reference Zhu, Pi, Su, Li, Zhang and Wu2021) that reduce the computational cost of numerical solutions of the Boltzmann equation significantly, it is still appreciable to describe a process with a set of partial differential equations in physical variables (or moments of the distribution function) as they provide a deeper insight to the underlying physics of the process and permit the deployment of classical mathematical techniques. A variety of approximation techniques, aiming to derive sets of partial differential equations that can describe processes at moderate Knudsen numbers, have been developed; see, e.g. Struchtrup (Reference Struchtrup2005) and Torrilhon (Reference Torrilhon2016) and references therein. Among others, the two well-known approximation techniques are the Chapman–Enskog (CE) expansion (Chapman & Cowling Reference Chapman and Cowling1970) and method of moments (Grad Reference Grad1949).
The CE expansion is suitable for near-equilibrium gas flows and relies on an asymptotic expansion of the distribution function around the equilibrium distribution function (which is a Maxwellian) in powers of a small parameter – typically the Knudsen number. The distribution function in the Boltzmann equation is replaced with this expanded distribution function, and coefficients of each power of the Knudsen number on both sides of the resulting equation are compared. The comparison of coefficients of 
${Kn}^0$ yields the Euler equations, of 
${Kn}^1$ results into the NSF equations, of 
${Kn}^2$ leads to the Burnett equations, of 
${Kn}^3$ gives the so-called super-Burnett equations, and so on (Struchtrup Reference Struchtrup2005; Bobylev Reference Bobylev2018). The Burnett (also the super-Burnett and beyond) equations are unstable (Bobylev Reference Bobylev1982), and due to lack of any coherent entropy inequality these equations lead to physically inadmissible solutions (Struchtrup & Nadler Reference Struchtrup and Nadler2020).
Grad – in a seminal work (Grad Reference Grad1949) – proposed an asymptotic solution procedure for the Boltzmann equation through the moment equations. The method is referred to as Grad's moment method and derives the governing equations for the macroscopic quantities, namely mass density, velocity, temperature, stress, heat flux and other higher-order moments, through the Boltzmann equation since these macroscopic quantities are related to the distribution function via its moments. The well-known Grad 13-moment equations, at least in the linearised case, are accompanied with a proper entropy inequality (Rana & Struchtrup Reference Rana and Struchtrup2016) and are stable.
Over the last few decades, significant progress has been made towards improving Grad's approach while keeping their advantageous features (see, e.g. Karlin et al. Reference Karlin, Gorban, Dukek and Nonnenmacher1998; Ohr Reference Ohr2001; Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2003; Struchtrup Reference Struchtrup2004; Cai, Fan & Li Reference Cai, Fan and Li2014; Karlin Reference Karlin2018; Cai & Wang Reference Cai and Wang2020). In particular, Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003) proposed a regularization procedure that derives a closed set of moment equations – referred to as the regularised 13-moment (R13) equations. The R13 equations, since their derivation, have paved the way to several developments, including numerical frameworks (see, e.g. Koellermeier, Schaerer & Torrilhon Reference Koellermeier, Schaerer and Torrilhon2014; Claydon et al. Reference Claydon, Shrestha, Rana, Sprittles and Lockerby2017; Cai & Wang Reference Cai and Wang2020). It can be shown via the order of magnitude approach (Struchtrup Reference Struchtrup2004) that the R13 equations possess the accuracy of order 
${Kn}^3$ whereas the Euler, NSF and Grad 13-moment equations have the accuracy of order 
${Kn}^0$, 
${Kn}^1$ and 
${Kn}^2$, respectively (Struchtrup Reference Struchtrup2004, Reference Struchtrup2005; Torrilhon Reference Torrilhon2016). It is generally accepted that to maintain accuracy, the number of moments ought to be increased with the Knudsen number. Aiming to increase the applicability of the regularization procedure to describe processes in the moderate transition regime, Gu & Emerson (Reference Gu and Emerson2009), by following the approach of Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003), derived the regularised 26-moment (R26) equations. The R26 equations are anticipated to possess the accuracy of order 
${Kn}^5$ in the bulk region (i.e. outside of the Knudsen layer) (Struchtrup Reference Struchtrup2004, Reference Struchtrup2005).
Boundary conditions complementing the Grad or regularised moment equations are typically obtained by employing the Maxwell accommodation model (Maxwell Reference Maxwell1879) along with Grad's closure at the boundary yielding the so-called macroscopic boundary conditions (MBC). The paper by Gu & Emerson (Reference Gu and Emerson2007) may be regarded as the pioneering work on wall boundary conditions for the R13 equations using the Maxwell accommodation model. Nevertheless, some inconsistencies and oscillations were reported in this work. To overcome these shortcomings, Torrilhon & Struchtrup (Reference Torrilhon and Struchtrup2008) proposed a theory of boundary conditions for the R13 equations using the Maxwell accommodation model but taking physical and mathematical requirements for the system into account. Moreover, following the approach of Torrilhon & Struchtrup (Reference Torrilhon and Struchtrup2008), the boundary conditions for the R26 equations – based on the Maxwell accommodation model – were obtained in Gu & Emerson (Reference Gu and Emerson2009). The boundary conditions derived through the Maxwell accommodation model are generically referred to as the Maxwell boundary conditions. Solutions of the R13 and R26 equations along with the Maxwell boundary conditions have been obtained for some benchmark problems (Taheri, Struchtrup & Torrilhon Reference Taheri, Struchtrup and Torrilhon2009; Gu, Emerson & Tang Reference Gu, Emerson and Tang2010; Struchtrup & Taheri Reference Struchtrup and Taheri2011; Gu & Emerson Reference Gu and Emerson2014; Claydon et al. Reference Claydon, Shrestha, Rana, Sprittles and Lockerby2017; Gu et al. Reference Gu, Barber, John and Emerson2019) in order to investigate rarefaction effects in gas at moderate Knudsen numbers. Recently, a hybrid approach combining the discrete velocity method with the MBC in the near-wall region, where the moment method is expected to be less accurate due to strong non-equilibrium effects, has been proposed by Yang et al. (Reference Yang, Gu, Wu, Emerson, Zhang and Tang2020). The approach not only yields more accurate solutions for wall-bounded flows but also reduces the computational cost in numerical simulations.
Thus far, analytical solutions for evaporation and condensation problems have been obtained from the linearised R13 (LR13) (Beckmann et al. Reference Beckmann, Rana, Torrilhon and Struchtrup2018) and linearised R26 (LR26) (Rana, Lockerby & Sprittles Reference Rana, Lockerby and Sprittles2018b) equations and compared with numerical solutions from kinetic theory. The R26 equations can capture Knudsen layers – a region of strong non-equilibrium within a few mean free paths of the boundary – more accurately than the R13 equations giving a good agreement with kinetic theory for the Knudsen number up to 1. The moment equations are now being implemented into state-of-the-art open-source softwares (Torrilhon & Sarna Reference Torrilhon and Sarna2017; Westerkamp & Torrilhon Reference Westerkamp and Torrilhon2019), providing a numerical platform to solve these equations for engineering applications.
Despite the success of the moment equations, there have been concerns over the accuracy of MBC obtained using Grad's closure at the boundary. Rana & Struchtrup (Reference Rana and Struchtrup2016) showed that the MBC for the LR13 equations violate the Onsager reciprocity relations (Onsager Reference Onsager1931a,Reference Onsagerb), and derived a new set of phenomenological boundary conditions (PBC) for them. For PBC, one writes the entropy generation at the boundary and finds the boundary conditions as phenomenological laws that guarantee positivity of the entropy inequality at the boundary, which in turn leads to stable numerical schemes and convergence of moments methods (Ringhofer, Schmeiser & Zwirchmayr Reference Ringhofer, Schmeiser and Zwirchmayr2001; Sarna & Torrilhon Reference Sarna and Torrilhon2018a). The PBC for the LR13 equations are thermodynamically admissible for all processes and comply with the Onsager reciprocity relations (Rana & Struchtrup Reference Rana and Struchtrup2016).
In this paper, we first determine a quadratic form without cross-product terms for the entropy density to prove the 
$H$-theorem (entropy inequality) for the LR26 equations, which in turn leads to the PBC for them. We then solve the LR26 equations along with the associated PBC to investigate a slow flow of a gas/vapour past an evaporating/non-evaporating spherical droplet analytically, and compare the results with those obtained from other continuum theories, existing experimental data and/or numerical results obtained from the linearised Boltzmann equation. Moreover, we estimate the drag force acting on the evaporating/non-evaporating spherical droplet with the LR26 equations. The results not only give an excellent match but – thanks to the analytic solutions – also provide a deeper understanding of the process involved. Through the examples considered in this paper, the present approach evidently provides the best ever results obtained via macroscopic models.
The remainder of the paper is organised as follows. The LR26 equations in the dimensionless form are presented in § 2. The 
$H$-theorem for the LR26 equations is proved in § 3. Physically admissible boundary conditions for the LR26 equations are derived in § 4, where different contributions to the entropy production rate at the interface are also expounded by exploiting the Curie principle and Onsager symmetry relations. Canonical boundary value problems are solved with the LR26 equations with the PBC and MBC in § 5 to assess the Onsager reciprocity relations and entropy generation at the interface. The problem of a steady gas flow past a spherical droplet with and without evaporation is investigated in § 6. Concluding remarks are made in § 7.
2. The linearised R26 equations
The treatment developed in this paper is based entirely on the linearised equations, which result from omitting the nonlinear terms (inertial terms, viscous heating terms, for instance) from the governing equations. This simplification is justified since many micro- and nano-fluid devices involve sufficiently slow flows. For succinctness, we shall present only the LR26 equations (and subsystems), which can be obtained from the original R26 equations derived in Gu & Emerson (Reference Gu and Emerson2009) in a straightforward way. We shall consider small deviations from a constant equilibrium state – given by a constant reference density 
$\hat {\rho }_{0}$, a constant reference temperature 
$\hat {T}_{0}$ and all other fields as zero. Thereby, all the equations can be linearised with respect to the reference state. Moreover, the equations will also be made dimensionless by introducing a length scale 
$\hat {L}_0$, a velocity scale 
$\sqrt {\hat {R} \hat {T}_{0}}$, a time scale 
$\hat {L}_0 / \sqrt {\hat {R} \hat {T}_{0}}$ and the dimensionless deviations in field variables from their respective reference states:
\begin{equation} \left. \begin{gathered}
\rho := \frac{\hat{\rho}-\hat{\rho}_{0}}{\hat{\rho}_{0}},
\quad v_{i} := \frac{\hat{v}_{i}}{\sqrt{\hat{R}
\hat{T}_{0}}}, \quad T :=
\frac{\hat{T}-\hat{T}_{0}}{\hat{T}_{0}},\\ \sigma_{ij} :=
\frac{\hat{\sigma}_{ij}}{\hat{\rho}_{0}\hat{R}
\hat{T}_{0}}, \quad q_{i} :=
\frac{\hat{q}_{i}}{\hat{\rho}_{0} \big(\hat{R}
\hat{T}_{0}\big)^{3/2}},\\ m_{ijk} :=
\frac{\hat{m}_{ijk}}{\hat{\rho}_{0} \big(\hat{R}
\hat{T}_{0}\big)^{3/2}}, \quad R_{ij} :=
\frac{\hat{R}_{ij}}{\hat{\rho}_{0} \big(\hat{R}
\hat{T}_{0}\big)^2}, \quad \varDelta :=
\frac{\hat{\varDelta}}{\hat{\rho}_{0} \big(\hat{R}
\hat{T}_{0}\big)^2}, \\ \varPhi_{ijkl} :=
\frac{\hat{\varPhi}_{ijkl}}{\hat{\rho}_{0} \big(\hat{R}
\hat{T}_{0}\big)^2}, \quad \varPsi_{ijk} :=
\frac{\hat{\varPsi}_{ijk}}{\hat{\rho}_{0} \big(\hat{R}
\hat{T}_{0}\big)^{5/2}}, \quad \varOmega_{i} :=
\frac{\hat{\varOmega}_{i}}{\hat{\rho}_{0} \big(\hat{R}
\hat{T}_{0}\big)^{5/2}}, \end{gathered} \right\}
\end{equation}
where 
$\hat {R}$ is the gas constant; 
$\rho$, 
$v_i$, 
$T$, 
$\sigma _{ij}$ and 
$q_i$ are the dimensionless deviations in the density, macroscopic velocity, temperature, stress tensor and heat flux from their respective reference state; 
$m_{ijk}$, 
$R_{ij}$, 
$\varDelta$, 
$\varPhi _{ijkl}$, 
$\varPsi _{ijk}$ and 
$\varOmega _i$ are dimensionless deviations in the higher-order moments. Note that the notation for the field variables are the same as those in Gu & Emerson (Reference Gu and Emerson2009) except for the quantities with hats here are with dimensions. As a consequence of the linearisation, the dimensionless deviation in the pressure 
$p$ is given by the linearised equation of state, i.e. 
$p = \rho + T$. The non-dimensionalisation also introduces a dimensionless parameter – the Knudsen number
\begin{equation} {Kn}= \frac{\hat{\mu}_{0}}{\hat{\rho}_{0} \sqrt{\hat{R} \hat{T}_{0}} L_{0}}, \end{equation}
where 
$\hat {\mu }_{0}$ is the viscosity of the gas at the reference state. Accordingly, the dimensionless and linearised conservation laws for the mass, momentum and energy read as
$$\begin{gather} \frac{\partial \rho}{\partial t} + \frac{\partial v_i}{\partial x_i}=0, \end{gather}$$
$$\begin{gather}\frac{\partial v_i} {\partial t} + \frac{\partial \sigma_{ij}}{\partial x_j} + \frac{\partial p}{\partial x_i} = F_i, \end{gather}$$
$$\begin{gather}\frac{3}{2} \frac{\partial T} {\partial t} + \frac{\partial q_i}{\partial x_i} + \frac{\partial v_i}{\partial x_i} = 0, \end{gather}$$
where 
$F_i$ is the dimensionless external force, 
$\partial (\cdot )/ \partial t$ denotes the dimensionless time derivative and 
$\partial (\cdot )/ \partial x_i$ denotes the dimensionless spatial derivative. In (2.3)–(2.5) and in what follows, indices denote vector and tensor components in Cartesian coordinates, and the Einstein summation convention is assumed over the repeated indices unless stated otherwise.
The higher moment system also includes the governing equations for the stress and heat flux, which in the linear-dimensionless form read as
$$\begin{gather} \frac{\partial \sigma_{ij}} {\partial t} + \frac{\partial m_{ijk}}{\partial x_k} + \frac{4}{5} \frac{\partial q_{\langle i}}{\partial x_{j\rangle}} ={-}\frac{1}{{Kn}} \left[\sigma_{ij} + 2 {Kn} \frac{\partial v_{\langle i}}{\partial x_{j\rangle}}\right], \end{gather}$$
$$\begin{gather}\frac{\partial q_i} {\partial t} +\frac{1}{2} \frac{\partial R_{ij}}{\partial x_j} +\frac{1}{6} \frac{\partial \varDelta}{\partial x_i} +\frac{\partial \sigma_{ij}}{\partial x_j} ={-}\frac{{Pr}}{{Kn}} \left[q_{i} + \frac{5}{2} \frac{{Kn}}{{Pr}} \frac{\partial T}{\partial x_i}\right], \end{gather}$$
where 
${Pr}$ is the Prandtl number and the angular brackets around indices denote the symmetric-tracefree part of a tensor (Struchtrup Reference Struchtrup2005). The right-hand sides of these balance equations contain NSF laws, which can also be obtained at the first order of a CE-like expansion performed on the stress and heat flux balance equations ((2.6) and (2.7)) in powers of the Knudsen number; see, e.g. Torrilhon & Struchtrup (Reference Torrilhon and Struchtrup2004) and Struchtrup (Reference Struchtrup2005). The first order of the CE-like expansion on these equations essentially amounts to setting their left-hand sides to zero, yielding the NSF laws.
The balance equations (2.3)–(2.7) form the so-called 13-moment system. However, this system is not closed since they contain additional higher moments 
$m_{ijk}$, 
$R_{ij}$ and 
$\varDelta$, which are defined as
$$\begin{gather} m_{ijk} =\int_{\mathbb{R}^{3}} C_{\langle i}C_{j}C_{k\rangle} f \, \mathrm{d}\boldsymbol{c}, \end{gather}$$
$$\begin{gather}R_{ij} =\int_{\mathbb{R}^{3}} C^2 C_{\langle i}C_{j\rangle }f\, \mathrm{d}\boldsymbol{c} - 7 \sigma_{ij}, \end{gather}$$
$$\begin{gather}\varDelta =\int_{\mathbb{R}^{3}}C^{4}f\, \mathrm{d}\boldsymbol{c} - 15, \end{gather}$$
where 
$C_{i}=c_{i}-v_{i}$ is the dimensionless peculiar velocity with 
$c_{i}$ being the dimensionless microscopic velocity. These higher moments are constructed in such a way that they vanish on being computed with the Grad 13-moment (G13) distribution function, i.e. 
$m_{ijk|{G13}}=R_{ij|{G13}}=\varDelta _{|{G13}} =0$.
The (linearised) 26-moment theory is comprised of the conservation laws (2.3)–(2.5), the stress balance equation (2.6), the heat flux balance equation (2.7) and the balance equations for 
$m_{ijk}$, 
$R_{ij}$ and 
$\varDelta$. The balance equations for 
$\hat {m}_{ijk}$, 
$\hat {R}_{ij}$ and 
$\hat {\varDelta }$ have already been obtained from the Boltzmann equation by Gu & Emerson (Reference Gu and Emerson2009); in the linear-dimensionless form, they read as
$$\begin{gather} \frac{\partial m_{ijk}} {\partial t} + \frac{\partial \varPhi_{ijkl}}{\partial x_l} +\frac{3}{7} \frac{\partial R_{\langle ij}}{\partial x_{k\rangle}} ={-}\frac{{Pr}_{m}}{{Kn}}\left[ m_{ijk}+3\frac{{Kn}}{{Pr}_{m}}\frac{\partial \sigma_{\langle ij}}{\partial x_{k\rangle}} \right], \end{gather}$$
$$\begin{gather}\frac{\partial R_{ij}} {\partial t} + \frac{\partial \varPsi_{ijk}}{\partial x_k} +2 \frac{\partial m_{ijk}}{\partial x_k} +\frac{2}{5} \frac{\partial \varOmega_{\langle i}}{\partial x_{j\rangle}} ={-}\frac{{Pr}_{R}}{{Kn}} \left[ R_{ij} + \frac{28}{5} \frac{{Kn}}{{Pr}_{R}}\frac{\partial q_{\langle i}}{\partial x_{j\rangle}} \right], \end{gather}$$
$$\begin{gather}\frac{\partial \varDelta} {\partial t} + \frac{\partial \varOmega_i}{\partial x_i} ={-}\frac{{Pr}_{\varDelta }}{{Kn}}\left[ \varDelta + 8 \frac{{Kn}}{{Pr}_{\varDelta}}\frac{\partial q_i}{\partial x_i}\right]. \end{gather}$$
The coefficients 
${Pr}$, 
${Pr}_{m}$, 
${Pr}_{R}$ and 
${Pr}_{\varDelta }$ in the above equations depend on the choice of intermolecular potential function appearing in the Boltzmann collision operator. For Maxwell molecules (MM), these parameters read as 
${Pr}=2/3$, 
${Pr}_{m}=3/2$, 
${Pr}_{R}=7/6$ and 
${Pr}_{\varDelta }=2/3$ (Gu & Emerson Reference Gu and Emerson2009).
For a third-order accuracy in 
${Kn}$, it is sufficient to consider only the terms on the right-hand sides of (2.11)–(2.13) and dropping all the terms on the left-hand sides of these equations; this gives the linear R13 constitutive relations (Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2003; Struchtrup Reference Struchtrup2005). Formally, the R13 constitutive relations stem from (2.11)–(2.13) on performing a CE-like expansion on these equations in powers of 
${Kn}$, i.e. by setting the left-hand sides of (2.11)–(2.13) to zero (Struchtrup Reference Struchtrup2004). Thus, the (dimensionless) LR13 equations consist of (2.3)–(2.7) closed with
\begin{equation} m_{ijk} ={-}3\frac{{Kn}}{{Pr}_{m}}\frac{\partial \sigma_{\langle ij}}{\partial x_{k\rangle}}, \quad R_{ij} ={-}\frac{28}{5} \frac{{Kn}}{{Pr}_{R}}\frac{\partial q_{\langle i}}{\partial x_{j\rangle}} \quad\text{and}\quad \varDelta ={-}8 \frac{{Kn}}{{Pr}_{\varDelta}}\frac{\partial q_i}{\partial x_i}. \end{equation}2.1. The R26 constitutive relations
The system of 26 moment equations (2.3)–(2.7) and (2.11)–(2.13) contains the additional moments 
$\varPhi _{ijkl}$, 
$\varPsi _{ijk}$ and 
$\varOmega _{i}$, i.e.
$$\begin{gather} \varPhi_{ijkl} = \int_{\mathbb{R}^{3}} C_{\langle i}C_{j}C_{k}C_{l\rangle }f\, \mathrm{d}\boldsymbol{c} , \end{gather}$$
$$\begin{gather}\varPsi_{ijk} = \int_{\mathbb{R}^{3}}C^{2}C_{\langle i}C_{j}C_{k\rangle }f\, \mathrm{d}\boldsymbol{c} -9 m_{ijk}, \end{gather}$$
$$\begin{gather}\varOmega_{i} = \int_{\mathbb{R}^{3}} C^{4}C_{i}\,f\, \mathrm{d}\boldsymbol{c} - 28 q_{i} , \end{gather}$$
which need to be specified in order to close the system of the 26 moment equations. The R26 constitutive relations, which provides a closure for the 26 moment equations, follow from a CE-like expansion on the evolution equations for the moments 
$\varPhi _{ijkl}$, 
$\varPsi _{ijk}$ and 
$\varOmega _{i}$ in a 45-moment theory (see Gu & Emerson Reference Gu and Emerson2009). The R26 constitutive relations read as (Gu & Emerson Reference Gu and Emerson2009)
\begin{equation} \left. \begin{aligned} \varPhi_{ijkl} & ={-}4\frac{{Kn}}{{Pr}_{\varPhi}}\frac{\partial m_{\langle ijk}}{\partial x_{l\rangle}}, \\ \varPsi_{ijk} & ={-}\frac{27}{7}\frac{{Kn}}{{Pr}_{\varPsi}}\frac{\partial R_{\langle ij}}{\partial x_{k\rangle}}, \\ \varOmega_{i} & ={-}\frac{7}{3}\frac{{Kn}}{{Pr}_{\varOmega}} \left(\frac{\partial \varDelta}{\partial x_i} + \frac{12}{7}\frac{\partial R_{ij}}{\partial x_j}\right), \end{aligned} \right\} \end{equation}
where the values of the parameters are 
${Pr}_{\varPhi }=2.097$, 
${Pr}_{\varPsi }=1.698$ and 
${Pr}_{\varOmega }=1$ for MM (Gu & Emerson Reference Gu and Emerson2009). The system of 26 moment equations (2.3)–(2.7) and (2.11)–(2.13) along with the constitutive relations (2.18) form the system of the LR26 equations.
3. 
$H$-theorem for the linearised R26 equations
The 
$H$-theorem plays a salient role in constructing constitutive relations in the bulk (de Groot & Mazur Reference de Groot and Mazur1962; Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2007) and in developing physically admissible boundary conditions (Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2007; Rana & Struchtrup Reference Rana and Struchtrup2016; Rana, Gupta & Struchtrup Reference Rana, Gupta and Struchtrup2018a; Sarna & Torrilhon Reference Sarna and Torrilhon2018a). Recently, Struchtrup & Nadler (Reference Struchtrup and Nadler2020) deduced that the Burnett equations, which do not comply with the 
$H$-theorem, lead to thermodynamically inadmissible solutions (producing work without a driving force). In this section, we shall – for the first time – show that the LR26 equations are accompanied with an entropy law given as
\begin{equation} \frac{\partial \eta}{\partial t} + \frac{\partial \varGamma_i}{\partial x_i} = \varSigma \geqslant 0, \end{equation}
where 
$\eta$ denotes the dimensionless entropy density, 
$\varGamma _{i}$ the dimensionless entropy flux and 
$\varSigma$ the dimensionless non-negative entropy generation rate.
The Boltzmann entropy density 
$\hat {\eta }$ and entropy flux 
$\hat {\varGamma }_i$ (in dimensional form) are given by (Müller & Ruggeri Reference Müller and Ruggeri1998)
$$\begin{gather} \hat{\eta}={-}\hat{k}_{B}\int \hat{f} \left( \ln \frac{\hat{f}}{\hat{y}}-1\right) \,\mathrm{d}\hat{\boldsymbol{c}}, \end{gather}$$
$$\begin{gather}\hat{\varGamma}_i ={-}\hat{k}_{B}\int \hat{f} \left( \ln \frac{\hat{f}}{\hat{y}}-1\right) \hat{c}_i \, \mathrm{d}\hat{\boldsymbol{c}}, \end{gather}$$
where 
$\hat {k}_{B}$ is the Boltzmann constant, 
$\hat {y}$ is a constant having dimensions of the distribution function 
$\hat {f}$ and 
$\hat {\boldsymbol {c}}$ is the microscopic velocity (with dimensions). Starting from the Boltzmann entropy density (3.2), it can be shown that the entropy density of a system of the linearised moment equations – obtained from the Boltzmann equation – is a quadratic form without cross-product terms (see, e.g. Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2007; Rana & Struchtrup Reference Rana and Struchtrup2016; Beckmann et al. Reference Beckmann, Rana, Torrilhon and Struchtrup2018; Rana et al. Reference Rana, Gupta and Struchtrup2018a; Sarna & Torrilhon Reference Sarna and Torrilhon2018b). It is worthwhile noting that the hyperbolic part of the LR26 equations is the same as the linearised Grad 26-moment (G26) equations (in the same way the hyperbolic part of the linearised NSF equations is the same as the linearised Euler equations); therefore, it is sufficient to use the G26 distribution function for computing the entropy density for the LR26 equations. Using the G26 distribution function for 
$\hat {f}$ in (3.2) and (3.3), we show that entropy density for the LR26 equations is a quadratic form without cross-product terms and that the entropy flux obtained with the G26 distribution function contains all the hyperbolic terms present in the expression of the entropy flux for the LR26 equations, although the calculations are relegated to Appendix A for better readability. The dimensionless entropy density for the LR26 equations is given by (see appendices A and B)
\begin{equation} \eta =\underline{a_{0}-\tfrac{1}{2}\rho ^{2}-\tfrac{1}{2} v^{2} -\tfrac{3}{4} T^{2}-\tfrac{1}{4}\sigma ^{2}-\tfrac{1}{5}q^{2}} - \tfrac{1}{12}m^{2} - \tfrac{1}{56}R^{2} - \tfrac{1}{240}\varDelta^{2}, \end{equation}
where 
$a_0$ is just a constant and the notation 
$A^{2}:=A_{i_{1}i_{2},\ldots , i_{k}}A_{i_{1}i_{2},\ldots , i_{k}}$ has been used to denote the contraction of a 
$k$-rank tensor 
$A_{i_{1}i_{2}\ldots i_{k}}$ with itself for brevity. The entropy for the LR26 equations (3.4) comprises of the contribution from the 13-moment theory (Rana & Struchtrup Reference Rana and Struchtrup2016) – given by the underlined terms in (3.4) – and additional contributions from the higher-order moments.
The dimensionless entropy density (3.4) can also be obtained in a more comprehensible way by assuming it to be a quadratic form without cross-product terms. The approach not only leads to the dimensionless entropy density (3.4) but also yields the complete entropy flux and entropy production rate for the LR26 equations alongside. The approach assumes the following form of dimensionless entropy density:
\begin{equation} \eta = a_0 + \frac{a_1}{2}\rho ^{2} + \frac{a_2}{2} v^{2} + \frac{a_3}{2} T^{2} + \frac{a_4}{2} \sigma^{2} + \frac{a_5}{2} q^{2} + \frac{a_6}{2} m^{2} + \frac{a_7}{2} R^{2} + \frac{a_8}{2}\varDelta^{2}. \end{equation}
Here 
$a_0,a_1,a_2,\ldots ,a_8$ are unknown coefficients, which are computed by taking the time derivative of the entropy density (3.5), then eliminating the time derivatives of the moments by means of the G26 equations, and finally comparing the resulting equation with (3.1). It turns out that one of the variables from 
$a_1, a_2, \ldots , a_8$ can be chosen arbitrarily. On taking 
$a_1 = -1$, the other coefficients 
$a_2, a_3, \ldots , a_8$ turn out to be such that (3.5) is the same as (3.4); see Appendix B for details.
For the computation of the entropy flux and entropy production rate for the LR26 equations, we therefore can start with the entropy density (3.4). Taking the time derivative of 
$\eta$ in (3.4), and substituting the time derivatives of the field variables 
$\rho$, 
$v_{i}$, 
$T$ from (2.3)–(2.5), the time derivative of 
$\sigma _{ij}$ from (2.6), the time derivative of 
$q_{i}$ from (2.7), and the time derivatives of 
$m_{ijk}$, 
$R_{ij}$, 
$\varDelta$ from (2.11)–(2.13), we obtain
$$\begin{align} \frac{\partial
\eta}{\partial t} &= \frac{\partial}{\partial x_i} \left(p
v_{i} + T q_{i} + v_{j}\sigma_{ij} +\frac{2}{5}\sigma_{ij}
q_j +\frac{1}{2} m_{ijk} \sigma_{jk}
+\frac{1}{5}R_{ij}q_{j} +\frac{1}{15}q_{i} \varDelta
\right.\nonumber\\ &\quad\left. +\frac{1}{14} m_{ijk}R_{jk}
+\frac{1}{6}\varPhi_{ijkl}m_{jkl}
+\frac{1}{28}\varPsi_{ijk}R_{jk}
+\frac{1}{70}R_{ij}\varOmega_{j} +\frac{1}{120}\varDelta
\varOmega_{i} \right) \nonumber\\
&\quad -\frac{1}{6}\varPhi_{ijkl} \frac{\partial m_{jkl}}{\partial
x_i} -\frac{1}{28}\varPsi_{ijk} \frac{\partial
R_{jk}}{\partial x_i} -\frac{1}{70}\varOmega_{j}
\frac{\partial R_{ij}}{\partial x_i}
-\frac{1}{120}\varOmega_{i} \frac{\partial
\varDelta}{\partial x_i} \nonumber\\
&\quad+\frac{1}{2}\frac{1}{{Kn}} \sigma^{2} +\frac{2}{5}
\frac{{Pr}}{{Kn}} q^{2} +\frac{1}{6}\frac{{Pr}_{m}}{{Kn}}
m^{2} +\frac{1}{28}\frac{{Pr}_{R}}{{Kn}} R^{2}
+\frac{1}{120}\frac{{Pr}_{\varDelta}}{{Kn}} \varDelta^{2}.
\end{align}$$
Comparison of (3.6) with (3.1) gives the expression for the entropy flux as
$$\begin{align} \varGamma_{i} = &- \left( p v_{i} + T q_{i} + v_{j}\sigma_{ij} +\tfrac{2}{5}\sigma_{ij} q_j +\tfrac{1}{2} m_{ijk} \sigma_{jk} +\tfrac{1}{5}R_{ij}q_{j} +\tfrac{1}{15}q_{i} \varDelta \right. \nonumber\\ &\left. +\tfrac{1}{14} m_{ijk}R_{jk} +\tfrac{1}{6}\varPhi_{ijkl}m_{jkl} +\tfrac{1}{28}\varPsi_{ijk}R_{jk} +\tfrac{1}{70}R_{ij}\varOmega_{j} +\tfrac{1}{120}\varDelta \varOmega_{i}\right), \end{align}$$and the entropy production rate as
\begin{align} \varSigma &= -
\frac{1}{6}\varPhi_{ijkl} \frac{\partial m_{jkl}}{\partial
x_i} -\frac{1}{28}\varPsi_{ijk} \frac{\partial
R_{jk}}{\partial x_i} -\frac{1}{70}\varOmega_{j}
\frac{\partial R_{ij}}{\partial x_i}
-\frac{1}{120}\varOmega_{i} \frac{\partial
\varDelta}{\partial x_i}\nonumber\\
&\quad +\frac{1}{2}\frac{1}{{Kn}} \sigma^{2} +\frac{2}{5}
\frac{{Pr}}{{Kn}} q^{2} +\frac{1}{6}\frac{{Pr}_{m}}{{Kn}}
m^{2} +\frac{1}{28}\frac{{Pr}_{R}}{{Kn}} R^{2}
+\frac{1}{120}\frac{{Pr}_{\varDelta}}{{Kn}} \varDelta^{2}.
\end{align}
Using the linearised R26 constitutive relations (2.18) in (3.8), the entropy production rate 
$\varSigma$ can be written as
\begin{align} \varSigma &=
\frac{1}{24} \frac{{Kn}}{{Pr}_{\varPhi}} \varPhi^2
+\frac{1}{108} \frac{{Kn}}{{Pr}_{\varPsi}} \varPsi^2 +
\frac{1}{280} \frac{{Kn}}{{Pr}_{\varOmega}} \varOmega^2 \nonumber\\ & \quad +\frac{1}{2}\frac{1}{{Kn}} \sigma^{2}
+\frac{2}{5} \frac{{Pr}}{{Kn}} q^{2}
+\frac{1}{6}\frac{{Pr}_{m}}{{Kn}} m^{2}
+\frac{1}{28}\frac{{Pr}_{R}}{{Kn}} R^{2}
+\frac{1}{120}\frac{{Pr}_{\varDelta}}{{Kn}} \varDelta^{2},
\end{align}
which is non-negative for 
${Pr}_{\varPsi }$, 
${Pr}_{\varPhi }$, 
${Pr}_{\varOmega }$, 
${Pr}_{m}$, 
${Pr}_{R}$, 
${Pr}_{\varDelta }$, 
$\Pr$, 
${Kn}\geq 0$. This completes the proof of the 
$H$-theorem for the linearised R26 equations.
4. Phenomenological boundary conditions for the LR26 equations
In this section, we shall demonstrate how the 
$H$-theorem can be used to generate a proper and complete set of boundary conditions for the LR26 equations. We shall derive the PBC for the LR26 equations from the entropy production rate at the interface exploiting the force-flux relationships.
Let there be two phases, namely the liquid (or solid) and the gaseous (or vapour), connected by a massless interface 
$I$ of negligible thickness and let there be no surface tension and surface energy at the interface. The entropy production rate at such an interface 
$\varSigma ^{I}$ is given by the difference between the entropy fluxes into and out of an interface (de Groot & Mazur Reference de Groot and Mazur1962; Rana & Struchtrup Reference Rana and Struchtrup2016; Beckmann et al. Reference Beckmann, Rana, Torrilhon and Struchtrup2018), i.e.
\begin{equation} \varSigma^{I}=\left[\varGamma_{i}^{{(gas)}} -\varGamma_{i}^{{(liquid)}}\right] n_{i}\geq 0, \end{equation}
where 
$n_{i}$ is the unit normal to the interface pointing from the liquid into the gas. Assuming that the gaseous (or vapour) phase is comprised of an ideal gas, the entropy flux of the gas (or vapour) medium is 
$\varGamma _{i}^{{(gas)}} = \varGamma _{i}$, given by (3.7), and assuming that the liquid phase is comprised of an incompressible liquid and can be described by the NSF equations, the entropy flux of the liquid phase is (Beckmann et al. Reference Beckmann, Rana, Torrilhon and Struchtrup2018)
\begin{equation} \varGamma_{i}^{{(liquid)}} ={-} p^{\ell} v_{i}^{\ell} - T^{\ell} q_{i}^{\ell} - v_{j}^{\ell}\sigma_{ij}^{\ell}, \end{equation}
where the superscript ‘
$\ell$’ on the variables are used to denote that they are the properties of the liquid. Substituting 
$\varGamma _{i}^{{(gas)}} = \varGamma _{i}$ from (3.7) and 
$\varGamma _{i}^{{(liquid)}}$ from (4.2) into (4.1), the entropy flux at the interface, after some algebra, turns out to be (see Appendix C for details)
$$\begin{align} \varSigma^{I}&={-}
\left(\mathcal{P} \mathcal{V}_i + \mathcal{T} q_{i} +
\mathscr{V}_j \sigma_{ij} +\tfrac{2}{5}\sigma_{ij} q_j
+\tfrac{1}{2} m_{ijk} \sigma_{jk} +\tfrac{1}{5}R_{ij}q_{j}
+\tfrac{1}{15}q_{i} \varDelta \right.\nonumber\\
&\left.\quad\, +\tfrac{1}{14} m_{ijk}R_{jk}
+\tfrac{1}{6}\varPhi_{ijkl}m_{jkl}
+\tfrac{1}{28}\varPsi_{ijk}R_{jk}
+\tfrac{1}{70}R_{ij}\varOmega_{j} +\tfrac{1}{120}
\varOmega_{i}\varDelta\right)n_{i},
\end{align}$$
where
\begin{equation} \mathcal{P} = p - p_{{sat}}, \quad \mathcal{T} = T - T^{\ell}, \quad \mathcal{V}_{i}=v_{i}-v_{i}^{I} \quad\text{and}\quad \mathscr{V}_j = v_j - v_j^{\ell}, \end{equation}
with 
$v_{i}^{I}$ being the velocity of the interface and 
$p_{{sat}} \equiv p_{{sat}} (T^{\ell })$ being the saturation pressure corresponding to the temperature 
$T^{\ell }$.
For further simplification, it is imperative to decompose the vectors and tensors into their components in the normal and tangential directions. Such a decomposition for symmetric-tracefree tensors of rank up to three is given in Appendix D. This decomposition allows one to write the entropy production at the wall (4.3) as a sum of the product of the fluxes with an odd degree in 
$n$ (i.e. 
$\mathcal {V}_n$, 
$q_{n}$, 
$\varOmega _{n}$, 
$m_{nnn}$, 
$\varPsi _{nnn}$, 
$\bar {\sigma }_{ni}$, 
$\bar {R}_{ni}$, 
$\bar {\varPhi }_{nnni}$, 
$\tilde {m}_{nij}$, 
$\tilde {\varPsi }_{nij}$, 
$\check {\varPhi }_{nijk}$) and the moments with an even degree in 
$n$ as
$$\begin{align} \varSigma^{I}
=&- \mathcal{V}_n (\mathcal{P}+\sigma_{nn}) -q_{n}\left(
\mathcal{T}+\tfrac{2}{5}\sigma_{nn}+\tfrac{1}{5}R_{nn}
+\tfrac{1}{15}\varDelta \right) \nonumber\\
&-m_{nnn}\left( \tfrac{3}{4}\sigma_{nn}+\tfrac{3}{28}
R_{nn}+\tfrac{5}{12}\varPhi_{nnnn}\right) -\varPsi_{nnn}
\left(\tfrac{3}{56} R_{nn}\right) \nonumber\\ &\quad
-\varOmega_{n}\left( \tfrac{1}{120} \varDelta +
\tfrac{1}{70} R_{nn} \right) -\bar{\sigma}_{ni} \left(
\bar{\mathscr{V}}_{i}+\tfrac{2}{5}\bar{q}_{i}
+\bar{m}_{nni}\right) \nonumber\\ &\quad -\bar{R}_{ni}\left(
\tfrac{1}{5}\bar{q}_{i} + \tfrac{1}{7} \bar{m}_{nni}
+\tfrac{1}{14} \bar{\varPsi}_{nni} +\tfrac{1}{70}
\bar{\varOmega}_{i}\right) - \bar{\varPhi}_{nnni}
\left(\tfrac{5}{8}\bar{m}_{nni}\right) \nonumber\\ &\quad
-\tilde{m}_{nij} \left(\tfrac{1}{2}
\tilde{\sigma}_{ij}+\tfrac{1}{14}\tilde{R}_{ij} +
\tfrac{1}{2} \tilde{\varPhi}_{nnij} \right)
-\tilde{\varPsi}_{nij} \left(\tfrac{1}{28}
\tilde{R}_{ij}\right) -\check{\varPhi}_{nijk}
\left(\tfrac{1}{6} \check{m}_{ijk}\right).
\end{align}$$
While writing (4.5), the relation 
$\sigma _{nn} \mathscr {V}_n \approx \sigma _{nn} \mathcal {V}_n$ in light of (C7) and the Curie principle – which states that only forces and fluxes of the same tensor type (scalars, vectors, 2-tensors, etc.) can be combined (de Groot & Mazur Reference de Groot and Mazur1962) – have been exploited. In the above equation, 
$\{ \mathcal {V}_n$, 
$q_{n}$, 
$m_{nnn}$, 
$\varPsi _{nnn}$, 
$\varOmega _{n}\}$ are the odd (in 
$n$) scalar fluxes, 
$\{\bar {\sigma }_{ni}$, 
$\bar {R}_{ni}$, 
$\bar {\varPhi }_{nnni}\}$ are the odd vector fluxes, 
$\{\tilde {m}_{nij}$, 
$\tilde {\varPsi }_{nij}\}$ are the odd tensor fluxes of rank two and 
$\check {\varPhi }_{nijk}$ is the odd tensor flux of rank three. The positive entropy production at the interface can be attained by taking the odd fluxes proportional to the even (in 
$n$) driving forces (written in brackets in (4.5)) – with proportionality constants being positive. This step to obtain the positive entropy production at the interface yields the required boundary conditions. Mathematically, the positive entropy production, and, hence, the thermodynamically consistent PBC at the interface are achieved by writing a linear force-flux relationship of the form
\begin{equation} \mathcal{J}_{r}={-} \mathcal{L}_{rs}\mathcal{F}_{s} \end{equation}
for scalar, vector and tensor fluxes of each rank, successively. Here, 
$\mathcal {J}_{r}$ is the 
$r^{\mathrm{th}}$ component of a vector 
$\boldsymbol {\mathcal {J}}$ containing the fluxes, 
$\mathcal {F}_{s}$ is the 
$s^{\mathrm{th}}$ component of a vector 
$\boldsymbol {\mathcal {F}}$ containing the forces, 
$\mathcal {L}_{rs}$ is the 
$(r,s)^\mathrm{th}$ element of an unknown matrix 
$\boldsymbol {\mathcal {L}}$ and Einstein summation is assumed over 
$s$ in (4.6). The unknown matrix 
$\boldsymbol {\mathcal {L}}$ consists of phenomenological coefficients as its elements. The matrix 
$\boldsymbol {\mathcal {L}}$ must be symmetric in order to comply with the Onsager reciprocity relations (Onsager Reference Onsager1931a,Reference Onsagerb) and must be positive semidefinite to ensure a non-negative entropy production.
4.1. Boundary conditions for the scalar fluxes
The scalar fluxes in (4.5) and their corresponding forces are collected in 
$\boldsymbol {\mathcal {J}}$ and 
$\boldsymbol {\mathcal {F}}$ as
\begin{equation}
\renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{J}}
= \begin{bmatrix} \mathcal{V}_n \\ q_{n} \\ m_{nnn} \\
\varPsi_{nnn} \\ \varOmega_{n} \end{bmatrix}
\quad\text{and}\quad \boldsymbol{\mathcal{F}} =
\begin{bmatrix} \mathcal{P}+\sigma_{nn} \\
\mathcal{T}+\frac{2}{5}\sigma_{nn}+\frac{1}{5}R_{nn}
+\frac{1}{15}\varDelta \\
\frac{3}{4}\sigma_{nn}+\frac{3}{28}R_{nn}
+\frac{5}{12}\varPhi_{nnnn} \\ \frac{3}{56} R_{nn} \\
\frac{1}{120} \varDelta + \frac{1}{70} R_{nn}
\end{bmatrix}.
\end{equation}
Substituting 
$\boldsymbol {\mathcal {J}}$ and 
$\boldsymbol {\mathcal {F}}$ from (4.7a,b) in (4.6), we obtain the PBC for the scalar fluxes. Nevertheless, the matrix 
$\boldsymbol {\mathcal {L}}$ is still unknown. The elements of the matrix 
$\boldsymbol {\mathcal {L}}$ can be obtained by comparing the coefficients of 
$\mathcal {P}$, 
$\mathcal {T}$, 
$\sigma _{nn}$, 
$R_{nn}$ and 
$\varDelta$ in the corresponding terms of the PBC and MBC for the scalar fluxes. The MBC for the scalar fluxes (obtained through the Maxwell accommodation model) read as (Rana et al. Reference Rana, Lockerby and Sprittles2018b)
\begin{equation} \left. \begin{aligned} \mathcal{V}_n & ={-} \varsigma \left(\mathcal{P} - \frac{1}{2} \mathcal{T} + \frac{1}{2} \sigma_{nn} - \frac{1}{28} R_{nn} - \frac{1}{120} \varDelta +b_1 \varPhi_{nnnn}\right), \\ q_n & ={-}\frac{\mathcal{V}_n}{2} - \varkappa \left(2 \mathcal{T} + \frac{1}{2} \sigma_{nn} + \frac{5}{28} R_{nn} + \frac{1}{15} \varDelta +b_2 \varPhi_{nnnn} \right), \\ m_{nnn} & ={-}\frac{2}{5} \mathcal{V}_n + \varkappa \left(\frac{2}{5} \mathcal{T} - \frac{7}{5} \sigma_{nn} - \frac{1}{14} R_{nn} + \frac{1}{75} \varDelta +b_3 \varPhi_{nnnn}\right), \\ \varPsi_{nnn} & = \frac{6}{5} \mathcal{V}_n + \varkappa \left(\frac{6}{5} \mathcal{T} + \frac{9}{5} \sigma_{nn} - \frac{93}{70} R_{nn} + \frac{1}{5} \varDelta +b_4 \varPhi_{nnnn}\right), \\ \varOmega_n & = 3 \mathcal{V}_n + \varkappa \left(8 \mathcal{T} + 2 \sigma_{nn} - R_{nn} - \frac{4}{3} \varDelta +b_5 \varPhi_{nnnn}\right), \end{aligned} \right\} \end{equation}where
\begin{gather} b_1 ={-}\tfrac{1}{24}, \quad b_2 ={-} \tfrac{1}{12}, \quad b_3 ={-} \tfrac{13}{30}, \quad b_4 ={-} \tfrac{1}{2}, \quad b_5 = 0, \end{gather}
\begin{gather} \varsigma = \frac{\vartheta}{2-\vartheta} \sqrt{\frac{2}{\rm \pi}} \quad\text{and} \quad \varkappa = \frac{\vartheta + \chi (1-\vartheta)}{2-\vartheta - \chi (1-\vartheta)} \sqrt{\frac{2}{\rm \pi}}, \end{gather}
with 
$\vartheta$ being the evaporation/condensation coefficient and 
$\chi$ being the accommodation coefficient. Clearly, for 
$\vartheta =0$, 
$\mathcal {V}_n=0$ from (4.8); hence, we obtain the boundary conditions for a non-evaporating boundary (Gu & Emerson Reference Gu and Emerson2009). It is worthwhile noting that the MBC (4.8) with (4.9a–e) cannot be expressed in the force-flux formalism (4.6), and we shall show later that this leads to violation of the Onsager reciprocity relations and/or the second law at the interface.
Comparison of the coefficients of 
$\mathcal {P}$, 
$\mathcal {T}$, 
$\sigma _{nn}$, 
$R_{nn}$ and 
$\varDelta$ in the PBC for the scalar fluxes (obtained on substituting 
$\boldsymbol {\mathcal {J}}$ and 
$\boldsymbol {\mathcal {F}}$ from (4.7a,b) in (4.6)) and MBC (4.8) yields
\begin{equation} \renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{L}} = \varsigma \begin{bmatrix} 1 & -\frac{1}{2} & -\frac{2}{5} & \frac{6}{5} & 3 \\ -\frac{1}{2} & \frac{1}{4} & \frac{1}{5} & -\frac{3}{5} & -\frac{3}{2} \\ -\frac{2}{5} & \frac{1}{5} & \frac{4}{25} & -\frac{12}{25} & -\frac{6}{5} \\ \frac{6}{5} & -\frac{3}{5} & -\frac{12}{25} & \frac{36}{25} & \frac{18}{5} \\ 3 & -\frac{3}{2} & -\frac{6}{5} & \frac{18}{5} & 9 \end{bmatrix} +\varkappa \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & -\frac{2}{5} & -\frac{6}{5} & -8 \\ 0 & -\frac{2}{5} & \frac{52}{25} & -\frac{44}{25} & \frac{8}{5} \\ 0 & -\frac{6}{5} & -\frac{44}{25} & \frac{916}{25} & -\frac{72}{5} \\ 0 & -8 & \frac{8}{5} & -\frac{72}{5} & 224 \end{bmatrix}. \end{equation}
Both the matrices on the right-hand side of (4.11) are symmetric and positive semidefinite. Therefore, the matrix 
$\boldsymbol {\mathcal {L}}$ in (4.11) is also symmetric and positive semidefinite as 
$\varsigma$ and 
$\varkappa$ are non-negative. Thus, the PBC for the scalar fluxes are obtained from (4.6) with 
$\boldsymbol {\mathcal {J}}$, 
$\boldsymbol {\mathcal {F}}$ and 
$\boldsymbol {\mathcal {L}}$ being given by (4.7a,b) and (4.11). On simplification, the PBC for the scalar fluxes are given by (4.8) with
\begin{equation} b_1 ={-}\tfrac{1}{6}, \quad b_2 ={-} \tfrac{1}{6}, \quad b_3 ={-} \tfrac{13}{15}, \quad b_4 = \tfrac{11}{15}, \quad b_5 ={-} \tfrac{2}{3}. \end{equation}4.2. Boundary conditions for the vector fluxes
The vector fluxes in (4.5) and their corresponding forces are collected in 
$\boldsymbol {\mathcal {J}}$ and 
$\boldsymbol {\mathcal {F}}$ as
\begin{equation} \renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{J}} = \begin{bmatrix} \bar{\sigma}_{ni} \\ \bar{R}_{ni} \\ \bar{\varPhi}_{nnni} \end{bmatrix} \quad\text{and}\quad \boldsymbol{\mathcal{F}} = \begin{bmatrix} \bar{\mathscr{V}}_{i}+\frac{2}{5}\bar{q}_{i} +\bar{m}_{nni} \\ \frac{1}{5}\bar{q}_{i} + \frac{1}{7} \bar{m}_{nni} +\frac{1}{14} \bar{\varPsi}_{nni} +\frac{1}{70} \bar{\varOmega}_{i} \\ \frac{5}{8}\bar{m}_{nni} \end{bmatrix}. \end{equation}
Substituting 
$\boldsymbol {\mathcal {J}}$ and 
$\boldsymbol {\mathcal {F}}$ from (4.13a,b) in (4.6), we obtain the PBC for the vector fluxes. Nevertheless, the matrix 
$\boldsymbol {\mathcal {L}}$ is still unknown. The elements of the matrix 
$\boldsymbol {\mathcal {L}}$ can be obtained by comparing the coefficients of 
$\bar {\mathscr {V}}_{i}$, 
$\bar {q}_{i}$ and 
$\bar {m}_{nni}$ in the corresponding terms of the PBC and MBC. The MBC for the vector fluxes (obtained through the Maxwell accommodation model) read as (Gu & Emerson Reference Gu and Emerson2009; Struchtrup et al. Reference Struchtrup, Beckmann, Rana and Frezzotti2017)
\begin{equation} \left. \begin{aligned} \bar{\sigma}_{ni} & ={-} \varkappa \left( \bar{\mathscr{V}}_{i} +\tfrac{1}{5}\bar{q}_{i} +\tfrac{1}{2}\bar{m}_{nni} + b_6 \bar{\varPsi}_{nni} + b_7 \bar{\varOmega}_{i}\right), \\ \bar{R}_{ni} & ={-} \varkappa \left( -\bar{\mathscr{V}}_{i} + \tfrac{11}{5} \bar{q}_{i} + \tfrac{1}{2} \bar{m}_{nni} + b_8 \bar{\varPsi}_{nni} + b_9 \bar{\varOmega}_{i}\right), \\ \bar{\varPhi}_{nnni} & ={-} \varkappa \left( -\tfrac{4}{7}\bar{\mathscr{V}}_{i} -\tfrac{12}{35} \bar{q}_{i}+\tfrac{9}{7} \bar{m}_{nni} + b_{10} \bar{\varPsi}_{nni} + b_{11} \bar{\varOmega}_{i}\right), \end{aligned} \right\} \end{equation}where
\begin{align} b_6 ={-}\tfrac{1}{36}, \quad b_7 ={-}\tfrac{1}{280}, \quad b_8 = \tfrac{5}{12}, \quad b_9 = \tfrac{5}{56}, \quad b_{10} = \tfrac{1}{14}, \quad b_{11} ={-} \tfrac{3}{490}. \end{align}
Comparison of the coefficients of 
$\bar {\mathscr {V}}_{i}, \bar {q}_{i}, \bar {m}_{nni}$ in the PBC for the vector fluxes (obtained on substituting 
$\boldsymbol {\mathcal {J}}$ and 
$\boldsymbol {\mathcal {F}}$ from (4.13a,b) in (4.6)) and MBC (4.14) yields
\begin{equation} \renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{L}} = \varkappa \begin{bmatrix} 1 & -1 & -\frac{4}{7} \\ -1 & 13 & -\frac{4}{7} \\ -\frac{4}{7} & -\frac{4}{7} & \frac{152}{49} \end{bmatrix}. \end{equation}
The matrix 
$\boldsymbol {\mathcal {L}}$ in (4.16) is symmetric and positive semidefinite as 
$\varkappa$ is non-negative. Thus, the PBC for the vector fluxes are obtained from (4.6) with 
$\boldsymbol {\mathcal {J}}$, 
$\boldsymbol {\mathcal {F}}$ and 
$\boldsymbol {\mathcal {L}}$ being given by (4.13a,b) and (4.16). On simplification, the PBC for the vector fluxes are given by (4.14) with
\begin{equation} b_6 ={-}\tfrac{1}{14}, \quad b_7 ={-}\tfrac{1}{70}, \quad b_8 = \tfrac{13}{14}, \quad b_9 = \tfrac{13}{70}, \quad b_{10} ={-} \tfrac{2}{49}, \quad b_{11} ={-} \tfrac{2}{245}. \end{equation}4.3. Boundary conditions for the tensor fluxes
4.3.1. Boundary conditions for the rank-2 tensor fluxes
The rank-2 tensor fluxes in (4.5) and their corresponding forces are collected in 
$\boldsymbol {\mathcal {J}}$ and 
$\boldsymbol {\mathcal {F}}$ as
\begin{equation} \renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{J}} = \begin{bmatrix} \tilde{m}_{n ij} \\ \tilde{\varPsi}_{n ij} \end{bmatrix} \quad\text{and}\quad \boldsymbol{\mathcal{F}} = \begin{bmatrix} \frac{1}{2} \tilde{\sigma}_{ij} + \frac{1}{14}\tilde{R}_{ij} + \frac{1}{2} \tilde{\varPhi}_{nn ij} \\ \frac{1}{28} \tilde{R}_{ij} \end{bmatrix}. \end{equation}
Substituting 
$\boldsymbol {\mathcal {J}}$ and 
$\boldsymbol {\mathcal {F}}$ from (4.18a,b) in (4.6), we obtain the PBC for the rank-2 tensor fluxes. Nevertheless, the matrix 
$\boldsymbol {\mathcal {L}}$ is still unknown. The elements of the matrix 
$\boldsymbol {\mathcal {L}}$ can be obtained by comparing the coefficients of 
$\tilde {\sigma }_{ij}$ and 
$\tilde {R}_{ij}$ in the corresponding terms of the PBC and MBC for the rank-2 tensor fluxes.
The MBC for the rank-2 tensor fluxes (obtained through the Maxwell accommodation model) read as
\begin{equation} \left. \begin{aligned} \tilde{m}_{nij} & ={-} \varkappa \left( \tilde{\sigma}_{ij} + \tfrac{1}{14}\tilde{R}_{ij} + b_{12} \tilde{\varPhi}_{nn ij} \right), \\ \tilde{\varPsi}_{nij} & = \varkappa \left( \tilde{\sigma}_{ij} - \tfrac{15}{14}\tilde{R}_{ij} + b_{13} \tilde{\varPhi}_{nn ij} \right), \end{aligned} \right\} \end{equation}where
\begin{equation} b_{12} = \tfrac{1}{2} \quad\text{and}\quad b_{13} ={-} \tfrac{1}{2}. \end{equation}
Comparison of the coefficients of 
$\tilde {\sigma }_{ij}$ and 
$\tilde {R}_{ij}$ in the PBC for the rank-2 tensor fluxes (obtained on substituting 
$\boldsymbol {\mathcal {J}}$ and 
$\boldsymbol {\mathcal {F}}$ from (4.18a,b) in (4.6)) and MBC (4.19) yields
\begin{equation} \renewcommand{\arraystretch}{1.5} \boldsymbol{\mathcal{L}} = \varkappa \begin{bmatrix} 2 & -2 \\ -2 & 34 \end{bmatrix}. \end{equation}
The matrix 
$\boldsymbol {\mathcal {L}}$ in (4.16) is symmetric and positive semidefinite as 
$\varkappa$ is non-negative. Thus, the PBC for the rank-2 tensor fluxes are obtained from (4.6) with 
$\boldsymbol {\mathcal {J}}$, 
$\boldsymbol {\mathcal {F}}$ and 
$\boldsymbol {\mathcal {L}}$ being given by (4.18a,b) and (4.21). On simplification, the PBC for the rank-2 tensor fluxes are given by (4.19) with
\begin{equation} b_{12} = 1 \quad\text{and}\quad b_{13} = 1. \end{equation}4.3.2. Boundary conditions for the rank-3 tensor fluxes
The rank-3 tensor flux in (4.5) and its corresponding forces are
\begin{equation} J_r = \check{\varPhi}_{n ijk} \quad\text{and}\quad \mathcal{F}_s = \tfrac{1}{6} \check{m}_{ijk}. \end{equation}
Substituting 
$\mathcal {J}_r$ and 
$\mathcal {F}_s$ from (4.23a,b) in (4.6), we obtain the PBC for the rank-3 tensor flux. Nevertheless, the coefficient 
$\mathcal {L}_{rs}$ in (4.6) is still unknown but can be obtained by comparing the coefficients of 
$\check {m}_{ijk}$ in the PBC and MBC for the rank-3 tensor flux.
The MBC for rank-3 tensor flux (obtained through the Maxwell accommodation model) read as
\begin{equation} \check{\varPhi}_{nijk} ={-}\varkappa \left(\check{m}_{ijk} + b_{14} \check{\varPsi}_{ijk}\right), \end{equation}where
\begin{equation} b_{14} = \tfrac{1}{18}. \end{equation}
Comparison of the coefficients of 
$\check {m}_{ijk}$ in the PBC for rank-3 tensor fluxes (obtained on substituting 
$\mathcal {J}_r$ and 
$\mathcal {F}_s$ from (4.23a,b) in (4.6)) and MBC (4.24) yields
\begin{equation} \mathcal{L}_{rs} = 6 \varkappa. \end{equation}
Thus, the PBC for the rank-3 tensor flux is obtained from (4.6) with 
$\mathcal {J}_r$, 
$\mathcal {F}_s$ and 
$\mathcal {L}_{rs}$ being given by (4.23a,b) and (4.26). On simplification, the PBC for the rank-3 tensor flux is given by (4.24) with
\begin{equation} b_{14} = 0. \end{equation}5. Onsager reciprocity relations in some canonical problems
In this section, we consider two canonical boundary value problems, namely (i) the Poiseuille and thermal transpiration gas flows and (ii) evaporation of a spherical liquid droplet, for which the Onsager reciprocity relations have been shown to hold for the linearised Boltzmann equation (LBE) (Chernyak & Margilevskiy Reference Chernyak and Margilevskiy1989; Wu et al. Reference Wu, Zhang, Liu, Zhang and Reese2017). We shall employ the LR26 equations with the PBC and MBC to these problems to demonstrate that the solutions of the LR26 equations with the MBC for these problems either do not satisfy the Onsager reciprocity relations or do not respect the second law of thermodynamics or do not comply with both. On the other hand, the solutions of the LR26 equations with the PBC satisfy the Onsager reciprocity relations and always respect the second law of thermodynamics. For comparison purpose, we shall also include the results obtained through the LR13 equations with the MBC and PBC.
5.1. Poiseuille flow and thermal transpiration flow
Let us first consider the steady state [
$\partial _t(\cdot )=0$] of a gas confined between two stationary (
$v^{w}_{i}=0$), isothermal (
$T^{w}=1$) and fully diffusive (
$\chi =1$) walls of a channel located at 
$y=\mp 1/2$; see figure 1. The flow is assumed to be fully developed and driven either by a uniform external force 
$F$ (Poiseuille flow) in the positive 
$x$-direction or by a constant temperature gradient parallel to the walls 
$\tau$ (thermal transpiration flow). The normals for the bottom and top walls are pointing in 
$\pm y$-directions, respectively.
Figure 1. Schematic of the Poiseuille flow (
$F=1$, 
$\tau =0$) and thermal transpiration flow (
$F=0$, 
$\tau =1$) in a channel.
The solution for the velocity and heat flux, for the problem under consideration, from the LR26 equations read as
$$\begin{gather} v_{x} = C_{1}+\underline{\frac{F}{2{Kn}}\left( \frac{1}{4}-y^{2}\right)} -\frac{2}{5} \sum_{i=1}^{4}\left( C_{i}^{K}-\frac{5}{2}\beta_{i}\right)\exp \left( \frac{y\gamma _{i}}{{Kn}}\right), \end{gather}$$
$$\begin{gather}q_{x} = \underline{-\frac{5{Kn}}{2\Pr }\tau} -\frac{{Kn}}{\Pr } F+\sum_{i=1}^{4}C_{i}^{K}\exp \left( \frac{y\gamma _{i}}{{Kn}}\right), \end{gather}$$
where 
$C_{1}$ and 
$C^{K}_{1}$, 
$C^{K}_{2}$, 
$C^{K}_{3}$ and 
$C^{K}_{4}$ are the constants of integration. Moreover, due to symmetry about the 
$x$-axis, 
$C^{K}_{2}=C^{K}_{1}$, 
$C^{K}_{4}=C^{K}_{3}$, 
$\beta _2 = \beta _1$ and 
$\beta _4 = \beta _3$. The unspecified three constants of integration, i.e. 
$C_1$, 
$C^{K}_{1}$ and 
$C^{K}_{3}$, are to be evaluated from the three (vector) boundary conditions at 
$y=-1/2$.
The underlined terms in (5.1) represent the classical solution of the NSF equations without velocity slip for the Poiseuille flow. The constant term 
$C_{1}$ in (5.1) accounts for the velocity slip and terms with exponential functions describe the Knudsen layer contributions to the velocity and heat flux. In classical hydrodynamics, the velocity profile for the Poiseuille flow is parabolic; the LR26 equations add velocity slip 
$C_{1}$ and Knudsen layer contributions. Similarly, the Fourier's contribution to the heat flux is given by the underlined term in (5.1b). The coefficients 
$\gamma _{i}$'s (
$i\in \{1,2,3,4\}$) in (5.1) determine the thickness of the Knudsen layer; hence, the LR26 equations describe the Knudsen layer as a superposition of four exponential functions. The coefficients 
$\beta _1$, 
$\beta _3$ and 
$\gamma _{i}$'s (
$i\in \{1,2,3,4\}$) for MM read as 
$\beta _1 = -0.24697$, 
$\beta _3 = 2.25494$, 
$\gamma _1 = -0.510199$, 
$\gamma _2 = 0.510199$, 
$\gamma _3 = -1.26497$ and 
$\gamma _4 = 1.26497$ (Gu et al. Reference Gu, Emerson and Tang2010; Gu & Emerson Reference Gu and Emerson2014).
The total mass flow and heat flow rates of the gas are, respectively, defined as
\begin{equation} M = \int_{{-}1/2}^{1/2}v_{x}\,\mathrm{d}y \quad\text{and}\quad Q = \int_{{-}1/2}^{1/2}q_{x}\,\mathrm{d}y. \end{equation}
We shall use the superscripts ‘
$P$’ and ‘
$T$’ for the Poiseuille (
$F=1$, 
$\tau =0$) and transpiration (
$F=0$, 
$\tau =1$) flows, respectively.
Figure 2(a) exhibits the mass flow rate, 
$M^{P}$, in the Poiseuille flow plotted over the Knudsen number. The results computed from the MBC applied to the LR26 equations with the MM collision model are compared with the LBE data from Wu et al. (Reference Wu, Zhang, Liu, Zhang and Reese2017). The Knudsen minimum is observed in the Poiseuille flow, in which the mass flow rate 
$M^{P}$ first decreases with the Knudsen number, attains a minimum value around a critical Knudsen number and then increases; see figure 2(a). The LR26 equations with PBC closely follow the mass flow rate profile from the LBE up to 
${Kn}\lesssim 5$.
Figure 2. Poiseuille flow and thermal transpiration flow: (a) the total mass flux 
$M^{P}$ in the Poiseuille flow and (b) violation of the Onsager reciprocity relation, 
$M^{T}=Q^{P}$, for the MBC in the case of the LR26 equations (main panel) and in the case of the LR13 equations (inset). The triangles in the left panel denote the results obtained from the LBE in Wu et al. (Reference Wu, Zhang, Liu, Zhang and Reese2017).
Due to the microscopic reversibility, the Onsager reciprocity relations should hold, which in this case state that 
$M^{T}=Q^{P}$ (Wu et al. Reference Wu, Zhang, Liu, Zhang and Reese2017). It is evident from figure 2(b) that the MBC do not comply with the Onsager reciprocity relations for the LR26 (and also for the LR13; shown in the inset) equations, while the PBC do.
Figure 3 illustrates the entropy generation rates at the interface computed with the LR26 equations (main panels) and the LR13 equations (insets) in the case of (a) the Poiseuille flow and (b) the transpiration flow. It is clear from figure 3(a) that the MBC for the LR26 equations (as well as for the LR13 equations) always violate the second law of thermodynamics (i.e. the entropy production 
$\varSigma$ becomes negative) for large Knudsen numbers. For the transpiration flow case (figure 3b), the violation of the second law of thermodynamics by the MBC for the LR26 equations occurs already at a relatively small Knudsen number (
${Kn}<1$), which is within the transition regime. Surprisingly, the LR13 equations with MBC in the transpiration case do not contravene with the second law of thermodynamics whilst they do in the Poiseuille flow case around 
${Kn}=1$; see the insets in figure 3. Nevertheless, the PBC always predict a positive entropy production in both the cases for both models in contrast to the conventional MBC.
Figure 3. The entropy generation rates at the boundary for (a) the Poiseuille flow and (b) the thermal transpiration flow. The main panels and insets exhibit the results computed with the LR26 and LR13 equations, respectively.
5.2. Evaporation of a spherical liquid droplet
We consider a liquid droplet of (fixed) initial radius 
$\hat {R}_0$ with a given interface temperature 
$\hat {T}_l$ and the corresponding saturation pressure 
$\hat {p}_{{sat}}$, immersed in its own vapour; see figure 4. The far-field (dimensional) temperature 
$\hat {T}_{\infty }$, pressure 
$\hat {p}_{\infty }$ and the initial radius of the droplet 
$\hat {R}_0$ are used for the purpose of non-dimensionalisation.
Figure 4. Schematic of an evaporating liquid droplet surrounded by its own vapour.
We consider two cases: (i) the flow is driven by a unit pressure difference while the temperature of the liquid is equal to the far-field temperature [i.e. 
$\varDelta p := (\hat {p}_{{sat}}-\hat {p}_{\infty })/\hat {p}_{\infty }=1$ and 
$\varDelta T := (\hat {T}_{l}-\hat {T}_{\infty })/\hat {T}_{\infty }=0$], and (ii) the flow is driven by a unit temperature difference while the saturation pressure is the same as the far-field pressure (i.e. 
$\varDelta p =0$ and 
$\varDelta T =1$). The problem has been investigated in detail by Rana et al. (Reference Rana, Lockerby and Sprittles2018b). We skip the details here for the sake of succinctness and focus mainly on looking at the Onsager reciprocity relations; for computational (and other) details, the reader is referred to Rana et al. (Reference Rana, Lockerby and Sprittles2018b). The Onsager reciprocity relations for this problem lead to 
$Q^{\varDelta p}=M^{\varDelta T}$ (Chernyak & Margilevskiy Reference Chernyak and Margilevskiy1989), where 
$Q^{\varDelta p}$ is the (dimensionless) heat flux in the pressure-driven process, 
$M^{\varDelta T}$ is the (dimensionless) mass flux in the temperature-driven process, and the superscripts ‘
$\varDelta p$’ and ‘
$\varDelta T$’ are used to denote the pressure- and temperature-driven processes, respectively.
The mass and heat fluxes in the pressure- and temperature-driven processes computed from the R13 and R26 theories using the PBC and MBC are illustrated in figure 5 for different values of the Knudsen number. The top and bottom rows of figure 5 exhibit the results for the pressure- and temperature-driven processes, respectively, while the left and right columns display the mass and heat fluxes, respectively. For comparison purpose, the results obtained by solving the LBE (Sone, Takata & Wakabayashi Reference Sone, Takata and Wakabayashi1994) are also included and they are depicted by symbols in the figure. In general, the results from both the R13 and R26 theories agree with those from the LBE for 
${Kn} \lesssim 0.5$. However, for the mass flux in the pressure-driven case (figure 5a), the R26 equations on using the PBC (blue solid lines) give an excellent match with the LBE results (symbols) for 
${Kn} \lesssim 0.5$ whereas the results from the other theories are comparatively higher than the LBE results even for small Knudsen numbers. Let us focus on panels (b,c) of figure 5, which illustrate the cross-effects, i.e. the heat flux in the pressure-driven case and mass flux in the temperature-driven case, respectively, and account for the verification of the Onsager reciprocity relations. The LBE predicts the same values the heat flux in the pressure-driven case and mass flux in the temperature-driven case, thanks to the Onsager reciprocity relations (Chernyak & Margilevskiy Reference Chernyak and Margilevskiy1989; Sone et al. Reference Sone, Takata and Wakabayashi1994). Evidently from panels (b,c) of figure 5, the heat flux in the pressure-driven case computed with the LR26/LR13 theory using the PBC (blue solid line for the LR26 equations and green solid line for the LR13 equations) seems to be the same as the mass flux in the temperature-driven case computed with the LR26/LR13 theory using the PBC (blue solid line for the LR26 equations and green solid line for the LR13 equations), whereas the heat flux in the pressure-driven case computed with the LR26/LR13 theory using the MBC (red dashed line for the LR26 equations and magenta dashed line for the LR13 equations) seems to be different from the mass flux in the temperature-driven case computed with the LR26/LR13 theory using the MBC (red dashed line for the LR26 equations and magenta dashed line for the LR13 equations). In order to have a clear picture of the Onsager reciprocity relations, the difference of 
$Q^{\varDelta p}$ and 
$M^{\varDelta T}$ is plotted in figure 6, which clearly shows that both the LR13 and LR26 equations on using the MBC for this problem violate the Onsager reciprocity relations while both the LR13 and LR26 equations on using the PBC respect the Onsager reciprocity relations for all Knudsen numbers.
Figure 5. Profiles of the mass and heat fluxes plotted over the Knudsen number for the problem of evaporation of a spherical liquid droplet: (a) mass flux in the pressure-driven case, (b) heat flux in the pressure-driven case, (c) mass flux in the temperature-driven case and (d) heat flux in the temperature-driven case. The LBE data (symbols) are taken from Sone et al. (Reference Sone, Takata and Wakabayashi1994).
Figure 6. Violation of the Onsager reciprocity relations for an evaporating droplet.
We have also computed the entropy generation rate at the interface for this problem and found (but not shown here for the sake of conciseness) that it remains positive for the R13 and R26 equations with the PBC in both pressure- and temperature-driven cases, also remains positive for the R13 equations with the MBC in the pressure-driven case but does not remain positive for the R13 equations with the MBC in the temperature-driven case for all Knudsen numbers. The R26 equations with the MBC yield a positive entropy generation rate for this problem.
6. Gas flow past an evaporating droplet: a new analytic solution
As an application of the above findings, we compute the analytic solution of the LR26 equations for a steady vapour passing over an evaporating spherical droplet. Furthermore, we also consider a special case of this problem – when there is no evaporation; this is essentially the problem of a streaming flow of a gas past a rigid sphere, which is a classical problem in fluid dynamics. The results presented in this section are determined using the PBC only, since we have established that they are the physically admissible boundary conditions. Nevertheless, we show in the supplementary material available at https://doi.org/10.1017/jfm.2021.622 that the difference between the results obtained using the PBC and MBC for this problem is negligibly small.
We consider a rarefied gas flow past a spherical droplet as illustrated in figure 7. It is convenient to study the problem in the spherical coordinate system 
$(r,\theta ,\phi )$ that relates to the Cartesian coordinate system 
$(x,y,z)$ via 
$(x,y,z) \equiv (r \sin {\theta } \cos {\phi }, r \sin {\theta } \sin {\phi }, r \cos {\theta })$, where 
$r\in [0,\infty )$, 
$\theta \in [0,{\rm \pi} ]$ and 
$\phi \in [0,2{\rm \pi} )$. The flow is assumed to be in the 
$(r,\theta )$-plane coming from the left as shown in the figure. Owing to spherical symmetry, the flow is independent of the direction 
$\phi$. Consequently, all the field variables can depend only on 
$r$ and 
$\theta$. All the vectors and tensors involved in the problem need to be converted to spherical coordinates. To this end, it is worthwhile noting that a symmetric-tracefree 
$k$-rank tensor has 
$2k+1$ independent components in three dimensions (Gupta Reference Gupta2020).
Figure 7. Schematic of a flow past a sphere.
For the considered problem, any vector 
$a_i \in \{v_i,q_i,\varOmega _i\}$ in the spherical coordinate system is given by
\begin{equation} a_i=(a_r, a_\theta, 0)^{\textrm{T}}. \end{equation}
A symmetric-tracefree second-rank tensor 
$a_{ij} \in \{\sigma _{ij}, R_{ij}\}$ in the spherical coordinate system is given by
\begin{equation} a_{ij}=\begin{bmatrix} a_{rr} & a_{r\theta} & 0 \\ a_{r\theta} & - \frac{1}{2} a_{rr} & 0 \\ 0 & 0 & - \frac{1}{2} a_{rr} \end{bmatrix}. \end{equation}
Note that a symmetric-tracefree second-rank tensor has five independent components; however, for the problem under consideration, it is clear from (6.2) that all components of 
$a_{ij}$ can be expressed in terms of only two components, namely 
$a_{rr}$ and 
$a_{r\theta }$. A symmetric-tracefree third-rank tensor has seven independent components; however, for the problem under consideration, all components of 
$a_{ijk}$ can again be expressed in terms of 
$a_{rrr}$ and 
$a_{rr\theta }$. The five independent components (other than 
$a_{rrr}$ and 
$a_{rr\theta }$) of 
$a_{ijk}$, for the problem under consideration, read as
\begin{equation} a_{r\theta\theta} ={-} \tfrac{1}{2} a_{rrr}, \quad a_{\theta\theta\theta} ={-} \tfrac{3}{4} a_{rr\theta}, \quad a_{rr\phi} = a_{\theta\theta\phi} = a_{r\theta\phi} = 0. \end{equation}
A symmetric-tracefree fourth-rank tensor 
$a_{ijkl}$ has nine independent components; however, for the problem under consideration, all components of 
$a_{ijkl}$ can again be expressed in terms of 
$a_{rrrr}$ and 
$a_{rrr\theta }$. The seven independent components (other than 
$a_{rrrr}$ and 
$a_{rrr\theta }$) of 
$a_{ijkl}$, for the problem under consideration, read as
\begin{equation} \left. \begin{gathered} a_{rr\theta\theta} ={-} \tfrac{1}{2} a_{rrrr}, \quad a_{\theta\theta\theta\theta} = \tfrac{3}{8} a_{rrrr}, \quad a_{r\theta\theta\theta} ={-} \tfrac{3}{4} a_{rrr\theta}, \\ a_{rrr\phi} = a_{rr\theta\phi} = a_{r\theta\theta\phi} = a_{\theta\theta\theta\phi} =0. \end{gathered} \right\} \end{equation}Accordingly, the conservation laws (2.3)–(2.5) for this problem read as
$$\begin{gather} \frac{\partial v_r}{\partial r} +\frac{2 v_r}{r} +\frac{\mathscr{D} v_\theta}{r} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial p}{\partial r} +\frac{\partial \sigma_{rr}}{\partial r} +\frac{3 \sigma_{rr}}{r} +\frac{\mathscr{D} \sigma_{r\theta}}{r} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \sigma_{r\theta}}{\partial r} +\frac{3 \sigma_{r\theta}}{r} -\frac{1}{2 r} \frac{\partial \sigma_{rr}}{\partial \theta} +\frac{1}{r} \frac{\partial p}{\partial \theta} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial q_r}{\partial r} +\frac{2 q_r}{r} +\frac{\mathscr{D} q_\theta}{r} = 0, \end{gather}$$
where 
$\mathscr {D} \equiv \cot {\theta } + \frac {\partial }{\partial \theta }$. The stress balance equation (2.6) reduces to
$$\begin{gather} \frac{\partial m_{rrr}}{\partial r} +\frac{4 m_{rrr}}{r} +\frac{4}{5} \frac{\partial q_r}{\partial r} +2 \frac{\partial v_r}{\partial r} +\frac{\mathscr{D} m_{rr\theta}}{r} ={-}\frac{1}{{Kn}} \sigma_{rr},\end{gather}$$
$$\begin{gather} \frac{\partial m_{rr\theta}}{\partial r} +\frac{4 m_{rr\theta}}{r} +\frac{2}{5} \left(\frac{\partial q_\theta}{\partial r} -\frac{q_\theta}{r} \right) +\frac{\partial v_\theta}{\partial r} -\frac{v_\theta}{r} \nonumber\\ -\frac{1}{2 r} \frac{\partial m_{rrr}}{\partial \theta} +\frac{1}{r} \frac{\partial v_r}{\partial \theta} +\frac{2}{5} \frac{1}{r} \frac{\partial q_r}{\partial \theta} ={-} \frac{1}{{Kn}} \sigma_{r\theta}. \end{gather}$$The heat flux balance equation (2.7) gives
\begin{gather} \frac{1}{2} \left(\frac{\partial R_{rr}}{\partial r} +\frac{\mathscr{D} R_{r\theta}}{r} +\frac{3 R_{rr}}{r}\right) +\frac{1}{6} \frac{\partial \varDelta}{\partial r} -\frac{\partial p}{\partial r} +\frac{5}{2} \frac{\partial T}{\partial r} ={-} \frac{{Pr}}{{Kn}} q_r, \end{gather}
\begin{gather} \frac{1}{2} \left(\frac{\partial R_{r\theta}}{\partial r} +\frac{3 R_{r\theta}}{r}\right) +\frac{1}{6} \frac{1}{r} \frac{\partial \varDelta}{\partial \theta} -\frac{1}{4 r} \frac{\partial R_{rr}}{\partial \theta} -\frac{1}{r}\frac{\partial p}{\partial \theta} +\frac{5}{2} \frac{1}{r} \frac{\partial T}{\partial \theta} ={-} \frac{{Pr}}{{Kn}} q_\theta. \end{gather}
The 
$m_{ijk}$ balance equation (2.11) simplifies to
$$\begin{gather} \frac{1}{r}\mathscr{D}
\left[-\frac{6}{5}\sigma_{r\theta}+\varPhi_{rrr\theta} -
\frac{6}{35}R_{r\theta}\right] +\frac{9}{5}\left(
\frac{\partial }{\partial r} -\frac{2}{r}\right)
\sigma_{rr} \nonumber\\ +\left( \frac{\partial }{\partial
r}+\frac{5}{r}\right) \varPhi_{rrrr} +\frac{9}{35}\left(
\frac{\partial }{\partial r}-\frac{2}{r}\right) R_{rr}
={-}\frac{{Pr}_{M}}{{Kn}}m_{rrr},\end{gather}$$
$$\begin{gather}
\frac{6}{5r}\frac{\partial \sigma_{rr}}{\partial \theta}
+\frac{6}{35r} \frac{\partial R_{rr}}{\partial \theta
}-\frac{1}{2r}\frac{\partial \varPhi_{rrrr}}{\partial
\theta} +\frac{8}{5}\left( \frac{\partial }{\partial
r}-\frac{2}{r}\right) \sigma_{r\theta} \nonumber\\
+\frac{8}{35} \left(
\frac{\partial }{\partial r}-\frac{2}{r}\right)
R_{r\theta}+\left( \frac{\partial }{\partial
r}+\frac{5}{r}\right) \varPhi_{rrr\theta
}={-}\frac{{Pr}_{M}}{{Kn}}m_{rr\theta}.
\end{gather}$$
The 
$R_{ij}$ balance equation (2.12) reduces to
$$\begin{gather} \frac{1}{r}\mathscr{D}
\left[ 2m_{rr\theta} -\frac{2}{15}\varOmega_{\theta
}+\psi_{rr\theta}-\frac{28}{15}q_\theta \right]
+\frac{56}{15}\left( \frac{\partial }{\partial
r}-\frac{1}{r}\right) q_r +2\left( \frac{\partial
}{\partial r}+\frac{4}{r}\right) m_{rrr} \nonumber\\
+\left( \frac{\partial }{\partial r}+\frac{4}{r}\right)
\psi_{rrr}+ \frac{4}{15}\left( \frac{\partial }{\partial
r}-\frac{1}{r}\right)
\varOmega_{r}={-}\frac{{Pr}_{R}}{{Kn}}R_{rr},\end{gather}$$
$$\begin{gather}
2\left(
\frac{\partial }{\partial r}+\frac{4}{r}\right)
m_{rr\theta} +\left( \frac{\partial }{\partial
r}+\frac{4}{r}\right) \psi_{rr\theta}+\frac{1}{5}\left(
\frac{\partial }{\partial r}-\frac{1}{r} \right) \varOmega
_{\theta} + \frac{14}{5}\left( \frac{\partial }{\partial
r}- \frac{1}{r}\right) q_\theta \nonumber\\
-\frac{1}{r}\frac{\partial
m_{rrr}}{\partial \theta }+\frac{14}{5r} \frac{\partial q_r
}{\partial \theta }-\frac{1}{2r}\frac{\partial
\psi_{rrr}}{\partial \theta }+\frac{2}{5r}\frac{\partial
\varOmega_{r}}{\partial \theta
}={-}\frac{{Pr}_{R}}{{Kn}}R_{r\theta }.
\end{gather}$$
The 
$\varDelta$ balance equation (2.13) gives
\begin{equation} \frac{1}{r}\mathscr{D} \left( 8 q_{\theta}+\varOmega _{\theta }\right) +\left( \frac{\partial }{\partial r}+\frac{2}{r}\right) \left( 8 q_{r}+\varOmega_{r}\right) ={-}\frac{{Pr}_{\varDelta }}{{Kn}}\varDelta. \end{equation}Finally, the LR26 closure (2.18) simplifies to
$$\begin{gather} \frac{4}{7}\left( \frac{\partial }{\partial r}-\frac{3}{r}\right) m_{rrr} -\frac{3}{7}\frac{\mathscr{D} m_{rr\theta}}{r} ={-}\frac{{Pr}_{\varPhi}}{4{Kn}}\varPhi_{rrrr} , \end{gather}$$
$$\begin{gather}\frac{15}{28}\left( \frac{\partial }{\partial r}-\frac{3}{r}\right) m_{rr\theta} +\frac{5}{14} \frac{1}{r} \frac{\partial m_{rrr}}{\partial \theta} ={-}\frac{{Pr}_{\varPhi}}{4{Kn}} \varPhi_{rrr\theta } , \end{gather}$$
$$\begin{gather}\frac{3}{5}\left( \frac{\partial }{\partial r}-\frac{2}{r}\right) R_{rr} -\frac{2}{5} \frac{\mathscr{D} R_{r\theta}}{r} ={-}\frac{7{Pr}_{\psi}}{27{Kn}}\psi_{rrr} , \end{gather}$$
$$\begin{gather}\frac{8}{15}\left( \frac{\partial }{\partial r}-\frac{2}{r}\right) R_{r\theta }+\frac{2}{5r}\frac{\partial R_{rr}}{\partial \theta } ={-}\frac{7{Pr}_{\psi}}{27{Kn}}\psi_{rr\theta} , \end{gather}$$and
$$\begin{gather} \frac{\partial \varDelta }{\partial r}+\frac{12}{7r}\mathscr{D} R_{r\theta}+\frac{12}{7} \left( \frac{\partial }{\partial r}+\frac{3}{r}\right) R_{rr} =-\frac{3 \textit{Pr}_\varOmega}{7 \textit{Kn}} \varOmega_r, \end{gather}$$
$$\begin{gather}\frac{1}{r}\frac{\partial \varDelta}{\partial \theta }+\frac{12}{7}\left( \frac{\partial R_{r\theta}}{\partial r}+\frac{3R_{r\theta }}{r}\right) -\frac{6}{7r}\frac{\partial R_{rr}}{\partial \theta } =-\frac{3 \textit{Pr}_\varOmega}{7 \textit{Kn}} \varOmega_\theta. \end{gather}$$The mathematical structure of the LR26 equation is the same as the LR13 equations, with which Torrilhon (Reference Torrilhon2010) obtained an analytic solution for the problem of gas flow past a rigid sphere. Following Torrilhon (Reference Torrilhon2010), we shall adopt the same solution method to solve (6.5)–(6.12). In the next subsection, we shall outline the steps involved in solving the LR26 equations for the considered problem. For more details, the reader is referred to Torrilhon (Reference Torrilhon2010).
6.1. Analytic solution of the LR26 equations
The analytic solutions of the LR26 equations (6.5)–(6.12) for the aforesaid problem can be obtained by converting these partial differential equations to ordinary differential equations via assuming an explicit dependency on the azimuthal angle 
$\theta$. To this end, the vectorial and tensorial components with an odd number of indices in 
$\theta$ are taken to be proportional to 
$\sin {\theta }$ while scalars and tensorial components with an even number of indices in 
$\theta$ are taken to be proportional to 
$\cos {\theta }$ (Torrilhon Reference Torrilhon2010). The components for vector fields, namely the velocity 
$\boldsymbol {v}(r,\theta )$ and the heat flux 
$\boldsymbol {q}(r,\theta )$, are taken to be
\begin{equation} \left. \begin{aligned} v_r(r, \theta) & = u_{\infty} [1+a_{1}(r)] \cos{\theta}, \quad v_{\theta}(r, \theta) ={-}u_{\infty}[1+a_{2}(r)] \sin{\theta}, \\ q_r(r, \theta) & = \alpha_1(r) \cos{\theta}, \quad q_\theta(r, \theta) = \alpha_2(r) \sin{\theta}. \end{aligned} \right\} \end{equation}
The scalars, namely 
$T(r,\theta )$, 
$p(r,\theta )$ and 
$\varDelta (r,\theta )$, are taken to be
\begin{equation} T(r,\theta) = c(r) \cos{\theta}, \quad p(r,\theta) = d(r) \cos{\theta}, \quad \varDelta(r,\theta) = \delta(r)\cos{\theta}. \end{equation}The components of the second- and third-rank tensors are taken to be
\begin{equation} \left. \begin{aligned} \sigma_{rr}(r, \theta) & = b_1(r)\cos{\theta}, \quad \sigma_{r\theta}(r, \theta) = b_2(r)\sin{\theta}, \\ R_{rr}(r, \theta) & = \beta_1(r)\cos{\theta}, \quad R_{r\theta}(r, \theta) = \beta_2(r) \sin{\theta}, \\ m_{rrr}(r, \theta) & = \kappa_1(r)\cos{\theta}, \quad m_{rr\theta}(r, \theta) = \kappa_2(r) \sin{\theta}. \end{aligned} \right\} \end{equation}
In (6.13)–(6.15), 
$a_{1}(r)$, 
$a_{2}(r)$, 
$\alpha _{1}(r)$, 
$\alpha _{2}(r)$, 
$c(r)$, 
$d(r)$, 
$\delta (r)$, 
$b_{1}(r)$, 
$b_{2}(r)$, 
$\beta _{1}(r)$, 
$\beta _{2}(r)$, 
$\kappa _{1}(r)$ and 
$\kappa _{2}(r)$ are thirteen unknowns in total, which need to be determined from thirteen equations (6.5)–(6.10).
Substitution of constitutive relations (6.11)–(6.12) along with solution ansatz (6.13)–(6.15) in (6.5)–(6.10) results into a set of eight first-order ordinary differential equations (emanating from (6.5)–(6.7)) and five second-order ordinary differential equations (from (6.8)–(6.10)). The latter include second derivatives of 
$\delta (r)$, 
$\beta _1(r)$, 
$\beta _2(r)$, 
$\kappa _1(r)$ and 
$\kappa _2(r)$ due to regularised constitutive relations (6.11)–(6.12). With this, we obtain an analytic solution for the unknowns 
$a_{1}(r)$, 
$a_{2}(r)$, 
$\alpha _{1}(r)$, 
$\alpha _{2}(r)$, 
$c(r)$, 
$d(r)$, 
$\delta (r)$, 
$b_{1}(r)$, 
$b_{2}(r)$, 
$\beta _{1}(r)$, 
$\beta _{2}(r)$, 
$\kappa _{1}(r)$ and 
$\kappa _{2}(r)$, which has been relegated to Appendix E for better readability. Expressions (E1)–(E13) contain eight integration constants 
$C_{1,2,3}$ and 
$K_{1,2,3,4, 5}$ that are determined from boundary conditions (4.8)–(4.14).
Let us investigate a few qualitative features of these solutions. These solutions are the sums of the bulk contributions (which are polynomials in 
$1/r$) and the Knudsen layer contributions (terms containing the exponential functions 
$\varTheta _i(r)$); see Appendix E. The bulk solutions contain three constants of integration 
$C_{1,2,3}$ while Knudsen layers introduce five additional constants – 
$K_{1,2}$ for vector variables and 
$K_{3,4,5}$ for scalar variables. The bulk part can further be divided into the classical NSF contribution and second-order contribution (in the Knudsen number) that is neglected in the classical NSF theory. In the hydrodynamic limit (i.e. Stokes flow without velocity slip and temperature jump), 
$C_1=-3$ and 
$C_2=1/2$ while 
$C_3=0$ (Lamb Reference Lamb1932). The last relation (
$C_3=0$) implies an isothermal flow condition. Hence, the non-zero bulk contributions to the temperature and heat flux stem from the second-order terms, i.e. the terms with 
$C_3$ in (E3a) and (E1c,d). The NSF equations with a velocity-slip boundary condition predict 
$C_1$ and 
$C_2$ as a function of the Knudsen number 
${Kn}$ and the accommodation coefficient 
$\chi$ (Torrilhon Reference Torrilhon2010).
For the Knudsen layer contributions in (E1)–(E13), the LR13 theory (Torrilhon Reference Torrilhon2010) approximates Knudsen layers with two exponential functions whilst the LR26 theory gives superposition of five exponential functions to describe the Knudsen layers. On the contrary, the NSF theory does not predict Knudsen layers at all. In the following, we shall show that the LR26 theory gives a remarkably improved flow description with added details to the Knudsen layers.
6.2. Gas flow past a rigid sphere
First, let us consider the classical problem of a gas flow past a rigid sphere. The same problem was investigated by Torrilhon (Reference Torrilhon2010) using the LR13 and linear NSF equations with velocity slip.
Figure 8 illustrates the profiles of the radial and polar components of the velocity, pressure and temperature as functions of the radial distance 
$r$ starting from the surface of the sphere 
$r=1$. The exact solutions for the problem under consideration obtained with the LR26 equations in the present work are compared with those obtained with the LR13 and linear NSF solutions in Torrilhon (Reference Torrilhon2010) and with the numerical solutions of the LBE for three values of the Knudsen number 
${\mathsf{Kn}} = 0.1, 0.4, 1$ given in Takata et al. (Reference Takata, Sone and Aoki1993). Note that the Knudsen number 
${\mathsf{Kn}}$ used in Takata et al. (Reference Takata, Sone and Aoki1993) relates to the Knudsen number 
${Kn}$ used in the present work via 
${Kn}= 1.27\,{\mathsf{Kn}}/\sqrt {2}$. Therefore, the Knudsen number values 
${\mathsf{Kn}} = 0.1, 0.4, 1$ in Takata et al. (Reference Takata, Sone and Aoki1993) correspond to the Knudsen number values 
${Kn} \approx 0.09, 0.36, 0.9$, respectively, used in the present work. As expected, for small values of the Knudsen number [red curves and symbols (
${Kn} = 0.09$)], all theories, including the NSF theory, give a good match with the LBE data for the velocities. However, for moderate values of the Knudsen number [blue curves and symbols (
${Kn} = 0.36$)], the results from the NSF theory show deviations from those from kinetic theory (LBE results) while the LR13 theory still gives reasonably good results. The mismatch between the solutions from the NSF equations and LBE is attributed to the lack of Knudsen layer description by the former. For even larger values of the Knudsen number [green curves and symbols (
${Kn} = 0.9$)], even the results from the LR13 theory start deviating from the LBE results. Notably, the LR13 theory in this case predicts 
$u_\theta / (u_\infty \sin {\theta }) >1$ near the surface of the sphere, which means that the slip velocity on the top of the sphere (
$\theta ={\rm \pi} /2$) is larger than the far-field velocity – a result which is clearly unphysical and contradicted by more accurate theories (LBE and LR26). Hence, one may conclude that for large values of the Knudsen number, the results from LR13 theory cannot be trusted. On the other hand, the LR26 theory offers an excellent match upto 
${Kn} \sim 1$. Interestingly, the LR26 theory also predicts 
$u_\theta / (u_\infty \sin {\theta }) > 1$ for 
${Kn}\gtrsim 5$, which deduces that the LR26 theory can predict reasonably good results for the Knudsen numbers 
$0<{Kn}\lesssim 5$.
Figure 8. Profiles of the (a) radial velocity, (b) polar velocity (c) pressure and (d) temperature of the gas past a rigid sphere as functions of the radial distance from the surface of the sphere for various Knudsen numbers. Results obtained from the LR26 theory (solid lines) are compared with those obtained from the LR13 (dashed lines) and NSF (dotted lines) theories (Torrilhon Reference Torrilhon2010) and with the numerical solution of the LBE (symbols) from Takata, Sone & Aoki (Reference Takata, Sone and Aoki1993).
Another remarkable feature in figure 8(d) is that the NSF theory predicts a non-negative temperature profile whereas all other theories, including the LBE, give a non-positive temperature profile. To demonstrate the consequences of this, the heat flux lines and temperature contours obtained from LR26, LR13 and NSF theories are illustrated in figure 9. It is evident from the figure that all the theories, including the NSF theory, do predict the thermal polarisation – an effect that, for the problem under consideration, shows different temperatures in the front and back sides of the sphere in the absence of nonlinear dissipations and external temperature differences (Aoki & Sone Reference Aoki and Sone1987; Takata et al. Reference Takata, Sone and Aoki1993; Torrilhon Reference Torrilhon2010). Nevertheless, the LR13 and LR26 theories show a hot region in the front side of the sphere and a cold region in the back side of the sphere while the temperature prediction by the NSF theory is the other way round. The temperature predictions by the LR13 and LR26 theories are in agreement with the findings of Takata et al. (Reference Takata, Sone and Aoki1993) based on the LBE and of Torrilhon (Reference Torrilhon2010) based on the LR13 equations, and also make sense because the temperature is likely to be more on the front side of the sphere due to compression of the gas at the front side of the sphere and is likely to be comparatively lesser on the back side of the sphere due to expansion of the gas at the back side of the sphere. Furthermore, the LR13 and LR26 theories predict that the heat flows from the cold region (back side of the sphere) to the hot region (front side of the sphere), which is a non-Fourier effect. The NSF theory cannot capture this effect and shows the heat flowing from the hot region to the cold region. Such a non-Fourier effect is well established in lid-driven cavity problems, (see, e.g. Rana, Torrilhon & Struchtrup Reference Rana, Torrilhon and Struchtrup2013, and references therein). By comparing the top and bottom rows of figure 9, it is also clear that the non-uniformity in temperature (i.e. the temperature difference between the temperatures at the front and back sides of the sphere) is more for a larger Knudsen number. This is also in agreement with the findings of Takata et al. (Reference Takata, Sone and Aoki1993).
Figure 9. Heat flux lines drawn on top of the temperature contours for 
${Kn} = 0.09$ (a–c) and 
${Kn} = 0.36$ (d–f) from the LR26, LR13 and NSF theories.
Figure 10 illustrates the drag force normalised with the Stokes drag, 
$F_{{Stokes}} = 6{\rm \pi} {Kn}u_{\infty }$, plotted over the Knudsen number. The experimental data from Goldberg (Reference Goldberg1954) (fitted by Bailey et al. Reference Bailey, Barber, Emerson, Lockerby and Reese2005) is included for comparison. As expected, all theories agree with the experimental data for small values of the Knudsen number. As the Knudsen number increases, the drag decreases due to slip, which in not embedded in the Stokes formula and, hence, the Stokes formula gives vanishing drag force for all values of the Knudsen number. Notably, the drag force from the NSF theory with velocity-slip boundary conditions attains a non-zero positive value as 
${Kn} \to \infty$ while the other theories (LR13 and LR26) predict a vanishing drag in this limit – an observation tantamount to experimental results. Furthermore, the LR26 theory provides an excellent prediction for the drag force for values of the Knudsen number beyond unity (upto 
${Kn} \approx 10$). On the contrary, the LR13 theory underpredicts the drag force significantly for 
${Kn} \gtrsim 0.2$.
Figure 10. Normalised drag force from the NSF (dot-dashed red line), R13 (dashed green line) and R26 (solid blue line) theories as a function of the Knudsen number. The circles denote the experimental data of Goldberg (Reference Goldberg1954) (fitted by Bailey et al. Reference Bailey, Barber, Emerson, Lockerby and Reese2005). The results are normalised with the Stokes drag 
$F_{{Stokes}}=6{\rm \pi} {Kn} u_\infty$.
6.3. Gas flow past an evaporating spherical droplet
Finally, we consider the case of an evaporating droplet placed in a saturated vapour (
$p_{{sat}} = p_\infty = 0$) flowing past the droplet with velocity 
$u_\infty$, i.e. phase change is driven by the flow field, not by the far field. The assumptions made in the previous section remain true here, i.e. the shape and size of the droplet remain fixed and the internal motion of the droplet is ignored. In addition, the evaporation coefficient 
$\vartheta$ is assumed to be unity.
Figures 11 again illustrates the profiles of the velocity components, pressure and temperature as functions of the radial direction 
$r$ starting from the surface of the sphere 
$r=1$. The exact solutions for the problem under consideration obtained with the LR26 equations (solid lines) and LR13 equations (dashed lines) are compared with those obtained with the numerical solutions of the LBE for three values of the Knudsen number 
${Kn} = 0.1, 0.4, 1$ given in Sone et al. (Reference Sone, Takata and Wakabayashi1994). Notice that the Knudsen number 
${Kn}$ is the same as that used in Takata et al. (Reference Takata, Sone and Aoki1993). Therefore, similarly to the above, these Knudsen number values correspond to 
${Kn} = 0.09, 0.26, 0.9$, respectively.
Figure 11. Profiles of the (a) radial velocity, (b) polar velocity (c) pressure and (d) temperature of the gas past an evaporating spherical droplet as functions of the radial distance from the surface of the droplet for various Knudsen numbers. Results obtained from the LR26 theory (solid lines) are compared with those obtained from the LR13 theory (dashed lines) and with the numerical solution of the LBE (symbols) from Sone et al. (Reference Sone, Takata and Wakabayashi1994).
As expected, when the Knudsen number is small [red lines and symbols (
${Kn} = 0.09$)], both the LR13 and LR26 theories (dashed and solid lines) agree well with the results from the LBE (symbols). However, as the Knudsen number increases, the differences between the results from the two theories become noticeable – in particular, for the radial velocity and temperature. At higher Knudsen numbers (
${Kn}\gtrsim 0.5$), the Knudsen layer description by the LR13 theory is inaccurate whilst the LR26 theory remains valid even for 
${Kn}= 1$.
The normalised drag force in the case of a gas flow past an evaporating droplet is illustrated in figure 12. The solutions obtained from the NSF, LR13 and LR26 theories are compared against the numerical solutions of the LBE given in Sone et al. (Reference Sone, Takata and Wakabayashi1994). In contrast to the non-evaporating droplet case, all theories in this case predict vanishing drag as 
${Kn}\to \infty$. For large values of the Knudsen number, the NSF theory overpredicts the drag while the LR13 theory underpredicts it. Again, the LR26 theory gives accurate predictions for the drag force even for large values of the Knudsen number (
${Kn} \lesssim 10$).
Figure 12. Normalised drag force from the LR13 (dashed green line) and LR26 (solid blue line) theories and from the numerical solution of the LBE (circles) from Sone et al. (Reference Sone, Takata and Wakabayashi1994) as a function of the Knudsen number. The results are normalised with the Stokes drag 
$F_{{Stokes}}=6{\rm \pi} {Kn} u_\infty$.
7. Conclusion
In this paper, we have addressed the question of thermodynamically admissible boundary conditions for the LR26 equations. It has been shown that these equations are accompanied with a proper entropy inequality with non-negative entropy generation (the 
$H$-theorem).
The entropy production rate at the interface has been written as a sum of the products of the fluxes with odd degree in 
$n$ and the requirement of the non-negative entropy generation at the wall has been imposed, which have eventually yielded the desired boundary conditions for the LR26 equations. Furthermore, considering the coupling between the cross-terms of the same tensorial structure (the Curie principle) a positive entropy production has been achieved by setting the unknown boundary values for the odd fluxes proportional to their driving forces. The phenomenological coefficients appearing in these boundary conditions have been obtained by comparing them with those obtained via the Maxwell accommodation model (with Grad's approximation) to the leading order, thus allowing us to formulate the PBC aligned with the Maxwell kinetic-wall model. We have demonstrated that the solutions of the LR26 equations with the MBC either do not satisfy the Onsager reciprocity relations or do not respect the second law of thermodynamics or do not comply with both. On the other hand, the solutions of the LR26 equations with the PBC satisfy the Onsager reciprocity relations and always respect the second law of thermodynamics.
As an application of the PBC, we have computed the analytic solution of the LR26 equations for a steady vapour passing over an evaporating spherical droplet. To the best of the authors’ knowledge, this is the first time this problem has been considered within the extended fluid dynamics framework. These analytic solution have been found to be useful in understanding some intriguing rarefaction effects, such as the Knudsen layers and thermal polarisation. A special case of this problem, when there is no evaporation, has also been considered, which is a classic problem in fluid dynamics. It has been found that the LR26 equations with the PBC are able to provide a more detailed and accurate description of the flow features emanating from rarefaction effects up to 
${Kn}\lesssim 1$. Furthermore, the analytic results for the drag force obtained from the LR26 theory with the PBC have been found to be in excellent agreement with the theory/experiments even for very large values of the Knudsen number – in the case of both an evaporating droplet and a non-evaporating droplet.
In § 4, it was noticed that matching the coefficients of the lower-order moments (i.e. moments involved in the G26 equations) in the PBC and MBC yields the matrices 
$\mathcal {L}$ to be symmetric and semi-positive definite. The reason behind it is presently unclear and will be explored in future work. Furthermore, it is worthwhile noting that the results in the present work are valid only for slow moving droplets, i.e. when the Reynolds number 
${Re}$ (
$={Ma}\sqrt {\gamma }/{Kn})$ is much less than unity. Here 
$\gamma$ (=5/3 for monatomic gases) is the ratio of specific heat capacities and 
${Ma}$ is the Mach number. For a relatively large value of the Reynolds number, the nonlinear convective terms play an important role in evaporation dynamics, which have been neglected in the present work. In classical fluid dynamics texts (see, e.g. Lamb Reference Lamb1932), Oseen correction terms are utilised to obtain solutions for the flow past various bodies at moderate Reynolds numbers [
$O({Re}) \sim 1$]. Oseen correction terms are linear in nature as oppose to the nonlinear Navier–Stokes operator, making these equations analytically tractable. The present work can be extended to flows with a small but finite Reynolds number by including the Oseen correction terms (with the partial inclusion of convective terms) in the future to study the combined effects of inertia and rarefaction on evaporation dynamics. Motivated by the present findings, the present work will also be extended to investigate unsteady flow problems and to incorporate coupling with liquid dynamics (see, e.g. Rana, Lockerby & Sprittles Reference Rana, Lockerby and Sprittles2019) elsewhere in the future. Thermophoresis on a sphere is a classic problem in rarefied gas flows, which has been studied by many researchers; see, e.g. Yamamoto & Ishihara (Reference Yamamoto and Ishihara1988), Takata, Aoki & Sone (Reference Takata, Aoki and Sone1994), Takata & Sone (Reference Takata and Sone1995), Beresnev & Chernyak (Reference Beresnev and Chernyak1995), Takata (Reference Takata2009), Padrino, Sprittles & Lockerby (Reference Padrino, Sprittles and Lockerby2019) and Kalempa & Sharipov (Reference Kalempa and Sharipov2020). The model developed in the present work can also be used to investigate thermophoresis on an evaporating droplet in the future.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2021.622.
Funding
A.S.R. acknowledges the funding from the BITS, Pilani (India) through the Research Initiation Grant. V.K.G. gratefully acknowledges financial support through the ‘MATRICS’ project MTR/2017/000693 funded by the SERB, India. This work has also been financially supported in the UK by EPSRC grants no. EP/N016602/1, EP/P031684/1, EP/S029966/1.
Declaration of interests
The authors report no conflict of interest.
Appendix A. The entropy density for the LR26 equations from the Boltzmann entropy density
The Boltzmann entropy density 
$\hat {\eta }$ and entropy flux 
$\hat {\varGamma }_i$ (in dimensional form) are given by (3.2) and (3.3). We replace 
$\hat {f}$ in (3.2) and (3.3) with the G26 distribution function (Struchtrup Reference Struchtrup2005; Gupta & Torrilhon Reference Gupta and Torrilhon2012; Gupta Reference Gupta2015)
\begin{equation} \hat{f}_{{G26}}=\hat{f}_{{M}} \left( 1+\varPhi \right) , \end{equation}where
\begin{equation} \hat{f}_{{M}}=\frac{\hat{\rho}}{(2{\rm \pi} \hat{\theta})^{3/2}}\exp \left( -\frac{\hat{C}^{2}}{2\hat{\theta}}\right), \end{equation}is the Maxwellian distribution function and
$$\begin{align} \varPhi
&=\frac{\hat{\sigma}_{ij}}{2}\frac{\hat{C}_{i}\hat{C}_{j}}{\hat{\rho}
\hat{\theta}^{2}}-\frac{\hat{q}_{i}\hat{C}_{i}}{\hat{\rho}\hat{\theta}^{2}}
\left( 1-\frac{\hat{C}^{2}}{5\hat{\theta}}\right)
+\frac{\hat{m}_{ijk}}{6
\hat{\rho}\hat{\theta}^{3}}\hat{C}_{i}\hat{C}_{j}\hat{C}_{k}\nonumber\\
&\quad +\frac{\hat{\varDelta}}{8\hat{\rho}\hat{\theta}^{2}}\left(
1-\frac{2\hat{C}^{2}}{3\hat{
\theta}}+\frac{\hat{C}^{4}}{15\hat{\theta}^{2}}\right)
-\frac{\hat{R}_{ij}}{4\hat{\rho}\hat{\theta}^{3}}\left(
1-\frac{\hat{C}^{2} }{7\hat{\theta}}\right)
\hat{C}_{i}\hat{C}_{j} \end{align}$$
is the dimensionless deviation of the distribution function from the equilibrium state. Here, 
$\hat {C}_{i}=\hat {c}_{i} -\hat {v}_{i}$ is the (dimensional) peculiar velocity with 
$\hat {c}_{i}$ and 
$\hat {v}_{i}$ being the (dimensional) microscopic and macroscopic velocities, respectively.
For the entropy density and entropy flux of the LR26 equations, the integrations in (3.2) and (3.3) need to be evaluated by expending 
$\hat {f}_{{G26}}$ up to second order in deviations. For this purpose, we introduce a smallness parameter 
$\varepsilon$, whose power in a term tells the order of the deviation. The parameter is introduced just for the bookkeeping purpose and will be replaced with unity at the end of computations. Each moment is replaced with its value in the equilibrium plus its deviation from its equilibrium value. We write
\begin{equation} \left. \begin{gathered} \hat{\rho} := \hat{\rho}_0 (1 + \varepsilon \rho), \quad \hat{v}_i := \varepsilon \sqrt{\hat{\theta}_0}\, v_i, \quad \hat{\theta} := \hat{\theta}_0 (1 + \varepsilon T), \quad \hat{\sigma}_{ij} := \varepsilon \hat{\rho}_0 \hat{\theta}_0 \sigma_{ij}, \\ \hat{q}_i := \varepsilon \hat{\rho}_0 \hat{\theta}_0^{3/2} q_i, \quad \hat{m}_{ijk} := \varepsilon \hat{\rho}_0 \hat{\theta}_0^{3/2} m_{ijk}, \quad \hat{\varDelta} := \varepsilon \hat{\rho}_0 \hat{\theta}_0^{2} \varDelta, \quad \hat{R}_{ij} := \varepsilon \hat{\rho}_0 \hat{\theta}_0^{2} R_{ij}, \end{gathered} \right\} \end{equation}
where 
$\rho$, 
$v_{i}$, 
$T$, 
$\sigma _{ij}$, 
$q_i$, 
$m_{ijk}$, 
$\varDelta$ and 
$R_{ij}$ are the dimensionless deviations of the respective quantities from their equilibrium values. The microscopic velocity 
$\hat {c}_i$ is also scaled to give the dimensionless microscopic velocity 
$\hat {c}_i = \sqrt {\hat {\theta }_0}\, c_i$.
A Taylor series expansion of 
$\hat {f}_{{M}}$, to second order in 
$\varepsilon$, yields
\begin{equation} \hat{f}_{{M}} \approx \hat{f}_{M}^{(0)} \left[1+\varepsilon f_{M}^{(1)}+\varepsilon^{2}f_{M}^{(2)}\right], \end{equation}where
\begin{equation} \hat{f}_{M}^{(0)}=\frac{\hat{\rho}_{0}}{(2{\rm \pi} \hat{\theta}_{0})^{3/2}} \exp \left( -\frac{c^{2}}{2}\right) \end{equation}is the Maxwellian distribution function in the equilibrium, and
$$\begin{gather} f_{M}^{(1)} =\rho -\frac{3T }{2} +\frac{T c^{2}}{2}+c_{j}v_{j}, \end{gather}$$
$$\begin{gather}f_{M}^{(2)} ={-}\frac{3\rho T }{2} + \frac{15T ^{2}}{8}-\frac{v^{2}}{2} +\left(\rho T - \frac{5T ^{2}}{2}+ v^2 +\frac{T ^{2}c^{2}}{4}\right) \frac{c^{2}}{2} +\left(\rho - \frac{5T }{2} + \frac{T c^{2}}{2}\right) c_{j}v_{j}. \end{gather}$$
Similarly, a Taylor series expansion of 
$\varPhi$, to second order in 
$\varepsilon$, yields
\begin{equation} \varPhi \approx \varepsilon \varPhi^{(1)}+\varepsilon^{2}\varPhi^{(2)}, \end{equation}where
$$\begin{align}
\varPhi^{(1)} &= \frac{\sigma_{kl}}{2} c_k
c_l - q_k c_k \left(1-\frac{c^2}{5}\right) +
\frac{m_{klr}}{6} c_k c_l c_r \nonumber\\ &\quad +
\frac{\varDelta}{8} \left(1 - \frac{2 c^2}{3} +
\frac{c^4}{15}\right) -\frac{R_{kl}}{4}
\left(1-\frac{c^2}{7}\right) c_k c_l,
\end{align}$$
\begin{align} \varPhi^{(2)} &={-} \frac{\sigma_{kl}}{2} \left[(\rho + 2 T ) c_k c_l + c_k v_l + c_l v_k \right] \nonumber\\ &\quad + q_k \left[(\rho + 2 T ) c_k - \frac{1}{5} (\rho + 3 T ) c^2 c_k + \left(1-\frac{c^2}{5}\right) v_k - \frac{2}{5} c_k c_l v_l\right] \nonumber\\ &\quad - \frac{m_{kl r}}{6} \left[(\rho + 3 T ) c_k c_l c_r + c_k c_l v_r + c_k c_r v_l+ c_l c_r v_k \right] \nonumber\\ &\quad -\frac{\varDelta}{8} \left[\rho + 2 T - \frac{2}{3} (\rho + 3 T ) c^2 + \frac{1}{15} (\rho + 4 T ) c^4 - \frac{4}{3} \left(1-\frac{c^2}{5}\right) c_k v_k\right] \nonumber\\ &\quad + \frac{R_{kl}}{4} \left[\left(1-\frac{c^2}{7} \right) \left\{(\rho + 3 T ) c_k c_l + c_k v_l + c_l v_k \right\} -\frac{T}{7} c^2 c_k c_l -\frac{2}{7} c_k c_l c_r v_r\right]. \end{align}
Therefore, to second order in 
$\varepsilon$, 
$\hat {f}_{{G26}}$ can be approximated as
\begin{equation} \hat{f}_{{G26}} \approx \hat{f}_{M}^{(0)}[1+\varepsilon \phi^{(1)} +\varepsilon^{2} \phi^{(2)}], \end{equation}where
\begin{equation} \phi^{(1)} = \varPhi^{(1)}+f_{M}^{(1)} \quad\text{and}\quad \phi^{(2)} = \varPhi^{(2)} + \varPhi^{(1)} f_{M}^{(1)} + f_{M}^{(2)}. \end{equation}
Let 
$\hat {y} = \hat {f}_{M}^{(0)}$. Hence, to second order in 
$\varepsilon$,
\begin{equation} \ln{\frac{\hat{f}_{{G26}}}{\hat{f}_{M}^{(0)}}}\approx \ln{\left( 1+\varepsilon \phi^{(1)} +\varepsilon^{2} \phi^{(2)}\right)} \approx \varepsilon \phi^{(1)} +\varepsilon^{2} \left(\phi^{(2)}-\frac{ \left\{ \phi^{(1)}\right\}^{2}}{2}\right). \end{equation}
From (3.2), the dimensionless entropy 
$\eta = \hat {\eta } / (\hat {k}_B \hat {\rho }_0)$, to second order in 
$\varepsilon$, is given by
\begin{align} \eta &={-} \frac{1}{\hat{\rho}_0} \int \hat{f}_{{G26}}\left( \ln{\frac{\hat{f}_{{G26}}}{\hat{f}_{M}^{(0)}}}-1\right) \,\mathrm{d}\hat{\boldsymbol{c}} \nonumber\\ &\approx{-} \frac{1}{\hat{\rho}_0}\int \hat{f}_{M}^{(0)}\left( 1 +\varepsilon \phi^{(1)} +\varepsilon^{2} \phi^{(2)}\right) \left[ \varepsilon \phi ^{(1)}+\varepsilon^{2}\left(\phi^{(2)}-\frac{ \left\{\phi^{(1)}\right\}^{2}}{2}\right) -1\right] (\hat{\theta}_0^{3/2} \,\mathrm{d}\boldsymbol{c}) \nonumber\\ &= \frac{1}{(2{\rm \pi})^{3/2}} \int \exp{\left(-\frac{c^{2}}{2}\right)} \left(1 - \varepsilon^{2}\frac{\left\{\phi^{(1)}\right\}^{2}}{2}\right) \,\mathrm{d}\boldsymbol{c}. \end{align}The integration over the velocity space is performed in a straightforward way, which yields
\begin{equation} \eta = 1 + \varepsilon^{2} \left(-\tfrac{1}{2}\rho ^{2}-\tfrac{1}{2} v^{2} -\tfrac{3}{4} T^{2}-\tfrac{1}{4}\sigma ^{2}-\tfrac{1}{5}q^{2} - \tfrac{1}{12}m^{2} - \tfrac{1}{56}R^{2} - \tfrac{1}{240}\varDelta^{2}\right). \end{equation}
Substituting 
$\varepsilon = 1$, we obtain (3.4) with 
$a_{0}=1$.
From (3.3), the dimensionless entropy flux 
$\varGamma _i = \hat {\varGamma }_i / (\hat {k}_B \hat {\rho }_0 \sqrt {\hat {\theta }_0})$, to second order in 
$\varepsilon$, is given by
\begin{align} \varGamma_i &={-} \frac{1}{\hat{\rho}_0 \sqrt{\hat{\theta}_0}} \int \hat{f}_{{G26}}\left( \ln{\frac{\hat{f}_{{G26}}}{\hat{f}_{M}^{(0)}}}-1\right) \hat{c}_i \,\mathrm{d}\hat{\boldsymbol{c}} \nonumber\\ &\approx{-} \frac{\hat{\theta}_0^{3/2}}{\hat{\rho}_0}\int \hat{f}_{M}^{(0)}\left( 1 +\varepsilon \phi^{(1)} +\varepsilon^{2} \phi^{(2)}\right) \left[ \varepsilon \phi ^{(1)}+\varepsilon^{2}\left(\phi^{(2)}-\frac{ \left\{\phi^{(1)}\right\}^{2}}{2}\right) -1\right] c_i \,\mathrm{d}\boldsymbol{c} \nonumber\\ &= \frac{1}{(2{\rm \pi})^{3/2}} \int \exp{\left(-\frac{c^{2}}{2}\right)} \left(1 - \varepsilon^{2}\frac{\left\{\phi^{(1)}\right\}^{2}}{2}\right) c_i \, \mathrm{d}\boldsymbol{c}. \end{align}Again, the integration over the velocity space is performed in a straightforward way, which yields
$$\begin{gather} \varGamma_i ={-} \varepsilon^{2} \left[(\rho + T ) v_{i} + T q_{i} + v_{j}\sigma_{ij} +\tfrac{2}{5}\sigma_{ij} q_j +\tfrac{1}{2} m_{ijk} \sigma_{jk} \right. \nonumber\\ \left. +\tfrac{1}{5}R_{ij}q_{j} +\tfrac{1}{15}q_{i} \varDelta +\tfrac{1}{14} m_{ijk}R_{jk} \right]. \end{gather}$$
Substituting 
$\varepsilon = 1$ and 
$\rho +T = p$, we obtain
\begin{align} \varGamma_i ={-} \left(p v_{i} + T q_{i} + v_{j}\sigma_{ij} +\tfrac{2}{5}\sigma_{ij} q_j +\tfrac{1}{2} m_{ijk} \sigma_{jk} +\tfrac{1}{5}R_{ij}q_{j} +\tfrac{1}{15}q_{i} \varDelta +\tfrac{1}{14} m_{ijk}R_{jk}\right). \end{align}
Comparing (A19) and (3.7), it is evident that except for the terms containing the higher-order moments, namely 
$\varPhi _{ijkl}$, 
$\varPsi _{ijk}$, 
$\varOmega _{j}$, 
$\varOmega _{i}$, all others terms in (A19) and (3.7) are exactly the same.
Appendix B. The entropy density for the linearised R26 equations
Let us assume that the entropy density for the LR26 equations is given by (3.5). Taking the partial derivative of (3.5) with respect to 
$t$, we obtain
$$\begin{align} \frac{\partial \eta}{\partial t} &= a_1 \rho \frac{\partial \rho}{\partial t} +a_2 v_i \frac{\partial v_i}{\partial t} + a_3 T \frac{\partial T}{\partial t} + a_4 \sigma_{ij} \frac{\partial \sigma_{ij}}{\partial t} + a_5 q_i \frac{\partial q_i}{\partial t} \nonumber\\ &\quad + a_6 m_{ijk} \frac{\partial m_{ijk}}{\partial t} + a_7 R_{ij} \frac{\partial R_{ij}}{\partial t} + a_8 \varDelta \frac{\partial \varDelta}{\partial t}. \end{align}$$
Replacing the time derivatives of field variables in the above equation using the G26 equations (2.3)–(2.7) (with 
$F_i = 0$ and 
$p \approx \rho + T$) and (2.11)–(2.13), the time derivative of 
$\eta$ can be simplified as follows:
\begin{align} \frac{\partial \eta}{\partial t} &={-}a_1 \rho \frac{\partial v_i}{\partial x_i} -a_2 v_i \left(\frac{\partial \sigma_{ij}}{\partial x_j} + \frac{\partial \rho}{\partial x_i} + \frac{\partial T}{\partial x_i}\right) -\frac{2}{3} a_3 T \left(\frac{\partial q_i}{\partial x_i} + \frac{\partial v_i}{\partial x_i}\right) \nonumber\\ &\quad- a_4 \sigma_{ij} \left(\frac{\partial m_{ijk}}{\partial x_k} + \frac{4}{5} \frac{\partial q_{\langle i}}{\partial x_{j\rangle}} +\frac{1}{{Kn}} \sigma_{ij} + 2 \frac{\partial v_{\langle i}}{\partial x_{j\rangle}}\right) \nonumber\\ &\quad- a_5 q_i \left(\frac{1}{2} \frac{\partial R_{ij}}{\partial x_j} +\frac{1}{6} \frac{\partial \varDelta}{\partial x_i} +\frac{\partial \sigma_{ij}}{\partial x_j} +\frac{{Pr}}{{Kn}} q_{i} + \frac{5}{2} \frac{\partial T}{\partial x_i} \right) \nonumber\\ &\quad- a_6 m_{ijk} \left(\frac{\partial \varPhi_{ijkl}}{\partial x_l} +\frac{3}{7} \frac{\partial R_{\langle ij}}{\partial x_{k\rangle}} + \frac{{Pr}_{m}}{{Kn}} m_{ijk} + 3 \frac{\partial \sigma_{\langle ij}}{\partial x_{k\rangle}} \right) \nonumber\\ &\quad- a_7 R_{ij} \left(\frac{\partial \varPsi_{ijk}}{\partial x_k} +2 \frac{\partial m_{ijk}}{\partial x_k} +\frac{2}{5} \frac{\partial \varOmega_{\langle i}}{\partial x_{j\rangle}} + \frac{{Pr}_{R}}{{Kn}} R_{ij} + \frac{28}{5} \frac{\partial q_{\langle i}}{\partial x_{j\rangle}}\right) \nonumber\\ &\quad- a_8 \varDelta \left(\frac{\partial \varOmega_i}{\partial x_i} + \frac{{Pr}_{\varDelta }}{{Kn}} \varDelta + 8 \frac{\partial q_i}{\partial x_i}\right) \end{align}or
\begin{align} \frac{\partial \eta}{\partial t} &={-} \left(a_1 \rho \frac{\partial v_i}{\partial x_i} + a_2 v_i \frac{\partial \rho}{\partial x_i} \right) - \left( a_2 v_i\frac{\partial T}{\partial x_i} + \frac{2}{3} a_3 T \frac{\partial v_i}{\partial x_i}\right) - \left(a_2 v_i \frac{\partial \sigma_{ij}}{\partial x_j} + 2 a_4 \sigma_{ij} \frac{\partial v_i}{\partial x_j}\right) \nonumber\\ &\quad- \left(\frac{2}{3} a_3 T \frac{\partial q_i}{\partial x_i} + \frac{5}{2} a_5 q_i \frac{\partial T}{\partial x_i} \right) - \left(a_4 \sigma_{ij} \frac{\partial m_{ijk}}{\partial x_k} + 3 a_6 m_{ijk} \frac{\partial \sigma_{ij}}{\partial x_k}\right) \nonumber\\ &\quad- \left(\frac{4}{5} a_4 \sigma_{ij} \frac{\partial q_i}{\partial x_j} + a_5 q_i \frac{\partial \sigma_{ij}}{\partial x_j}\right) - \left(\frac{1}{2} a_5 q_i \frac{\partial R_{ij}}{\partial x_j} + \frac{28}{5} a_7 R_{ij} \frac{\partial q_i}{\partial x_j}\right) \nonumber\\ &\quad+\left(\frac{1}{6} a_5 q_i \frac{\partial \varDelta}{\partial x_i} + 8 a_8 \varDelta \frac{\partial q_i}{\partial x_i}\right) - \left(\frac{3}{7} a_6 m_{ijk} \frac{\partial R_{ij}}{\partial x_k} + 2 a_7 R_{ij} \frac{\partial m_{ijk}}{\partial x_k} \right) \nonumber\\ &\quad- a_6 m_{ijk} \frac{\partial \varPhi_{ijkl}}{\partial x_l} - a_7 R_{ij} \frac{\partial \varPsi_{ijk}}{\partial x_k} - \frac{2}{5} a_7 R_{ij} \frac{\partial \varOmega_{\langle i}}{\partial x_{j\rangle}} - a_8 \varDelta \frac{\partial \varOmega_i}{\partial x_i} \nonumber\\ &\quad- \frac{a_4}{{Kn}} \sigma^2 - a_5 \frac{{Pr}}{{Kn}} q^2 - a_6 \frac{{Pr}_{m}}{{Kn}} m^2 - a_7 \frac{{Pr}_{R}}{{Kn}} R^2 - a_8 \frac{{Pr}_{\varDelta}}{{Kn}} \varDelta^2. \end{align}The goal is to write the above equation in the form of (3.1), which can be accomplished by taking
\begin{equation} \left. \begin{gathered} a_1 = a_2, \quad a_2 = \tfrac{2}{3} a_3, \quad a_2 = 2 a_4, \quad \tfrac{2}{3} a_3 = \tfrac{5}{2} a_5, \quad a_4 = 3 a_6, \\ \tfrac{4}{5} a_4 = a_5, \quad \tfrac{1}{2} a_5 = \tfrac{28}{5} a_7, \quad \tfrac{1}{6} a_5 = 8 a_8, \quad \tfrac{3}{7} a_6 = 2 a_7. \end{gathered} \right\} \end{equation}
These are nine equations in eight variables. Nevertheless, they are consistent, and only seven of them are linearly independent. Thus, on taking 
$a_1=-1$, the other coefficients are
\begin{align} a_2 ={-}1, \quad a_3 ={-}\tfrac{3}{2}, \quad a_4 ={-}\tfrac{1}{2}, \quad a_5 ={-}\tfrac{2}{5}, \quad a_6 ={-}\tfrac{1}{6}, \quad a_7 ={-}\tfrac{1}{28}, \quad a_8 ={-}\tfrac{1}{120}. \end{align}
Note that the constant 
$a_1$ has been taken as negative in order to get the non-negative entropy production rate. Consequently, the entropy (3.5) becomes (3.4).
Appendix C. Entropy flux at the interface
The entropy flux at the interface is given by (4.1). Substituting 
$\varGamma _{i}^{{(gas)}} = \varGamma _{i}$ from (3.7) and 
$\varGamma _{i}^{{(liquid)}}$ from (4.2) into (4.1), the entropy flux at the interface reads as
\begin{align} \varSigma^{I}&={-} \underline{\left(p v_{i} - p^{\ell} v_{i}^{\ell} + T q_{i} - T^{\ell} q_{i}^{\ell} + v_{j}\sigma_{ij} - v_{j}^{\ell}\sigma_{ij}^{\ell} \right) n_i} \nonumber\\ &\quad- \left(\tfrac{2}{5}\sigma_{ij} q_j +\tfrac{1}{2} m_{ijk} \sigma_{jk} +\tfrac{1}{5}R_{ij}q_{j} +\tfrac{1}{15}q_{i} \varDelta +\tfrac{1}{14} m_{ijk}R_{jk} \right. \nonumber\\ &\quad\left.+\tfrac{1}{6}\varPhi_{ijkl}m_{jkl} +\tfrac{1}{28}\varPsi_{ijk}R_{jk} +\tfrac{1}{70}R_{ij}\varOmega_{j} +\tfrac{1}{120} \varOmega_{i} \varDelta\right)n_{i}. \end{align}
Let us refer the underlined term in (C1) by 
$\underline {\varSigma }^{I}$, i.e.
$$\begin{align}
\underline{\varSigma}^{I} &= (p v_{i} - p^{\ell}
v_{i}^{\ell} + T q_{i} - T^{\ell} q_{i}^{\ell} +
v_{j}\sigma_{ij} - v_{j}^{\ell}\sigma_{ij}^{\ell}) n_i
\nonumber\\ &=p v_{i} n_i - p^{\ell} v_{i}^{\ell} n_i + T
q_{i} n_i - T^{\ell} q_{i}^{\ell} n_i + v_i \sigma_{ij} n_j
- v_i^{\ell}\sigma_{ij}^{\ell} n_j .
\end{align}$$
Equation (C2) is simplified using the conservation laws at the interface. For a simple interface (assumed here), the conservation laws for the mass, momentum and energy apprise the continuity of the normal fluxes (Struchtrup et al. Reference Struchtrup, Beckmann, Rana and Frezzotti2017). With the gas (or vapour) properties on the left-hand sides and the liquid properties on the right-hand sides, the conservation laws at the interface in linear-dimensionless form read as (Müller Reference Müller1985; Struchtrup et al. Reference Struchtrup, Beckmann, Rana and Frezzotti2017)
$$\begin{gather} {(}v_i -v_i^{I}) n_i =\rho^{\ell} (v_i^{\ell}-v_i^{I}) n_i, \end{gather}$$
$$\begin{gather} {(}p \delta_{ij}+\sigma_{ij}) n_j = ( p^{\ell} \delta_{ij}+\sigma_{ij}^{\ell}) n_j , \end{gather}$$
$$\begin{gather} {[}h_0 (v_i - v_i^{I}) + q_i] n_i = [\rho^{\ell} h_0^{\ell} ( v_i^{\ell} - v_i^{I}) + q_i^{\ell}] n_i, \end{gather}$$
where 
$v_i^{I}$ is the velocity of the interface, and 
$h_0$ and 
$h_0^{\ell }$ are the dimensionless enthalpies of the vapour and liquid, respectively, at temperature 
$\hat {T}_0$. The variables in (C3)–(C5) have been non-dimensionalised using 
$\hat {\rho }_0$, 
$\hat {\theta }_0$ and their appropriate combinations. Substituting 
$\sigma _{ij}^{\ell } n_j$ and 
$q_i^{\ell } n_i$ from (C4) and (C5) in (C2), the underlined term in (C1) becomes
$$\begin{align} \underline{\varSigma}^{I} &=p (v_i - v_i^{\ell}) n_i + (T- T^{\ell}) q_{i} n_i - T^{\ell} (h_0 - h_0^{\ell}) (v_i -v_i^{I}) n_i + (v_j - v_j^{\ell}) \sigma_{ij} n_i \nonumber\\ &=p \mathscr{V}_i n_i + \mathcal{T} q_{i} n_i - T^{\ell} (h_0 - h_0^{\ell}) \mathcal{V}_i n_i + \mathscr{V}_j \sigma_{ij} n_i, \end{align}$$
where 
$\mathscr {V}_i = v_i - v_i^{\ell }$, 
$\mathcal {T} = T- T^{\ell }$, 
$\mathcal {V}_i = v_i -v_i^{I}$. On using (C3), it turns out that
\begin{equation} \mathscr{V}_i \, n_i =(v_i - v_i^{\ell}) n_i = [\mathcal{V}_i - (v_i^{\ell} - v_i^{I})] n_i = \left(1 - \frac{1}{\rho^{\ell}}\right) \mathcal{V}_i n_i \approx \mathcal{V}_i \, n_i \end{equation}
because the density of a gas/vapour is usually much smaller than the density of a liquid/solid, and, hence, 
$\rho ^{\ell } \gg 1$ or 
$1/\rho ^{\ell } \ll 1$. The enthalpy term in (C6) is eliminated using the definition of the (dimensionless) saturation pressure 
$p_{{sat}}$ at temperature 
$T^{\ell }$ (Rana et al. Reference Rana, Lockerby and Sprittles2019),
\begin{equation} p_{{sat}} \equiv p_{{sat}} (T^{\ell}) = T^{\ell} (h_0 - h_0^{\ell}). \end{equation}Using (C7) and (C8), (C6) reduces to
\begin{equation} \underline{\varSigma}^{I} =p \mathcal{V}_i \, n_i + \mathcal{T} q_{i} n_i - p_{{sat}} \mathcal{V}_i \, n_i + \mathscr{V}_j \sigma_{ij} n_i = \left(\mathcal{P} \mathcal{V}_i + \mathcal{T} q_{i} + \mathscr{V}_j \sigma_{ij}\right) n_i, \end{equation}
where 
$\mathcal {P} = p - p_{{sat}}$. Replacing the underlined terms in (C1) with the value of 
$\underline {\varSigma }^{I}$ obtained in (C9), the entropy flux at the interface is given by (4.3).
Appendix D. Decomposition of symmetric-tracefree tensors
In the following, we decompose vectors and symmetric-tracefree tensors of rank two and three into their components in the direction of a given normal 
$\boldsymbol {n}$ and in the directions perpendicular to it (mutually perpendicular tangential directions, let us say 
$\boldsymbol {t}_1$ and 
$\boldsymbol {t}_2$). Obviously, 
$n_i {t_1}_i = n_i {t_2}_i = {t_1}_i {t_2}_i = 0$. When it is not required to distinguish between 
$\boldsymbol {t}_1$ and 
$\boldsymbol {t}_2$, both are referred to with a common notation 
$\boldsymbol {t}$ for simplicity, and 
$\boldsymbol {t}$ could either be 
$\boldsymbol {t}_1$ and 
$\boldsymbol {t}_2$.
(i) A vector
$a_i$ is decomposed as
(D 1)where\begin{equation} a_{i} = a_{n} n_{i} + \bar{a}_{i}, \end{equation}$a_{n} n_{i}$ with $a_n := a_k n_k$ is the normal component of $a_i$ while $\bar {a}_{i} := a_{i} - a_{n} n_{i}$ is the tangential component of $a_i$. By definition, $\bar {a}_{i}$ is such that $\bar {a}_{i} n_i = 0$. Furthermore, multiplication of (D1) with $t_i \in \{{t_1}_i,{t_2}_i\}$ gives
(D 2)where\begin{equation} a_t=\bar{a}_t, \end{equation}$a_t:=a_i t_i$, $\bar {a}_t:=\bar {a}_{i} t_i$.(ii) A symmetric-tracefree rank-2 tensor
$A_{ij}$ is decomposed as (Rana & Struchtrup Reference Rana and Struchtrup2016; Rana et al. Reference Rana, Gupta and Struchtrup2018a)
(D 3)where Einstein summation never applies to the indices ‘\begin{equation} A_{ij} = \tfrac{3}{2} A_{nn} \, n_{\langle i}n_{j \rangle} + \bar{A}_{ni} n_{j} + \bar{A}_{nj} n_{i} + \tilde{A}_{ij}, \end{equation}$n$’, and
(D 4a–c)It is easy to verify that\begin{equation} \left. \begin{aligned} A_{nn} & := A_{ij} n_i n_j, \\ \bar{A}_{ni} & := A_{ij} n_{j} - A_{nn} n_{i}, \\ \tilde{A}_{ij} & := A_{ij} - \tfrac{3}{2} A_{nn} \, n_{\langle i}n_{j \rangle} - \bar{A}_{ni} n_{j} - \bar{A}_{nj} n_{i}. \end{aligned} \right\} \end{equation}$\bar {A}_{ni} n_{i} = 0$ and $\tilde {A}_{kk} = \tilde {A}_{ij} n_{i} = \tilde {A}_{ij} n_{j} = 0$. Furthermore, multiplication of (D4b) with $t_i \in \{{t_1}_i,{t_2}_i\}$ gives
(D 5)where\begin{equation} \bar{A}_{nt} = A_{nt}, \end{equation}$\bar {A}_{nt} := \bar {A}_{ni} t_i$, $A_{nt}=A_{tn}:= A_{ij} t_i n_{j}$; multiplication of (D4c) with $t_i t_j$ gives
(D 6)where\begin{equation} \tilde{A}_{tt} = A_{tt} + \tfrac{1}{2} A_{nn}, \end{equation}$\tilde {A}_{tt} := \tilde {A}_{ij} t_i t_j$ and $A_{tt} := A_{ij} t_i t_j$; multiplication of (D4c) with ${t_1}_i {t_2}_j$ gives
(D 7)where\begin{equation} \tilde{A}_{{t_1}{t_2}} = A_{{t_1}{t_2}}, \end{equation}$\tilde {A}_{{t_1}{t_2}}: = \tilde {A}_{ij} {t_1}_i {t_2}_j$ and $A_{{t_1}{t_2}} := A_{ij} {t_1}_i {t_2}_j$.(iii) A symmetric-tracefree rank-3 tensor
$A_{ijk}$ is decomposed as
(D 8)where\begin{equation} A_{ijk} = \tfrac{5}{2} A_{nnn} n_{\langle i}n_{j}n_{k\rangle} +\tfrac{15}{4}\bar{A}_{nn\langle i}n_{j}n_{k\rangle} + 3\tilde{A}_{n\langle ij}n_{k\rangle} + \check{A}_{ijk}. \end{equation}(D 9a–d)It is again easy to verify that\begin{equation} \left. \begin{aligned} A_{nnn} & := A_{ijk} n_{i} n_{j}n_{k}, \\ \bar{A}_{nni} & := A_{ijk}n_{j}n_{k}-A_{nnn}n_{i}, \\ \tilde{A}_{nij} & := A_{ijk} n_{k} - \tfrac{3}{2} A_{nnn} n_{\langle i}n_{j\rangle} - \bar{A}_{nni} n_{j} - \bar{A}_{nnj} n_{i}, \\ \check{A}_{ijk} & := A_{ijk} - \tfrac{5}{2} A_{nnn} n_{\langle i}n_{j}n_{k\rangle} - \tfrac{15}{4}\bar{A}_{nn\langle i}n_{j}n_{k\rangle} - 3\tilde{A}_{n\langle ij}n_{k\rangle}. \end{aligned} \right\} \end{equation}$\bar {A}_{nni}n_{i}=\tilde {A}_{nij}n_{j}=\tilde {A}_{nii}=\check {A}_{ijk}n_{k}=\check {A}_{ikk}=0$. Furthermore, multiplication of (D9b) with $t_i \in \{{t_1}_i,{t_2}_i\}$ gives
(D 10)where\begin{equation} \bar{A}_{nnt} = A_{nnt}, \end{equation}$\bar {A}_{nnt} := \bar {A}_{nni} t_i$, $A_{nnt}=A_{tnn}:= A_{ijk}t_i n_{j}n_{k}$; multiplication of (D9c) with $t_i t_j$ gives
(D 11)where\begin{equation} \tilde{A}_{ntt} = A_{ntt} + \tfrac{1}{2} A_{nnn}, \end{equation}$\tilde {A}_{ntt}: = \tilde {A}_{nij} t_i t_j$ and $A_{ntt} = A_{ttn} := A_{ijk} t_i t_j n_{k}$; multiplication of (D9c) with ${t_1}_i {t_2}_j$ gives
(D 12)where\begin{equation} \tilde{A}_{n{t_1}{t_2}} = A_{n{t_1}{t_2}}, \end{equation}$\tilde {A}_{n {t_1}{t_2}}: = \tilde {A}_{nij} {t_1}_i {t_2}_j$ and $A_{n{t_1}{t_2}} = A_{{t_1}{t_2}n} := A_{ijk} {t_1}_i {t_2}_j n_{k}$; multiplication of (D9d) with $t_i t_j t_k$ gives
(D 13)where\begin{equation} \check{A}_{ttt} = A_{ttt} + \tfrac{3}{4} A_{n nt}, \end{equation}$\check {A}_{ttt} := \check {A}_{ijk} t_i t_j t_k$ and $A_{ttt} := A_{ijk} t_i t_j t_k$; multiplication of (D9d) with ${t_1}_i {t_1}_j {t_2}_k$ gives
(D 14)where\begin{equation} \check{A}_{{t_1} {t_1} {t_2}} = A_{{t_1} {t_1} {t_2}} + \tfrac{1}{4} A_{n n{t_2}}, \end{equation}$\check {A}_{{t_1} {t_1} {t_2}} := \check {A}_{ijk} {t_1}_i {t_1}_j {t_2}_k$ and $A_{{t_1} {t_1} {t_2}} := A_{ijk} {t_1}_i {t_1}_j {t_2}_k$.
Appendix E. Analytic solution for the unknowns in ansatz (6.13)–(6.15)
The analytic solution for the unknowns appearing in (6.13)–(6.15) are as follows:
\begin{align} a_{1}(r) &=
\frac{C_1}{2r}+\frac{C_2}{3r^3}-K_1 \varTheta_1(r)
\left(\frac{0.975701 {Kn}^3}{r^3}+\frac{0.497885
{Kn}^2}{r^2}\right) \nonumber\\ &\quad - K_2 \varTheta_2(r)
\left(\frac{1.4893 {Kn}^3}{r^3}+\frac{1.88528
{Kn}^2}{r^2}\right),\\ a_{2}(r) &=
\frac{C_1}{4r}-\frac{C_2}{6r^3}+K_1 \varTheta_1(r)
\left(\frac{0.487851 {Kn}^3}{r^3}+\frac{0.248943
{Kn}^2}{r^2}+\frac{0.127032 {Kn}}{r}\right) \nonumber\\
&\quad +K_2 \varTheta_2(r) \left(\frac{0.74465
{Kn}^3}{r^3}+\frac{0.942641 {Kn}^2}{r^2}+\frac{1.19327
{Kn}}{r}\right), \\ \alpha_{1}(r) &= \frac{{Kn}^2
}{r^3}\left(\frac{C_3}{6{Kn}}-\frac{3 C_1}{2}\right) +K_1
\varTheta_1(r) \left(\frac{1.50825
{Kn}^3}{r^3}+\frac{0.769634 {Kn}^2}{r^2}\right) \nonumber\\
&\quad -K_2 \varTheta_2(r) \left(\frac{0.80299
{Kn}^3}{r^3}+\frac{1.01649 {Kn}^2}{r^2}\right), \\
\alpha_{2}(r) &= \frac{{Kn}^2
}{r^3}\left(\frac{3C_1}{4}-\frac{C_3}{12 {Kn}}\right)
\nonumber\\ &\quad-K_1 \varTheta_1(r) \left(\frac{0.754123
{Kn}^3}{r^3}+\frac{0.384817 {Kn}^2}{r^2}+\frac{0.196366
{Kn}}{r}\right)\nonumber\\
&\quad+ K_2 \varTheta_2(r) \left(\frac{0.401495
{Kn}^3}{r^3}+\frac{0.508246 {Kn}^2}{r^2}+\frac{0.643381
{Kn}}{r}\right),
\end{align}
\begin{align} c(r) &= \frac{C_3}{45
r^2}-K_3 \varTheta_3(r) \left(\frac{0.0172036
{Kn}^2}{r^2}+\frac{0.0200113 {Kn}}{r}\right) \nonumber\\
&\quad -K_4 \varTheta_4(r) \left(\frac{0.0629047
{Kn}^2}{r^2}+\frac{0.0426083 {Kn}}{r}\right) \nonumber\\
&\quad -K_5 \varTheta_5(r) \left(\frac{0.443094
{Kn}^2}{r^2}+\frac{0.200538 {Kn}}{r}\right),\\
d(r) &=\frac{C_1{Kn}}{2r^2}-K_3\varTheta_3(r)
\left(\frac{0.931108 {Kn}^2}{r^2}+\frac{1.08307
{Kn}}{r}\right) \nonumber\\ &\quad-K_4 \varTheta_4(r)
\left(\frac{0.254173 {Kn}^2}{r^2}+\frac{0.172163
{Kn}}{r}\right) \nonumber\\ &\quad-K_5 \varTheta_5(r)
\left(\frac{0.0501203 {Kn}^2}{r^2}+\frac{0.0226838
{Kn}}{r}\right), \\ \delta(r) &=K_3 \varTheta_3(r)
\left(\frac{14.0055 {Kn}^2}{r^2}+\frac{16.2913
{Kn}}{r}\right)\nonumber\\
&\quad -K_4 \varTheta_4(r) \left(\frac{1.71119
{Kn}^2}{r^2}+\frac{1.15907 {Kn}}{r}\right) \nonumber\\[2pt]
&\quad+K_5 \varTheta_5(r)\left(\frac{3.75265 {Kn}^2}{r^2}+\frac{1.6984 {Kn}}{r}\right),
\end{align}
\begin{align} b_1(r) &=\frac{C_1
{Kn}}{r^2}+\frac{{Kn}^4}{r^4}\left(-\frac{10
C_1}{{Kn}}+\frac{2 C_2}{{Kn}^3}+\frac{2
C_3}{5{Kn}^2}\right) \nonumber\\[2pt] &\quad+K_3 \varTheta_3(r)
\left(\frac{6.19341 {Kn}^4}{r^4}+\frac{7.20421
{Kn}^3}{r^3}+\frac{3.72443 {Kn}^2}{r^2}+\frac{1.08307
{Kn}}{r}\right) \nonumber\\[2pt] &\quad+K_4 \varTheta_4(r)
\left(\frac{4.98597 {Kn}^4}{r^4}+\frac{3.37723
{Kn}^3}{r^3}+\frac{1.01669 {Kn}^2}{r^2}+\frac{0.172163
{Kn}}{r}\right) \nonumber\\[2pt] &\quad+K_5 \varTheta_5(r)
\left(\frac{2.20218 {Kn}^4}{r^4}+\frac{0.996676
{Kn}^3}{r^3}+\frac{0.200481 {Kn}^2}{r^2}+\frac{0.0226838
{Kn}}{r}\right), \\[2pt] b_2(r)
&=\frac{{Kn}^4}{r^4}\left(-\frac{5
C_1}{{Kn}}+\frac{C_2}{{Kn}^3}+\frac{\text{C3}}{5{Kn}^2}\right)
\nonumber\\[2pt] &\quad+K_3 \varTheta_3(r)\left(\frac{3.09671
{Kn}^4}{r^4}+\frac{3.6021 {Kn}^3}{r^3}+\frac{1.39666
{Kn}^2}{r^2}\right) \nonumber\\[2pt] &\quad+K_4 \varTheta_4(r)
\left(\frac{2.49299 {Kn}^4}{r^4}+\frac{1.68862
{Kn}^3}{r^3}+\frac{0.381259 {Kn}^2}{r^2}\right) \nonumber\\[2pt]
&\quad+K_5 \varTheta_5(r) \left(\frac{1.10109
{Kn}^4}{r^4}+\frac{0.498338 {Kn}^3}{r^3}+\frac{0.0751805
{Kn}^2}{r^2}\right), \\[2pt] \beta_1(r)
&=\frac{{Kn}^4}{r^4}\left(\frac{12C_3}{5{Kn}^2}-\frac{228
C_1}{7{Kn}}\right) \nonumber\\[2pt] &\quad+K_1 \varTheta_1(r)
\left(\frac{46.338 {Kn}^4}{r^4}+\frac{23.6455
{Kn}^3}{r^3}+\frac{4.02199 {Kn}^2}{r^2}\right) \nonumber\\[2pt]
&\quad-K_2 \varTheta_2(r) \left(\frac{4.00877
{Kn}^4}{r^4}+\frac{5.07465 {Kn}^3}{r^3}+\frac{2.14131
{Kn}^2}{r^2}\right) \nonumber\\ &\quad-K_3
\varTheta_3(r)\left(\frac{42.868 {Kn}^4}{r^4}+\frac{49.8643
{Kn}^3}{r^3}+\frac{25.7788 {Kn}^2}{r^2}+\frac{7.49652
{Kn}}{r}\right) \nonumber\\ &\quad+K_4\varTheta_4(r)
\left(\frac{7.38704 {Kn}^4}{r^4}+\frac{5.00359
{Kn}^3}{r^3}+\frac{1.5063 {Kn}^2}{r^2}+\frac{0.255071
{Kn}}{r}\right)\nonumber\\ &\quad+K_5\varTheta_5(r)
\left(\frac{37.9774 {Kn}^4}{r^4}+\frac{17.1881
{Kn}^3}{r^3}+\frac{3.45738 {Kn}^2}{r^2}+\frac{0.391191
{Kn}}{r}\right), \\ \beta_2(r)
&=\frac{{Kn}^4}{r^4}\left(\frac{6C_3}{5{Kn}^2}-\frac{114
C_1}{7{Kn}}\right) \nonumber\\ &\quad+K_1\varTheta_1(r)
\left(\frac{23.169 {Kn}^4}{r^4}+\frac{11.8228
{Kn}^3}{r^3}+\frac{3.01649 {Kn}^2}{r^2}+\frac{0.51309
{Kn}}{r}\right) \nonumber\\
&\quad-K_2\varTheta_2(r)\left(\frac{2.00439
{Kn}^4}{r^4}+\frac{2.53732 {Kn}^3}{r^3}+\frac{1.60598
{Kn}^2}{r^2}+\frac{0.677662 {Kn}}{r}\right) \nonumber\\
&\quad-K_3 \varTheta_3(r)\left(\frac{21.434
{Kn}^4}{r^4}+\frac{24.9321 {Kn}^3}{r^3}+\frac{9.66706
{Kn}^2}{r^2}\right) \nonumber\\
&\quad+K_4\varTheta_4(r)\left(\frac{3.69352
{Kn}^4}{r^4}+\frac{2.50179 {Kn}^3}{r^3}+\frac{0.564861
{Kn}^2}{r^2}\right) \nonumber\\
&\quad+K_5\varTheta_5(r) \left(\frac{18.9887 {Kn}^4}{r^4}+\frac{8.59403 {Kn}^3}{r^3}+\frac{1.29652 {Kn}^2}{r^2}\right),
\end{align}
In (E1)–(E13), the coefficients
\begin{equation} \left. \begin{gathered} \varTheta_1(r) = \exp\left({-\frac{0.510285 (r-1)}{{Kn}}}\right), \quad \varTheta_2(r)=\exp\left({-\frac{1.26588 (r-1)}{{Kn}}}\right), \\ \quad \varTheta_3(r) = \exp\left({-\frac{1.16321 (r-1)}{{Kn}}}\right), \quad \varTheta_4(r) = \exp\left({-\frac{0.677347 (r-1)}{{Kn}}}\right)\\ \quad\text{and}\quad \varTheta_5(r) = \exp\left({-\frac{0.452587 (r-1)}{{Kn}}}\right) \end{gathered} \right\} \end{equation}
describe the Knudsen layer functions that vanish as 
$r \to \infty$.
