Thermodynamic Analysis of Bistability in Rayleigh–Bénard Convection

Bistability is often encountered in association with dissipative systems far from equilibrium, such as biological, physical, and chemical phenomena. There have been various attempts to theoretically analyze the bistabilities of dissipative systems. However, there is no universal theoretical approach to determine the development of a bistable system far from equilibrium. This study shows that thermodynamic analysis based on entropy production can be used to predict the transition point in the bistable region during Rayleigh–Bénard convection using the experimental relationship between the thermodynamic flux and driving force. The bistable region is characterized by two distinct features: the flux of the second state is higher than that of the first state, and the entropy production of the second state is lower than that of the first state. This thermodynamic interpretation provides new insights that can be used to predict bistable behaviors in various dissipative systems.


Introduction
Many attempts have been made to identify a universal function whose extremum determines the development of a system far from equilibrium. Entropy production characterizes systems during nonequilibrium processes [1][2][3][4][5], and its extrema may be used for determining the system behavior [6][7][8][9][10][11][12][13]. Recently, theoretical and experimental investigations have suggested the maximization of entropy production during nonequilibrium processes (the so-called maximum entropy-production principle, MEPP) [10][11][12][14][15][16]. According to MEPP, when a nonequilibrium system transitions from one state to another, it is characterized by the highest rate of entropy production. Analysis based on MEPP can thus be used to predict the transition point between two nonequilibrium states, such as those observed in the morphological changes of crystal growth, mode changes in droplet oscillation, and pattern changes in thermal convection [12,[17][18][19]. For example, in the case of friction phenomena in the flow field, fluid velocity is treated as thermodynamic flux, and the transition point is predicted. For a pendant droplet that changes in the oscillation mode induced by the solutal Marangoni effect with viscous dissipation, the transition point of the oscillation mode is predicted from the intersection of the entropy-production curves determined from the velocity of the oscillating droplet, which is considered as thermodynamic flux [18]. However, dissipative systems far from equilibrium frequently include solutions with several linearly stable branches, i.e., bistable behavior. In such cases, the selected solution depends on the initial conditions, and variational principles based on MEPP would not be required. The prediction of bistable behavior in various dissipative systems is considered an unsolved problem when using variational principles and overshadows the universality of MEPP. This problem is addressed in this study by examining these predictions for a situation involving bistable behavior, i.e., where hexagonal and roll flow patterns coexist during Rayleigh-Bénard convection.

MEPP
The state of a nonequilibrium system is characterized by the thermodynamic flux J expressed as a function of the driving force of the entire system F, which is proportional to differences in the temperature, concentration, pressure, etc. [12,[17][18][19]: where a i and b i are constant coefficients. Assume that there is a transition from state α to β as F increases. Each state is described by a local entropy-production σ curve, which is characterized by the product of the local thermodynamic force X and thermodynamic flux J [1]. Thus, the relationship between σ and F can be expressed as [12,[17][18][19] According to MEPP, a nonequilibrium system develops in a manner that maximizes its entropy production under the given binding conditions [10,14]. The transition point between states α and β corresponds to the intersection of the two σ curves A and not that of the two J lines B. Although the coefficients of F in σ i are not strictly equal to a i in Equation (1), there is no change in the intersection. Thus, for simplicity, we use the same coefficients. Before intersection A, the entropy-production curve of state α lies above that of state β, whereas after intersection A, the converse is true. It should be noted that when entropy production is expressed as a function of X, σ(X) cannot distinguish between the nonequilibrium states because all σ(X) are present on a single quadratic curve σ(X) = LX 2 [19], where L is the phenological coefficient [1]. Furthermore, we cannot predict the transition point nor even understand whether another kind of nonequilibrium state exists in the system when the thermodynamic flux expressed as a function of the driving force can be described only on a single line. Even at equilibrium, if we know the thermodynamic properties of water only in the liquid phase at 1 atm, we cannot predict that water will undergo a phase transition to the gas phase at 100 • C. If we know the dependence of the chemical potential of water in both phases on temperature, we can predict the boiling point of water from the intersection of their two chemical-potential curves.
As shown in Figure 1, we find that intersection A represents the transition point between the nonequilibrium states. However, we wish to understand intersection B of the thermodynamic-flux lines in terms of the physical behavior of the nonequilibrium system. It transpires that this intersection is a starting point for the coexistence of two different nonequilibrium states. When F > F B , state β begins to manifest because the thermodynamic flux of state β is higher than that of state α. However, the system mainly comprises state α because σ(F) of state α remains greater than that of state β. Therefore, states α and β coexist until F < F A . The thermodynamic flux of state β increases continuously with F, so the proportion of state β will increase until it represents a major part of the system. Once F > F A , the system consists only of state β.

MEPP 43
The state of a nonequilibrium system is characterized by the thermodynamic flux expressed as a function 44 of the driving force of the entire system , which is proportional to differences in the temperature, 45 concentration, pressure, etc. [12,[17][18][19]: 46 where and are constant coefficients. 47 Assume that there is a transition from state to as increases. 48

51
Each state is described by a local entropy-production curve, which is characterized by the product of 52 the local thermodynamic force and thermodynamic flux [1]. Thus, the relationship between and F can 53 be expressed as [12,[17][18][19]

Thermodynamic Analysis of Rayleigh-Bénard Convection
To verify this MEPP prediction, precise experimental data are required to calculate the relationship between the thermodynamic flux and driving force of a dissipative system exhibiting bistability. One example is reported in the literature [20], where hexagonal and roll patterns coexisted during Rayleigh-Bénard convection, subject to external temporal modulation of the reduced Rayleigh number ε 0 ≡ ∆T/∆T c − 1, where ∆T is the temperature difference between the bottom and top plates of the water-filled container and ∆T c is the critical temperature difference for the onset of convection without modulation. Here, ε has the form where the time t and frequency ω are scaled according to the vertical thermal-diffusion time and δ is the amplitude of modulation. The 13-mode Lorenz model proposed by Hohenberg and Swift predicts a positive threshold-shift change in the convection onset from ε 0 = 0 to ε 0 = ε c [21]. Above ε c , the roll patterns that appear through supercritical bifurcation are unstable to the hexagonal patterns (reproduced in the inset of Figure 2a) [21]. This continues as ε 0 increases until ε 0 = ε R , beyond which the roll patterns are stable. Hexagonal patterns manifest through subcritical bifurcation, first becoming stable at ε 0 = ε A < ε c , which continues until ε 0 = ε B > ε R , where they are unstable to roll patterns. For ε R < ε < ε B , both hexagonal and roll patterns are stable. Meyer reported that the bistable region for the range of loop ε B − ε R was approximately two orders of magnitude larger than that for the loop between ε A and ε C [22]. Therefore, it is easy to resolve the bistable region experimentally.

107
Heat flux conv is a dimensionless quantity given by the ratio of the convective heat flux to the heat flux 108 conducted through the fluid. When > 0 , convective flow occurs and conv becomes positive. Meyer, 109 Cannell, and Ahlers performed experimental observations of Rayleigh-Bénard convection and the heat flux, as 110 shown in Figure 2a [20]. They detected the bistable region where both hexagonal and roll patterns were stable. 111 The model quantitatively and qualitatively agrees with the experimental results in the pure-hexagonal ( c ≤ 112 0 < R ), bistable ( R ≤ 0 < B ), and pure-roll ( 0 > B ) regions for = 1.97, where c = 0.132, R = 0.244, 113 and B = 0.350. 114 Let us begin by predicting the transition point between heat conduction and heat convection using MEPP 115 with static measurements without modulation. This is easy to predict because of the linear relationship 116 between conv and 0 in each state.  [20]. They detected the bistable region where both hexagonal and roll patterns were stable. The model quantitatively and qualitatively agrees with the experimental results in the pure-hexagonal (ε c ≤ ε 0 < ε R ), bistable (ε R ≤ ε 0 < ε B ), and pure-roll (ε 0 > ε B ) regions for δ = 1.97, where ε c = 0.132, ε R = 0.244, and ε B = 0.350.
Let us begin by predicting the transition point between heat conduction and heat convection using MEPP with static measurements without modulation. This is easy to predict because of the linear relationship between J conv and ε 0 in each state.
Here, ε 0 corresponds to the driving force of the entire system F. The values of J conv may be linearly fitted as a function of ε 0 in each state. The heat flux for heat conduction is obviously zero. The best-fit line for J conv of heat convection with ε 0 > 0 is J conv = 1.29 ε 0 − 9.42 × 10 −4 . The entropy production of the heat-conduction and heat-convection regions is easily calculated from Equation (2), where σ = 0 and σ = 1.29 ε 0 − 9.42 × 10 −4 2 , respectively. Thus, we obtain the intersection ε 0 = 9.42 × 10 −4 ≈ 0. It is clear that the transition occurs at ε 0 = 0, because of the definition of heat flux wherein the positive value of J conv signifies the occurrence of heat flux produced by only convective flow; however, it is important that this is predicted using MEPP.
Next, we analyze Rayleigh-Bénard convection subject to external temporal modulation on the MEPP basis. The transition points cannot be distinguished easily because the values of J conv show rounded bifurcations due to modulation from heat conduction to convection; for the heat-convection Entropy 2020, 22, 800 5 of 8 region (ε 0 > ε C ), they align approximately along a single line that changes only slightly in slope. The derivative of J conv with respect to ε 0 enables recognizing the point of change in the slope (Figure 3a). fitted as a function of 0 in each state. The heat flux for heat conduction is obviously zero. The best-fit line for 125 conv of heat convection with 0 > 0 is conv = 1.29( 0 − 9.42 × 10 −4 ) . The entropy production of the 126 heat-conduction and heat-convection regions is easily calculated from Equation (2), where = 0 and = 127 1.29( 0 − 9.42 × 10 −4 ) 2 , respectively. Thus, we obtain the intersection 0 = 9.42 × 10 −4 ≈ 0. It is clear that the 128 transition occurs at 0 = 0, because of the definition of heat flux wherein the positive value of conv signifies 129 the occurrence of heat flux produced by only convective flow; however, it is important that this is predicted 130 using MEPP. 131 Next, we analyze Rayleigh-Bénard convection subject to external temporal modulation on the MEPP 132 basis. The transition points cannot be distinguished easily because the values of conv show rounded 133 bifurcations due to modulation from heat conduction to convection; for the heat-convection region ( 0 > C ), 134 they align approximately along a single line that changes only slightly in slope. The derivative of conv with 135 respect to 0 enables recognizing the point of change in the slope (Figure 3a). For 0.13 < ε 0 < 0.21, dJ conv /dε 0 increases monotonically with ε 0 , whereas for ε 0 > 0.21, it remains approximately constant. Thus, we can divide the heat-convection region in two based on the change in dJ conv /dε 0 . The values of J conv in the two regions can be linearly fitted as functions of ε 0 , yielding J conv = 1.108(ε 0 − 0.116) and J conv = 1.197(ε 0 − 0.125) (Figure 3b). The former corresponds to the hexagonal convection and the latter to roll convection. Using the method described above yields two curves for entropy production of the hexagonal and roll flows, respectively: σ = 1.108(ε 0 − 0.116) 2 and σ = 1.197(ε 0 − 0.125) 2 . These curves intersect at one point, ε 0 = 0.354, even though the two curves appear to overlap.
This intersection corresponds to a transition point between the hexagons and rolls, and it is in good agreement with ε B = 0.350, where the rolls become stable to the hexagonal patterns in both the experimental and theoretical results. Furthermore, the intersection of the two J conv lines occurs at ε 0 = 0.238. The value of ε 0 is close to a starting point for the bistable region, ε R = 0.244, where both hexagonal and roll patterns are stable in the experimental and theoretical results. The actual patterns obtained by Meyer et al. are shown in Figure 4.
For ε 0 > ε R , the roll patterns gradually overlap with the stable hexagonal patterns, and coexisting patterns persist until ε 0 = ε B . As shown in Figure 3b, the J conv of the rolls lies above that of the hexagons, whereas σ of the rolls lies below that of the hexagons. For ε 0 > ε B , J conv and σ for the rolls lie above those for the hexagons, and the actual patterns show that only the rolls are stable.
The experimental values of J conv represent the total heat flux produced by the three directional components of heat flux, i.e., x, y, and z. The vertical component, which is aligned with the temperature difference ∆T between the bottom and top plates, accounts for a large proportion of the total heat flux; hence, it obscures the slight change produced by the horizontal components [19]. Thus, the experimental results do not show the jump in J conv , as indicated by the schematic in Figure 1. If the effect of the vertical component is removed from the total flux and the heat flux is analyzed only in the direction perpendicular to ∆T, the jump in J conv can be observed, and more accurate transition points may be obtained.

176
For > , the roll patterns gradually overlap with the stable hexagonal patterns, and coexisting patterns 177 persist until = . As shown in Figure 3b, the of the rolls lies above that of the hexagons, whereas 178 of the rolls lies below that of the hexagons. For > , and for the rolls lie above those for the 179 hexagons, and the actual patterns show that only the rolls are stable. 180 The experimental values of represent the total heat flux produced by the three directional 181 components of heat flux, i.e., x, y, and z. The vertical component, which is aligned with the temperature 182 difference between the bottom and top plates, accounts for a large proportion of the total heat flux; hence, 183 it obscures the slight change produced by the horizontal components [19]. Thus, the experimental results do 184 not show the jump in , as indicated by the schematic in Figure 1. If the effect of the vertical component is 185 removed from the total flux and the heat flux is analyzed only in the direction perpendicular to , the jump 186 in can be observed, and more accurate transition points may be obtained 187

Discussion 188
The state with the highest entropy is the state where intensive variables are uniform in the entire system. 189 However, the dissipative system applied by external forces, such as difference in temperature, concentration, 190 and pressure fields, never develops into a uniform state. Under the given binding conditions, in order for the 191 dissipative system to approach more rapidly the state with a uniform field, heat, molecules, and the fluid 192 momentum must transfer faster. In Rayleigh-Bénard convection, the heat transfer of the system changes from 193 heat conduction to heat convection in a hexagonal pattern, and then to heat convection in a roll pattern, and 194 finally to turbulence as the difference in temperature increases. The state of heat transfer changes such that the 195 system more rapidly approaches the uniform temperature field. Therefore, the system should realize the state 196 with higher damping and higher energy consumption by changing the flow state to approach the uniform field 197 more rapidly.

198
The accuracy of the transition point predicted from our thermodynamic approach depends on the precision 199 of the data on the relationship between the thermodynamic flux and driving force. Thus, the predicted 200 transition point necessarily involves some uncertainty. However, as shown in Figures 2 and 3, in Rayleigh-201 Bénard convection, the state with the highest entropy production is more stable than a state with lower entropy 202 production. Bistability occurs when a new state has higher thermodynamic flux and the existing state has higher 203

Discussion
The state with the highest entropy is the state where intensive variables are uniform in the entire system. However, the dissipative system applied by external forces, such as difference in temperature, concentration, and pressure fields, never develops into a uniform state. Under the given binding conditions, in order for the dissipative system to approach more rapidly the state with a uniform field, heat, molecules, and the fluid momentum must transfer faster. In Rayleigh-Bénard convection, the heat transfer of the system changes from heat conduction to heat convection in a hexagonal pattern, and then to heat convection in a roll pattern, and finally to turbulence as the difference in temperature increases. The state of heat transfer changes such that the system more rapidly approaches the uniform temperature field. Therefore, the system should realize the state with higher damping and higher energy consumption by changing the flow state to approach the uniform field more rapidly.
The accuracy of the transition point predicted from our thermodynamic approach depends on the precision of the data on the relationship between the thermodynamic flux and driving force. Thus, the predicted transition point necessarily involves some uncertainty. However, as shown in Figures 2 and 3, in Rayleigh-Bénard convection, the state with the highest entropy production is more stable than a state with lower entropy production. Bistability occurs when a new state has higher thermodynamic flux and the existing state has higher entropy production. It is a significant achievement that Rayleigh-Bénard convection cannot develop against the above rules.

Conclusions
This study shows that thermodynamic analysis based on entropy production can be used to predict the transition point in the bistable region, i.e., the region where hexagonal and roll flow patterns coexist during Rayleigh-Bénard convection. This addresses a gap in our understanding of MEPP with respect to dissipative systems far from equilibrium. These systems frequently require solutions with several linearly stable branches, and in such cases, the selected solution should depend on the initial conditions. This work may be used to predict bistable dissipative-system behavior in a wide range of applications, such as biological switching [23][24][25], optical switching [26,27], and chemical switching [28][29][30].