Quasi-linear irreversible thermodynamics of a low-temperature-differential kinematic Stirling heat engine

Low-temperature-differential (LTD) Stirling heat engines are able to operate with a small temperature difference between low-temperature heat reservoirs that exist in our daily lives, and thus they are considered to be an important sustainable energy technology. The author recently proposed a nonlinear dynamics model of an LTD kinematic Stirling heat engine to study the rotational mechanism of the engine [Y. Izumida, EPL \textbf{121}, 50004 (2018)]. It was shown that the rotational motion of the engine is described by a stable limit cycle of the dynamical equations, and disappears via a homoclinic bifurcation, with the temperature difference being the bifurcation parameter. This paper presents our study of the nonequilibrium thermodynamics analysis of this engine model. Therefore, we introduce a load torque into the model so that the engine does work against the load torque under a given temperature difference. We demonstrate that the engine begins to rotate by exhibiting a homoclinic bifurcation at the bifurcation point with respect to the load torque. It is found that there is a quasi-linear response regime sufficiently far from the bifurcation points where the averaged fluxes show a linear dependency on the external forces. Significantly, it is found that the response coefficients of the quasi-linear relations are symmetric, which is similar to Onsager symmetry in linear irreversible thermodynamics. Based on these relations, we formulate the maximum efficiency of the engine. We also elucidate that the symmetry of the quasi-linear response coefficients emerges by reflecting the (anti-)reciprocity of the Onsager kinetic coefficients in the relaxation dynamics of the engine in the vicinity of an equilibrium state. The present study is expected to pave the way for developing nonequilibrium thermodynamics of autonomous heat engines described as a nonlinear dynamical system.


I. INTRODUCTION
The development of heat engines that operate with small temperature differences and at low friction is an important task in heat engine technology. This task has been undertaken by low-temperature-differential (LTD) Stirling heat engines [1][2][3]. These heat engines were invented by Kolin in the 1980s and subsequently developed primarily by Kolin and Senft [3]. An LTD Stirling heat engine can operate with a small temperature difference between low-temperature heat reservoirs that are available in everyday life, e.g., between the warmth of our hand and the coldness of air temperature. Thus, it is considered to be an important sustainable energy technology.
Appropriate mathematical modeling plays an important role in describing and understanding the dynamics of LTD Stirling engines [4,5]. The author recently proposed a nonlinear dynamics model of an LTD kinematic Stirling engine to elucidate the rotational mechanism of the engine [6]. In this model, the engine was described as a driven nonlinear pendulum powered by the temperature difference, which obeys simple dynamical equations with only a few dynamical degrees of freedom. The rotational motion of the engine was described as a stable limit cycle of the dynamical equations sustained by the temperature difference. Moreover, it was shown that the limit cycle disappears via a homoclinic bifurcation [7], with the temperature difference being the bifurcation parameter. The model was recently used to explain the experimental results on an LTD kinematic Stirling engine [8]. It was demonstrated that the core dynamics of the engine are * izumida@k.u-tokyo.ac.jp captured by the simple dynamical equations with some modifications that are associated with a few fitting parameters. Although the study in [6] elucidated the rotational mechanism of the engine based on nonlinear dynamics, the thermodynamic performance of the LTD Stirling engine, such as its thermodynamic efficiency, has not yet been formulated. Apart from the present performance of the LTD Stirling engine, it is of interest to formulate its maximum thermodynamic efficiency based on the minimal model. Consequently, we need to develop a nonequilibrium thermodynamic theory of the LTD Stirling engine described as a nonlinear dynamical system.
In this paper, we develop the nonequilibrium thermodynamics of the LTD kinematic Stirling engine model that was previously introduced [6]. In particular, our goal is to find relevant thermodynamic relations that describe the rotational state of this thermodynamic nonlinear pendulum model, which may be compared to Onsager relations [9,10] used in describing linear irreversible heat engines [11][12][13]. We formulate the thermodynamic efficiency of the LTD kinematic Stirling engine model based on these relations.
The remainder of this paper is organized as follows. In Sec. II, we introduce the LTD kinematic Stirling engine model [6]. In Sec. III and Sec. IV, we investigate stationary and rotational states of the engine, respectively, based on the dynamical equations. The formal analytical expressions of the thermodynamic fluxes (angular velocity and heat flux) are derived for the rotational state. In particular, the quasi-linear response regime is identified for the rotational state where the thermodynamic fluxes and forces show linear dependency (quasi-linear relations), though this regime is not connected to an equilibrium state. In Sec. V, we formulate the thermodynamic efficiency of the engine using the quasi-linear relations The displacer that advances π 2 in phase serves to transfer the gas into one side of the cylinder and makes the gas in contact with the bottom and top heat reservoirs. The motive force of the rotation is the temperature difference between the bottom and top heat reservoirs with temperature T b and T t , respectively. (b) Schematic illustration of the pressure-volume diagram of an ideal Stirling cycle (solid outer cycle) consisting of the four thermodynamic processes (see the main text) and the kinematic Stirling engine (dotted inner cycle).
for which the coefficients turn out to be symmetric. In Sec. VI, we elucidate the origin of the symmetric coefficients in terms of (anti-)reciprocity of the Onsager kinetic coefficients inherited in the relaxation dynamics of the engine. We summarize this study in Sec. VII.

A. Setup
We use the same model as in our previous study [6], but with a slight extension to add a load torque, which enables the thermodynamic efficiency to be studied. Because the model was previously explained in detail [6], we introduce it here in a simplified but self-contained manner.
The LTD kinematic Stirling engine, regarded as a γ-type Stirling engine [3], utilizes two connected cylinders (one large and one small) with two movable pistons of different types in these cylinders ( Fig. 1 (a)). The working substance of the engine is a gas that is confined to the cylinders. Heat reservoirs at temperatures T b and T t , such as a warm palm and the cold air surrounding it, are attached to the bottom and top surfaces of the large cylinder, respectively, where we define temperature difference ∆T ≡ T b − T t and averaged temperature T eq ≡ T b +T t 2 for later use. The piston that reciprocates in the large cylinder is a displacer. The motion of the displacer serves to transfer the gas into one side of the cylinders through a gap between the displacer piston and wall of the large cylinder, such that the gas comes into contact with the top and bottom heat reservoirs alternately. In contrast, the small cylinder is fitted with a power piston at the top, and its reciprocating motion constitutes a motive part of the engine. Each piston is connected to a crank with a radius r through a connecting rod, and the reciprocating motion of the power piston is converted into rotational motion via the crank (piston-crank mechanism). The phase angle of the crank connected to the power piston is θ (mod 2π), whereas that of the crank connected to the displacer is fixed as θ + π 2 so that it advances in π 2 . The phase angle θ increases as it rotates clockwise and θ = 0 at the lowest height of the power piston. The cranks are attached to a flywheel with a large moment of inertia I to smoothen the rotation; the engine can continue to maintain rotation by overcoming θ = 0, known as top dead center (TDC), and θ = π, known as bottom dead center (BDC), at which the reciprocating motion of the piston is not transmitted to the crank.
The phase angle θ is one of the dynamical variables that expresses the mechanical degree of freedom of this engine model. The other dynamical variable, as a thermodynamic degree of freedom, is the temperature T of the gas. We assume an ideal gas with f internal degrees of freedom as the working substance, for which the equation of state pV = nRT holds. Here, p and V are the pressure and volume of the gas, respectively, and n and R are the amount of substance and gas constant, respectively. The volume V is calculated as the sum of the volume of the large cylinder excluding the volume of the displacer, V d , and that of the small cylinder, V p (θ): where σ d and σ p are the surface areas of the large and small cylinders, respectively, and is the height of the power piston measured from the lowest position corresponding to θ = 0. An ideal Stirling engine cycle repeats an (I) isochoric heating process, (II) isothermal expansion process, (III) isochoric cooling process, and (IV) isothermal compression process [1,2], whose pressure-volume diagram is shown in Fig. 1 (b). Conversely, the pressure-volume diagram of an LTD Stirling engine is presented as a circular shape as shown in Fig. 1 (b), which is observed in the experiments on LTD kinematic Stirling engines [8,14]. While the above thermodynamic processes of the ideal cycle become vaguer and may not be fully discriminated from each other for an LTD Stirling engine, they can operate autonomously without being controlled by external agents. Therefore, in Sec. II B, we introduce the dynamical equations of our engine model [6].

B. Dynamical equations
The set of equations that describe our LTD kinematic Stirling engine constitute the equation of motion of the power piston, equation of motion of the crank, and time-evolution equation of the gas temperature given as the energy conservation law (the first law of thermodynamics): Here, m p and Γ p in Eq. (3) are the mass and friction coefficient of the power piston, respectively. Further, F int in Eqs. (3)-(5) is the action-reaction force between the power piston and crank [5,15]. Γ and T load in Eq. (4) are the friction coefficient of the crank and load torque acting on the crank, respectively. p air in Eqs. (3) and (5) is the atmospheric pressure acting on the power piston. The rate of internal energy change of the gas on the left-hand side of Eq. (5) is equated to the heat fluxes and work flux on the right-hand side. The heat fluxes from the bottom and top surfaces of the large cylinder obey the Fourier law ( Fig. 1 (a)): where G is the thermal conductance associated with the heat transfer between the gas and surface of the large cylinder, and is a function that controls the coupling between the gas and bottom or top heat reservoir depending on the phase angle [6]. The role of the displacer piston transferring the gas into one side of the cylinders is represented by the function Eq. (8). Thus, we can revise where T eff (θ) is the effective temperature that periodically changes depending on the phase angle θ as We can thus consider the gas as though it were in contact with a single heat reservoir, the temperature of which dynamically oscillates in a sinusoidal manner between T b at θ = π 2 (χ b ( π 2 ) = 1 and χ t ( π 2 ) = 0) and T t at θ = 3π 2 (χ b ( 3π 2 ) = 0 and χ t ( 3π 2 ) = 1), which loosely approximates the ideal Stirling thermodynamic cycle [6].
We assume that the mass of the power piston and friction coefficient in Eq. (3) are negligible, as m p = Γ p = 0. We then obtain F int = σ p nRT V(θ) − p air from Eq. (3). By inserting this into Eqs. (4) and (5), and noting Eq. (9), we obtain Subsequently, Eqs. (12) and (13) are expressed in terms of the three-dimensional dynamical system as where ω denotes the angular velocity. By assuming a timescale separation between the crank and gas dynamics, we can make the adiabatic approximation dT dt = 0, by regarding T as a fast variable and θ and ω as slow variables [16]. By formally substituting dT dt = 0 into Eq. (16) and solving it with respect to T , we have the adiabatic approximated solution which is determined by the slow variables θ and ω of the crank (see Appendix A for detailed derivation). This approximation indicates that the motion of the piston and crank is considered as an externally controlled parameter for the gas, rather than being dynamically determined by the coupled equations in Eqs. (12) and (13) involving the gas dynamics. By substituting Eq. (17) into Eq. (15), we obtain the following twodimensional dynamical system: These dynamical equations describe the engine as a nonlinear pendulum driven by the temperature difference; the term that is proportional to sin 2 θ∆T included in the first term (rotational torque) on the right-hand side of Eq. (19) represents this driving force. Equations (18) and (19) (or Eqs. (14)- (16) before the adiabatic approximation) are the basic dynamical equations of our LTD kinematic Stirling engine model. The stationary and rotational states of the engine are described as a fixed point and stable limit cycle of Eqs. (18) and (19), respectively, which coexist depending on the parameters [6].
The stability of the fixed points on the thermodynamic and dead center branches are investigated by checking the determinant ∆ and trace T , calculated from the linearized equations of Eqs. (18) and (19) respectively, where θ * is given as the solution of Eq. (23). Figure 2 (a) shows the four branches for T load = 0, where the solid and dashed curves denote the stable fixed point and unstable fixed point (saddle point), respectively. For the given parameters, we have θ eq1 = π 4 and θ eq2 = 7π 4 , and we thus have the thermodynamic branches θ th1 and θ th2 satisfying θ th1 | (∆T,T load )=(0,0) = π 4 and θ th2 | (∆T,T load )=(0,0) = 7π 4 . In the vicinity of the equilibrium state θ eq , the thermodynamic branch θ th can be expanded as θ th ≃ θ eq +a 1Tload +a 2 ∆T , where a i are the expansion coefficients to be determined. By substituting this expansion into Eq. (23), we obtain For the present case of θ eq1 = π 4 and θ eq2 = 7π 4 , we can easily ∆T .
The linear response lines of θ th from the original equilibrium value θ eq are shown in Fig. 2 (a).

B. Heat fluxes at stationary states
For the non-vanishing ∆T , the engine conducts heat from the hot heat reservoir to the cold heat reservoir at the stationary states. The heat flux from each heat reservoir into the gas at the stationary state is given by with G 4 cos 2 θ * being an effective thermal conductance that depends on θ * . Figure 2 nRT eq nRTeq I on the stable thermodynamic branches θ * = θ th1 , θ th2 corresponding to those in Fig. 2 (a), where we can approximate J Q b (θ th1 ) and J Q b (θ th2 ) as in the vicinity of the equilibrium state, by using in Eq. (29).

IV. ROTATIONAL STATE
A. Numerical calculations of time-averaged angular velocity and heat fluxes We investigate the stable limit cycle of Eqs. (18) and (19) representing the rotational state of the engine. Denoting one cycle period of the stable limit cycle by τ, we define the timeaveraged angular velocity and heat fluxes as respectively, where · · · ≡ 1 τ τ 0 · · · dt denotes a time average and T (θ, ω) in Eq. (35) is given by Eq. (17).
In Fig. 3 (a), we present the ω -T load curve of the stable limit cycle. See also Fig. 3 (b) for the corresponding thermodynamic and dead center branches.
For sufficiently smallT load > 0, the engine is able to rotate against the load torque, producing positive work ( ω > 0). As T load increases, the engine stops rotating atT ′ load ≃ 7.0125 × 10 −5 , which is the bifurcation point of the stable limit cycle. AsT load further increases and exceeds the bifurcation point T ′′ load ≃ 9.9027 × 10 −5 , the stable limit cycle appears again; the engine is able to rotate again but in the same direction as the applied load torque ( ω < 0). ω shows the linear dependency withT load as it deviates sufficiently from the bifurcation points.
The above bifurcations are homoclinic bifurcations [7]. To illustrate this for the bifurcation atT ′ load , we show the orbit of the stable limit cycle on the phase plane in Fig. 4 (a) and the periodτ = 2π ω in Fig. 4 (b) in the vicinity ofT ′ load . In Fig. 4 (a), we can see that the orbit of the stable limit cycle closely passes the saddle point on the BDC branch in Fig. 3 (b) by taking a long time. At the bifurcation point, part of the orbit touches the saddle point and the stable limit cycle disappears, forming a homoclinic orbit [7]. Thus, although the dead center branch is not connected to the equilibrium state, the saddle point on the branch plays an important role in the homoclinic bifurcation of the limit cycle. As characteristics of the homoclinic bifurcation, the period of the limit cycle exhibits slow divergence according to the theoretical predictioñ τ ∝ − log(T ′ load −T load ) [7], which is confirmed in Fig. 4 (b). This slow divergence indicates a steep change in the angular velocity ω = 2π τ near the bifurcation points, as shown in Fig. 3 (a).
We show theT load dependence of J Q b in Fig. 5 (a). J Q b shows the linear dependency withT load as it deviates sufficiently from the bifurcation points in the same manner as ω in Fig. 3 (a). This linear dependency will be further investigated in Sec. IV D. AsT load approaches the bifurcation points, we find that J Q b deviates from the linear line and slowly converges to a constant value. This behavior is associated with the homoclinic bifurcation, which will be clarified in Sec. IV C.  24) and (25)). There are one or two stable fixed points, depending on the value ofT load .

B. Derivation of formal analytical expressions
We derive formal analytical expressions of the timeaveraged fluxes ω and J Q b for a small temperature difference and load torque, to explain their behaviors as we have seen in Sec. IV A. We first derive a formal analytical expression of ω using Eqs. (18) and (19). We assume that (θ, ω) is the stable limit cycle with period τ of Eqs. (18) and (19). Then, time-averaging both sides of Eq. (19) yields Note that the inertia term on the left-hand side has vanished as dω We can then approximate Eq. (17) as assuming that |∆T | and |ω| are sufficiently small. By using Eq. (37), we can rewrite the first term (rotational torque term) on the right-hand side of Eq. (36) as From Eqs. (36) and (38), we obtain where · · · θ ≡ 1 2π 2π 0 · · · dθ denotes a phase average. This formal analytical expression states that the averaged angular velocity is determined by the time average of the rotational torque and load torque.
We next derive a formal analytical expression of the time-averaged heat flux J Q b . Under the approximation of Eq. (37), the heat flux J Q b in Eq. (35) is approximated as By using τ = 2π ω and noting that ω = dθ dt , we obtain where ω is given in Eq. (39). The first term on the right-hand side of Eq. (41) is the time-averaged heat flux formally obeying the Fourier law, with G 4 cos 2 θ being the time-averaged thermal conductance. However, this is not similar to the heat leakage at the stationary state in Eq. (29) because of its strong correlation with the engine's rotational motion through the time-averaged thermal conductance. The second term on the right-hand side of Eq. (41) represents the heat transfer in proportion to the averaged angular velocity, which is also caused by the engine's rotational motion.

C. Near the bifurcation point
Near the bifurcation point, the orbit of the limit cycle stays in proximity to the saddle point almost all the time (Fig. 4 (a)). Thus, the effective thermal conductance G 4 cos 2 θ in Eq. (41) is approximated as G 4 cos 2 θ ≃ G 4 cos 2 θ H + a ω , where θ H of the saddle point (θ H , 0) on the BDC branch in Fig. 3 (b) is evaluated at the homoclinic bifurcation points and a is a coefficient that needs to be determined numerically. Equation (41) can be approximated in the vicinity of the bifurcation points as In Fig. 5 (b), Eq. (42) is compared with the numerical results for the bifurcation pointT ′ load . They are in good agreement and the linear decreasing from the constant value is confirmed.

D. Quasi-linear response regime
The angular velocity ω shows a linear dependency on T load as it deviates from the bifurcation point to a sufficient extent (Fig. 3 (a)). It also shows a similar linear dependency with respect to ∆T [6]. We call a regime with this linear dependency a quasi-linear response regime. In this regime, we may approximate ω by a constant value Ω as ω ≃ Ω by assuming that the periodic variation around the constant value is sufficiently small. Under this assumption, one cycle period is approximated as dt ≃ dθ Ω and thus τ = τ 0 dt ≃ 2π 0 dθ Ω = 2π Ω . Thus, the rotational torque component in Eq. (39) is approximated as where we used In Fig. 3 (a), the theoretical line and numerical calculations are compared, which are in good agreement. Next, we consider the heat flux Eq. (41) in the quasi-linear response regime using the result of Eq. (44). By approximating dt ≃ dθ Ω as above, we have G 4 cos 2 θ ≃ G 4 cos 2 θ θ ≃ G 8 . Then, the heat flux in Eq. (41) is approximated as where we used 2π 0 sin θ V(θ) dθ = 0. The theoretical line and numerical calculations show a good agreement (Fig. 5 (a)). Note that J Q t ≃ − J Q b , which can be confirmed by repeating the same calculations as J Q b . The theoretical expressions Eqs. (44) and (45) will be used for developing a theory of the thermodynamic efficiency of the engine in the quasi-linear response regime in Sec. V.

A. Definition of power and thermodynamic efficiency
We define the power and thermodynamic efficiency of the LTD Stirling engine. The instantaneous power produced by the gas, which is the second term on the right-hand side of Eq. (13), can be rewritten as where we used Eq. (19) from the second line to the third line. We can interpret each term in Eq. (46) as follows. The first term is the rotational kinetic energy change of the crank, and the second, third, and last terms represent the work carried out against the atmospheric pressure, frictional torque, and load torque, respectively. By using Eq. (46), we define the cycle-averaged power as where we used τ 0 d dt I 2 ω 2 dt = 0 and τ 0 p air dV dt dt = 0. The power P, which is referred to as the indicated power [15], defined as the closed area of the pressure-volume diagram of an engine, was decomposed into that carried out against the friction torque P fric and that carried out against the load torque P load , referred to as the brake power [15]. The former is eventually dissipated into the surrounding air as heat. By timeaveraging the energy conservation equation Eq. (13), we have J Q b + J Q t = P load + P fric . The thermodynamic efficiency η is then defined as the ratio of the input heat flux from the hot heat reservoir converted into the available power exerted against the load torque (brake power) P load [15]. For ∆T > 0, it is explicitly given as In Fig. 6 (a) and (b), we present the numerical results of thẽ T load dependence of the (nondimensionalized) brake power P load = P load nRT eq nRTeq I and the efficiency η, respectively. We can see that the values at which the maximum efficiency and maximum power are realized are close, which is characteristic of heat engines operating with non-negligible heat leakage (the first term on the right-hand side of Eq. (41) for the present model) [17]. When the maximum efficiency is located in the quasi-linear response regime, we can obtain its theoretical value, as we will show in Sec. V B.

B. Thermodynamic theory of an LTD kinematic Stirling heat engine in a quasi-linear response regime
Before constructing a thermodynamic theory of the LTD kinematic Stirling engine, we first review the theory for conventional linear irreversible heat engines [11][12][13].
For generic heat engines, the entropy production rate of the total systemσ (a heat engine and heat reservoirs at temperatures T 0 and T 1 ) is given bẏ where Q 0 (Q 1 ) is the heat flowing into the working substance from the heat reservoir at T 0 (T 1 ), and we used P load =Q 0 +Q 1 (the energy conservation law). Hereafter, the dot refers to quantities per unit time for steady-state heat engines or quantities averaged over one cycle period for cyclic heat engines. We can express P load as P load = Fẋ using an external force F and its conjugate fluxẋ. By taking the limit of a small temperature difference and small external force, we can approximatė σ asσ where the temperature difference and the averaged temperature are given as ∆T = T 0 − T 1 and T eq = T 0 +T 1 2 , respectively, for the present setup. Here, we defined the thermodynamic forces F i and their conjugate fluxes J i as and In linear irreversible thermodynamics, we assume the following linear relations between the thermodynamic fluxes and forces as where L i j are the Onsager coefficients with reciprocity L 12 = L 21 [9,10]. The use of Eqs. (53) and (54) enable us to rewrite Eq. (50) asσ Fromσ ≥ 0 for the arbitrary F 1 and F 2 (the second law of thermodynamics), we obtain the following restrictions on the Onsager coefficients L i j : Here, we define the coupling-strength parameter q as which should satisfy |q| ≤ 1 from the last inequality in Eq. (56). The meaning of q can be elucidated by rewriting the heat flux in Eq. (54) by using J 1 instead of F 1 as The case of |q| = 1 is an ideal condition known as the tightcoupling condition for which the heat flux J 2 is in proportion to the motion flux J 1 . For the non-tight-coupling case |q| 1, the non-vanishing heat leakage L 22 (1 − q 2 )F 2 arises from the simultaneous contact between the two heat reservoirs on the engine, which decreases the thermodynamic performance of the engine, as will be shown below. The power and thermodynamic efficiency are written using the thermodynamic fluxes and forces in Eqs. (53) and (54) as respectively, where we assume F 2 > 0. It is more convenient to express them in terms of J 1 instead of F 1 as using Eqs. (53) and (58). For the tight-coupling case |q| = 1, the quasistatic limit J 1 → 0 yields the vanishing power P load → 0 and the Carnot efficiency η → F 2 T eq = ∆T For the non-tight-coupling case |q| 1, J 1 that yields the maximum efficiency is obtained as the solution of ∂η ∂J 1 = 0 as which takes a finite value unlike the quasistatic limit J 1 → 0 for the tight-coupling case |q| = 1. The maximum efficiency then reads [12,13] which is a monotonic function of q.
Thus far, we have reviewed the theory for conventional linear irreversible heat engines. Returning to our model of the LTD kinematic Stirling engine, the linear response relations such as Eqs. (53) and (54) expanded from an equilibrium state with F 1 = 0 and F 2 = 0 do not hold. This is because the rotational state described as the limit cycle is not connected to the equilibrium state, and the linear dependency in Eqs. (44) and (45) holds only when the external forces deviate sufficiently far from the bifurcation points. Nevertheless, we can formally write the linear relations applied to these quasi-linear response regimes in terms of the thermodynamic fluxes and forces.
We identify each quantity used in the theory of the linear irreversible heat engines as T 0 = T b , T 1 = T t ,ẋ = Ω, F = T load , Q 0 = J Q b , andQ 1 = J Q t − P fric . Using these quantities, we can write the entropy production rate of the LTD kinematic Stirling engine in the quasi-linear response regime aṡ where we used P load = J Q b + J Q t − P fric (the energy conservation law), and the thermodynamic fluxes and forces are related via the following linear relations: where L ′ i j are the quasi-linear response coefficients. The prime notation is used to demonstrate that they are defined for the quasi-linear response regime. The use of the definitions of the thermodynamic fluxes and forces, and Eqs. (44) and (45), makes it possible to identify the quasi-linear response coefficients L ′ i j as Here, we can confirm that a symmetric relation holds as L ′ 12 = L ′ 21 . In Fig. 7, the quasi-linear response coefficients in Eq. (69) and the symmetric relation are numerically confirmed. At this point, the origin of this symmetry is not questioned and it will be elucidated in Sec. VI in terms of the (anti-)reciprocity of the Onsager kinetic coefficients. Because the quasi-linear response relations Eqs. (67) and (68) with the symmetric relation formally take the same form as the conventional Onsager relations Eqs. (53) and (54), the thermodynamic theory developed using Eqs. (53) and (54) are also applied to the quasi-linear response regime.
In the present case, the coupling-strength parameter q in Eq. (57) is calculated from the quasi-linear response coefficients in Eq. (69) as Notably, the coupling strength depends on three major (nondimensionalized) physical parameters of the model:σ,G, and Γ. Thus, the maximum efficiency is given by Eq. (64) with the coupling strength q in Eq. (70) being the single figure of merit. In Fig. 6 (a) and (b), we compare the numerical results of the efficiency and power with the theoretical results Eqs. (62) and (61) using L ′ i j in Eq. (69). We can find that the theory approximates the numerical results well. Although we used Eq. (61) for the calculations of P load , it should be consistent with the expression P load = J Q b + J Q t − P fric (the energy conservation law), which was used in the derivation of Eq. (66). See Appendix B for a detailed demonstration of the equivalence of these two expressions. The maximum efficiency in Eq. (64) using q ≃ 0.17513 calculated for the present parameters also approximates the numerical result well (Fig. 6  (a)).
The simple formula Eq. (64) using Eq. (70) may provide a new guiding principle for designing efficient LTD kinematic Stirling engines. By noting we obtain q → 1 √ 2 as the upper bound of q in Eq. (70) as σ → 0 andGΓ → 0 withGΓ ≪σ 2 being satisfied. Within this limit, η max in Eq. (64) is given as This is the upper bound that the present model in the quasilinear response regime can attain. We note that η max of the present model cannot attain the Carnot efficiency achieved by the ideal Stirling cycle because it lacks a regenerator. We can also obtain lim q→ 1 as the upper bound of the efficiency at maximum power in Eq. (65) that the present model in the quasi-linear response regime can attain.

VI. ORIGIN OF SYMMETRIC RELATION
The symmetric relation L ′ 12 = L ′ 21 in Eq. (69) is reminiscent of the Onsager reciprocity in linear irreversible thermodynamics, whereas the rotational state of the engine described as the limit cycle may not be described as a linear response regime. Here, we explain the origin of the symmetry in terms of (anti-)reciprocity of the Onsager kinetic coefficients [9,10,19] in the original three-dimensional dynamical model Eqs. (14)-(16) before the adiabatic elimination.

A. Relaxation dynamics towards the equilibrium state
Let us consider that a mesoscopic LTD Stirling heat engine specified by (θ, p θ , U) is in thermal equilibrium with a heat reservoir, where p θ ≡ Iω is the angular momentum of the crank and U = f 2 nRT is the internal energy of the gas. The engine may be perturbed from the equilibrium state (θ eq , 0, U eq ) by thermal fluctuation and relaxes to the original equilibrium state, where U eq = f 2 nRT eq . By linearizing Eqs. (14)-(16) with ∆T = 0 and T load = 0 around the equilibrium value, we obtain the following linear relaxation equations: These are rewritten as (k, l = θ, p, U) where x θ ≡ δθ, x p ≡ δp θ , and x U ≡ δU are the thermodynamic variables that express variation (or fluctuation) from the equilibrium state, and λ kl are the linear relaxation coefficients.
Next, we expresess Eq. (78) (Eqs. (75)-(77)) as where X l are the conjugate thermodynamic forces to be determined, and γ kl are the Onsager kinetic coefficients. Following the methods in [19], we introduce δH as the change of the crank's Hamiltonian from the vanishing value at the equilibrium state (θ eq , 0, U eq ) as We then define X θ and X p as the thermodynamic forces for the mechanical degrees of freedom as where we can interpret X θ as a restoring force and X p as an inertial force. Under these thermodynamic forces, we can easily find which satisfy the Onsager's anti-reciprocal relation γ θp = −γ pθ . We note that the anti-reciprocity is fundamentally derived from the fact that x θ is a time-reversely symmetric quantity, whereas x p is an anti-symmetric quantity under time reversal of microscopic dynamics [19]. We also find γ U p = γ U p (θ eq ) = nRT 2 eq rσ p sin θ eq V(θ eq ) .
Once X θ and X p have been determined as above, X U , the other thermodynamic force of the thermodynamic degree of freedom, can be uniquely determined such that it satisfies the Onsager symmetry principle [19]. Because we want to have the anti-reciprocal relation γ U p = −γ pU for x U as a time-reversely symmetric quantity, we naturally choose X U as which determines γ pU = γ pU (θ eq ) = − nRT 2 eq rσ p sin θ eq V(θ eq ) , The other kinetic coefficients vanish as γ θθ = γ θU = γ Uθ = 0. We note that the thermodynamic variables x l and the forces X k are linearly related from Eqs. (78) and (79) as where The only non-vanishing components of β kl are the diagonal elements as β θθ = nRr 2 σ 2 p sin 2 θ eq V 2 (θ eq ) , (94) The total entropy variation (the heat engine and reservoir) around the maximum, equilibrium value is then approximated as the following quadratic form [19]: This variation is consistent with equilibrium statistical mechanics; fluctuation of temperature and volume of a system under isothermal and isobaric conditions obeys the following probability distribution [19]: where C V and k B are the constant-volume heat capacity of the system and the Boltzmann constant, respectively. The momentum of the system also fluctuates around its equilibrium value δp θ = 0 according to the Maxwell distribution: According to Einstein's fluctuation formula, the total weight W = w(δT, δV)w(δp θ ) is given using entropy variation δS as V 2 , and δV = ∂V ∂θ δθ = rσ p sin θδθ in Eq. (96), we find that δS in Eq. (98) agrees with that in Eq. (95). We note that we can define the thermodynamic force using δS as The instantaneous entropy production rate is thus given as to which the terms with the anti-reciprocal coefficients of γ kl do not contribute. Due to the anti-reciprocal component, the relaxation dynamics in the vicinity of the equilibrium state show damping oscillation toward the equilibrium state [20].
B. Expression of quasi-linear response coefficients using the Onsager kinetic coefficients Equations (75)-(77) describe the relaxation dynamics when the engine slightly deviates from the equilibrium state. For a nonequilibrium condition in which the externally sustained thermodynamic forces ∆T 0 and T load 0 are applied, the situation can drastically change. The engine can show rotational motion, and an engine under this state cannot be regarded as being in the linear response regime as we have seen in Sec. IV. Nevertheless, we will examine how the antireciprocity of the Onsager kinetic coefficients γ kl included in the relaxation dynamics is inherited by the symmetric response coefficients L ′ i j . We recall that we have adiabatically eliminated T from the three-dimensional dynamical model Eqs. (14)-(16) by assuming that the dynamics of the gas are fully subject to those of the crank, and obtained the two-dimensional dynamical model Eqs. (18) and (19). The quasi-linear relations in Sec. IV have been formulated for the rotational state of the two-dimensional dynamical model with ∆T 0 and T load 0. We now rewrite the thermodynamic fluxes J 1 = Ω = p θ I and J 2 = J Q b of the two-dimensional model in a form that highlights the relation to the Onsager kinetic coefficients γ kl included in the relaxation dynamics of the three-dimensional dynamical model for ∆T = 0 and T load = 0.
For the rotational state realized under the nonequilibrium condition, the engine is largely perturbed from the equilibrium state (θ eq , 0, U eq ). This is, however, only with respect to the phase angle θ. Because the deviations of p θ and U (or equivalently T ) from their equilibrium values are small even for the rotational state for small |∆T | and |T load |, we write p θ ≃ δp θ and T ≃ T eq + δT . We can thus expand Eq. (16) in terms of δp θ and δT around their equilibrium value (p θ , T ) = (0, T eq ), with θ being held fixed as an arbitrary value as Equivalently, we have where we denote by γ kl (θ) the Onsager kinetic coefficients with θ eq being formally replaced with θ of the stable limit cycle (θ, ω). For ∆T = 0 and T load = 0, no stable limit cycle exists and Eq. (102) recovers the relaxation dynamics Eq. (77) with θ = θ eq . The adiabatic approximation solution T (θ, δp θ ) = T eq + δT (θ, δp θ ) of Eq. (101) satisfying dT dt = dδT dt = 0 is given as Equivalently, from Eq. (102), we have the adiabatic approximation solution as using the Onsager kinetic coefficients. We next expand Eq. (15), which describes the rotational state in terms of δp θ and δT around their equilibrium value (p θ , T ) = (0, T eq ), with θ being held fixed as an arbitrary value as Equivalently, we can rewrite Eq. (105) as in terms of the thermodynamic forces. By putting X U in Eq. (104) into that in Eq. (106), we obtain Alternatively, by noting that T b − T (θ, δp θ ) = ∆T − δT (θ, δp θ ) and using Eq. (103), we can also rewrite the instantaneous heat flux We assume that the angular velocity δp θ I in the quasi-linear response regime is a constant as δp θ I = T eq X p ≃ Ω, in a similar manner as we have assumed in Sec. IV D. By taking a time average of Eqs. (107) and (108), and repeating essentially the same calculations as in Sec. IV D, we have By putting Eq. (109) into Eq. (110), we obtain Finally, from Eqs. (109) and (111), the quasi-linear response coefficients are found to be which are given using the phase averages · · · θ of the quantities that include the Onsager kinetic coefficients using θ instead of θ eq . By performing the phase averages in Eq. (112), we can confirm that Eq. (112) agrees with Eq. (69). From the expression of Eq. (112), we immediately notice that the symmetric relation L ′ 12 = L ′ 21 in the adiabatically eliminated model holds as a consequence of the anti-reciprocal relation of the Onsager kinetic coefficients γ pU (θ) = −γ U p (θ) included in the three-dimensional dynamical model before the adiabatic elimination. Recalling that the anti-reciprocity of the Onsager kinetic coefficients reflects the time-reversal symmetry of the underlying microscopic dynamics [19], the present symmetric relation may also be attributed to the time-reversal symmetry. Although the anti-reciprocal terms do not contribute to the instantaneous entropy production rate Eq. (100) during the relaxation dynamics, they can contribute to the entropy production rate averaged over one cycle period for the rotational state in the quasi-linear response regime through L ′ i j (Eq. (55)). In-terestingly, the restrictions on L ′ i j in Eq. (56) imposed by the second law of thermodynamics are also assured by this antireciprocity.

VII. SUMMARY AND DISCUSSION
This paper presented the nonequilibrium thermodynamics of a nonlinear dynamics model of an LTD kinematic Stirling heat engine [6]. The two-dimensional dynamical equations describing the crank of the engine were derived from the original three-dimensional dynamical equations based on the adiabatic elimination of the gas dynamics. By using the twodimensional dynamical equations, we investigated the stationary and rotational states, which are the fixed points and stable limit cycle of the equations, respectively. In particular, we focused on the regime near the bifurcation points and the quasilinear response regime sufficiently far from the bifurcation points of the latter state. The formal analytical expressions of the averaged angular velocity and heat fluxes (thermodynamic fluxes) as a function of temperature difference and load torque (thermodynamic forces) were derived to explain these regimes in the rotational state. In the quasi-linear response regime, it was found that the thermodynamic fluxes and forces are described by the linear relations with symmetric coefficients. Based on the linear relations, we obtained the maximum efficiency formula in terms of the coupling-strength parameter as the single figure of merit. We also demonstrated that the symmetric coefficients are considered as a consequence of the anti-reciprocal relation of the Onsager kinetic coefficients in the relaxation dynamics before the adiabatic elimination.
Irrespective of whether the engine operates with external agents, such as conventional heat engines, or autonomously, such as in the present model, we analyzed their thermodynamic performance on an equal footing based on the linear relations [11][12][13]. When the adiabatic elimination is not valid, the dynamics of the gas and piston-crank system are not separated, and they constitute a dynamical system as a whole. When the adiabatic elimination is valid, the dynamics of the gas are completely subject to those of the piston-crank system. This yields explicit separation between the system and external agents. In this sense, there may not be much difference between conventional periodically driven heat engines operated by external agents and the present LTD kinematic Stirling engine, although the dynamics of the external operator itself in the latter case obeys the equations of motion. However, for the present self-sustained engine we have observed that the emergence of the symmetric coefficients can be explained based on the property of the relaxation dynamics towards the equilibrium state before the adiabatic elimination. This demonstrates the importance of modeling an autonomous heat engine as a dynamical system with mechanical and thermodynamic degrees of freedom. The symmetric relation can be experimentally verified in principle. It is of interest to investigate the similarities and differences between the present emergent symmetry in the quasi-linear response regime and other various symmetries found in periodically driven heat engines operated by external agents in the linear response regime [21][22][23][24][25][26].
Our theory may be useful for predicting the possible future design of efficient LTD Stirling heat engines by employing our maximum efficiency formula. In this sense, although our theory is expected to describe existing LTD Stirling heat engines, it could also be used to describe more advanced engines in the near future.