{'fraction_non_alphanumeric': 0.035530720366904765, 'fraction_numerical': 0.013337199188170484, 'mean_word_length': 4.484244001156403, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interfacial instability induced by lateral vapor pressure fluctuation in bounded thin liquid-vapor layers

We study an instability of thin liquid-vapor layers bounded by rigid parallel walls from both below and above. In this system, the interfacial instability is induced by lateral vapor pressure fluctuation, which is in turn attributed to the effect of phase change: evaporation occurs at a hotter portion of the interface and condensation at a colder one. The high vapor pressure pushes the interface downward and the low one pulls it upward. A set of equations describing the temporal evolution of the interface of the liquid-vapor layers is derived. This model neglects the effect of mass loss or gain at the interface and guarantees the mass conservation of the liquid layer. The result of linear stability analysis of the model shows that the presence of the pressure dependence of the local saturation temperature mitigates the growth of long-wave disturbances. The thinner vapor layer enhances the vapor pressure effect. We find the stability criterion, which suggests that only slight temperature gradients are sufficient to overcome the gravitational effect for a water/vapor system. The same holds for the Rayleigh-Taylor unstable case, with a possibility that the vapor pressure effect may be weakened if the accommodation coefficient is below a certain critical value.


I. INTRODUCTION
Thin liquid films have been intensively studied over the last decades. Many contributions have been devoted to them, owing to their technological importance and wide industrial applications. Their rich interfacial behaviors originate from combinations of various effects such as capillarity, intermolecular forces, thermocapillarity and gravity. One important effect among them is evaporation or condensation. It was often incorporated in the studies of thin films, 1,2 and dewetting patterns resulting from drying of films were analyzed. 3,4,5 Many past studies on evaporating or condensing liquid films have ignored the dynamics of the gas above the liquid film, assuming the infinitely deep gas phase. 1,2,6 However, in this study we consider a situation where the gas phase is bounded by rigid parallel wall from above and has a finite depth comparable with the liquid one. The full linear stability analyses of this system were performed in several papers. 7,8,9,10,11 Despite the apparent simplicity of the configuration, this system includes a free surface and interfacial boundary conditions involving phase change, and therefore is very sophisticated. In order to simplify this problem, we apply long-wave approximation to both layers.
The advantage of the use of long-wave or lubrication approximation is the reduction of dimensionality: a one-dimensional (two-dimensional) film evolution equation can be derived in a two-dimensional (three-dimensional) system. Normally, only the dynamics of the liquid is considered, leading to a one-sided model. However, if the ambient gas layer is thin enough, a two-layer model would better describe the system. This was demonstrated by VanHook et al., 12 who developed a two-layer theory to reproduce their experimental results. They showed that their two-layer model better predicts the onset of instability in their experiment than the corresponding one-layer model and also correctly describes the formation of localized elevations. In their approach, only the heat conduction in the gas phase is taken into account, and the gas dynamics is ignored because the viscosity of the gas is much less than that of the liquid. Later, Merkt et al. 13 presented an evolution equation of the interface of two viscous fluid layers in the same geometry. Their model allows for the shear stress induced by the motion of the upper layer and therefore is reduced to the single layer equation in the limit of small viscosity of the upper layer. Although their goal is the observation of pattern formation in the long-time regime, the two-layer models have been applied to the cases of Rayleigh-Taylor instability 14,15 and ultrathin films. 16,17 Nevertheless, in the two-layer systems mentioned above there is no phase transformation at the interface. Here, we construct a two-layer theory for liquid-vapor layers which undergo phase change, using long-wave approximation. Note that the application of long-wave theory to the vapor phase was made in the study of film boiling. 18 If we take into account the effect of the mass flux across the interface, an instability peculiar to this system is expected, even for the presence of large disparity in viscosity and density between liquid and vapor; see The liquid film is initially in equilibrium with its vapor layer. If the liquid side is heated or the vapor side is cooled, evaporation occurs at a hotter portion of the interface and condensation at a colder one. Since the vapor layer is bounded, the vapor pressure becomes higher in the evaporating region and lower in the condensing one. According to this lateral vapor pressure gradient, the higher vapor pressure pushes the interface downward and the lower one pulls it upward. Then, the surface deflection is amplified. To our knowledge, this pressure-induced instability mechanism has not been considered in the past, because the uniform ambient vapor pressure has been assumed in previous studies of evaporating or condensing liquid films. 1,2,6 To derive the model, we require the interfacial boundary conditions such that the mass transfer occurs between the phases. We follow those of the earlier studies on evaporating or condensing liquid films 1,2 except the thermodynamic relation at the interface. They used the linearized equation where the mass flux through the interface is proportional to the difference between the interfacial temperature and its saturation value corresponding to the surrounding vapor pressure, based on kinetic theory. 19 This relation cannot be directly applied to the present problem, because the local saturation temperature varies in the lateral direction, depending on the vapor pressure. Hence, we must modify the relation to take this effect into account. Fortunately, this can be easily done in the thermodynamic framework.
For instance, Ajaev and Homsy 20 and Wayner 21 used a nonequilibrium thermodynamic relation including the saturation temperature variation due to capillarity and disjoining pressure. This effect was later included in the model of evaporating or condensing thin liquid films. 22 However, they assumed the vapor pressure to be constant. In our two-layer model, this relation should be extended in accordance with the vapor pressure variation.
Thus, one of the purposes of this work is to investigate the effect of the vapor pressure dependence of the local saturation temperature. Note that this effect was not considered in the film boiling case, 18 although the lateral vapor pressure variation, which drives the motion of the vapor, certainly exists in the boiling film. It is worthwhile noting that a somewhat similar motivation to ours is found in Ref. 23, where the vapor concentration in the ambient gas phase above the liquid film fluctuates and thereby the mass flux varies along the interface. However, they neglect the bulk gas dynamics itself and consider only the diffusion of the vapor.
Here, we start with a more general nonequilibrium thermodynamic law, which reduces to, in the linear domain, a proportional connection between the interfacial mass flux and the difference of chemical potential in each phase. 24 From this law, we can naturally derive a thermodynamic relation similar to that of Ajaev and Homsy and Wayner. Moreover, the condition of local thermodynamic equilibrium, adopted in some papers on two-phase problem, 8,10,11,25 is recovered by taking the appropriate limit of the derived relation. Therefore, the thermodynamic relation used here is also the extension of the interfacial equilibrium condition into nonequilibrium states. We note that more general formulation taking into account the nonequilibrium effect contains a temperature discontinuity at the liquid-vapor interface during evaporation or condensation, as was done in Ref. 26. However, this temperature jump may be neglected unless the phase change occurs too rapidly. 24 In the derivation of our model, we manipulate the mass flux balance equation at the liquidvapor interface. In the literature, the effect of mass loss or gain at the liquid surface due to evaporation or condensation has been included into the model through this equation by assuming that the vapor speed is much larger than the liquid one because the vapor density is much smaller than the liquid one. 1,2 However, in this paper we show another interpretation of this equation as a consequence of order estimate. If we assume that the degree of the disparity in density between two phases is much greater than that in viscosity, which is valid for most substances far from the critical point, we can find that the liquid velocity is not balanced by the vapor one in the mass balance equation. Instead, it is balanced by the interface velocity, which indicates that the effect of mass loss or gain can be neglected. This order estimate leads to the approximation of the mass balance equation, where the total mass of the liquid is conserved and the effect of the mass flux affects only the dynamics of the vapor. Under this approximation, we derive the model where the conservation of the total liquid mass is guaranteed, whereas the effect of evaporation or condensation remains in the vapor dynamics.
The paper is organized as follows. The model is formulated in Sec. II, where the interfacial boundary condition and scaling peculiar to this system mentioned above are introduced.
Linear stability results are presented in Sec. III, including analyses of the Rayleigh-Taylor instability and the effect of degree of nonequilibrium on it. Section IV summarizes the results and future work.

II. FORMULATION
For simplicity we consider a two-dimensional system as in Fig. 1, where the horizontal bilayers, liquid and vapor of the same substance, are confined by rigid parallel walls from both below and above. We assume that the initial equilibrium thicknesses of the liquid and vapor layers, d l and d v , are small enough to ignore the buoyancy effect. The temperatures of the liquid-side and vapor-side plates are controlled at T 1 and T 2 . The x axis is taken to be parallel to the walls, and the z axis perpendicular to them. The z = 0 plane corresponds to the boundary between the liquid and the liquid-side plate. The position of the liquid-vapor interface is described by z = h(x, t). Gravity acts in the negative direction of the z axis.
We assume that the equations of continuity for incompressible fluids and of the momentum and energy balance hold in each phase. They are given by, respectively, Here, v β = (u β , w β ), p β and T β are velocity, pressure and temperature fields, respectively, in the β = {v, l} phase, where v denotes the vapor and l the liquid. The differential operator is ∇ ≡ (∂ x , ∂ z ). The coefficients ρ β , η β and κ β denote the density, dynamic viscosity and thermal diffusivity in the β phase, respectively, which are assumed to be constant in each phase. In Eq. (1b), g is the gravitational acceleration and e z the unit vector in the z direction.

B. Boundary conditions
At the walls (z = 0 and z = d l + d g ), we impose no-slip boundary conditions. Along with the temperature conditions prescribed above, they read At the liquid-vapor interface z = h(x, t), the mass flux J must be conserved: Here, n is the unit normal vector directed toward the vapor, and v I represents the interface velocity, which satisfies the kinematic condition We assume the continuity of the tangential velocity along the interface, where t is the unit tangent vector to the interface, The interfacial stress and energy balance equations read, respectively, 1,2 where E β , σ, H, L and λ β are the rate-of-strain tensor in the β phase, the surface tension, the mean curvature of the interface the latent heat and the thermal conductivity in the β phase, respectively. In Eq. (8), we ignore the thermocapillary (Marangoni) term and assume σ as well as L and λ β to be constant for simplicity. The Marangoni effect on two-phase surfaces has often been neglected in the literature. 5,6,7,18 The recent investigations on the linearized systems of liquid-vapor layers 8,10,11 and drop or bubble 25 , show that the Marangoni effect has little significance in pure two-phase coexisting states, because the liquid-vapor interface becomes almost isothermal owing to the large entropy difference between the two phases. In the Appendix, we include the thermocapillary term in our model and examine its effect. We numerically find that the thermocapillarity makes little contribution to linear stability of a stationary state of the model at least in the physical situations considered here. The projection of Eq. (8) on the normal and tangent to the interface yields, respectively, where Eq. (6) was used in the second equation. Assuming the moderate phase change rate, the continuity of the temperature at the interface holds: Finally, in order to close the system we require an additional boundary condition, which relates to the interfacial thermodynamic state. In this study, we adopt the linearized phenomenological law such that the mass flux across the interface is proportional to deviation from local thermodynamic equilibrium: 24 Here, µ β is the chemical potential in the β phase, which is a function of the pressure in the corresponding phase and the temperature at the interface. A proportionality coefficient K will be later specified by analogy with the kinetic theory. We now expand the chemical potentials into Taylor series in this equation around their initial equilibrium value µ 0 with respect to the variations of the pressure and the temperature, where T sat (p 0 ) is the saturation temperature at the initial equilibrium pressure p 0 . Using the Gibbs-Duhem relation for a one-component system, we obtain for each phase where s β is the entropy density of the β phase. Substituting Eq. (15) into Eq. (13) yields where ∆s ≡ s v −s l is the entropy difference between the phases and related to the latent heat by L = T sat (p 0 )∆s. If we neglect the pressure terms on the right hand side of Eq. (16), we recover the usual kinetic relation. 1,2,6,18,19 On the other hand, in the limitK → ∞, Eq. (16) reduces to the condition of local thermodynamic equilibrium, used in several phase-boundary problems. 8,10,11,25 Therefore, Eq. (16) is an intermediate relation connecting the two different interfacial conditions appearing in the studies of two-phase systems with phase change.

C. Dimensionless equations and parameters
In order to nondimensionalize the above equations, we scale lengths, time, velocities, pressures, temperatures and mass flux by d l , where ∆T is the initial temperature difference across the liquid layer. We find ∆T together with T sat (p 0 ), solving Eq. (1c) for both phases with the boundary conditions (2), (9) and (12) in the equilibrium steady state (J = 0), as follows: Here, the dimensionless parameters λ and d have been introduced. The definitions of dimensionless parameters appearing in this paper are presented in Table I. Furthermore, we define the dimensionless pressure and temperature such that their initial equilibrium values at the interface, p 0 and T sat (p 0 ), correspond to 0 in their new variables. In the following, we show the resulting nondimensionalized equations.
First, the governing equations of the liquid layer (1) become and those of the vapor layer The boundary conditions at the walls (2) reduce to and those at the interface (3), (5), (6), (11), (9), (12) and (16), respectively, In Eq. (21h) we have introduced the dimensionless parameter K, instead ofK in Eq. (16), which has the dimension. The parameters K andK are related through K = d l LK/λ l .
The value of the parameter K defined in Table I is determined by the comparison with the kinetic theory (Hertz-Knudsen law). In its definition, α is the accommodation coefficient, m is the molecular mass of the fluid and k B is the Boltzmann constant.

D. Long-wave asymptotics
We apply the long-wave approximation to both layers. 2 Letting a small parameter ǫ be d l /Λ, where Λ represents the characteristic lateral length scale, new space and time variables are introduced as This rescaling indicates that the physical quantities vary much slower in the horizontal direction than in the vertical one. Assuming that ǫ ≪ 1, we expand the velocities, the pressures and the mass flux in powers of ǫ as follows: Here, we required w β /u β = O(ǫ) based on the continuity equations. To take the pressure In order for all the physical effects to appear in the leading-order equations, some of the dimensionless parameters shown in Table I are scaled by ǫ as where we set Θ → 0 to neglect the molecular kinetic energy and viscous dissipation terms in Eq. (21f). The tildes denote the quantities of order O(1), which will be used instead of the original dimensionless parameters. Here we assumed that the evaporation number  (21a)). Since the density ratio ρ is much smaller than the dynamic viscosity one η for most substances far from the critical point, we set ρ/η = O(ǫ). Moreover, we took the thermal conductivity ratio λ as O(1) in the scalings (24), although the real value of λ is very small because in general the thermal conductivity of gas is much lower than that of liquid. We cannot take the limit λ → 0 because the temperature boundary condition at the vapor-side wall (20b) contains the factor λ −1 . We substitute these scalings into the previous dimensionless equations and take the limit ǫ → 0, so that only the leading-order terms in ǫ are left in the equations. Hereafter, we shall omit the primes, the tildes and the subscripts 0, unless otherwise stated.
The leading-order governing equations are for the liquid layer (0 < z < h) from Eqs. (18), and for the vapor layer (h < z < 1 + d) from Eqs. (19). Whereas the wall boundary conditions (20) remain unchanged, those at the interface (z = h) from Eqs. (21) result in Note that from Eq. (21a) with Eq. (21b) we obtain the two decoupled equations (27a) and (27b), because the liquid and vapor velocities in Eq. (21a) do not balance under the scalings (24). In other words, since we have scaled the density ratio as O(ǫ 2 ), the terms representing mass loss or gain become the next order in ǫ in Eq. (27a), and hence are discarded there.
The next order mass balance gives Eq. (27b), where the effect of the interface velocity in Eq. (21a) has been neglected compared to that of the vapor velocity. Therefore, Eq. (27a) implies that the total mass of the liquid layer is conserved at the leading order in ǫ, while the effect of the mass flux affects only the dynamics of the vapor layer through Eq. (27b).
In addition, the vapor recoil term in the normal stress balance (21d) and the liquid pressure term in the thermodynamic relation (21h) have disappeared with these scalings.
The origin of the decoupling of the mass flux balance equation (21a) is more specifically explained as follows. As is mentioned below Eq. (23), |v l |/|v v | is the same order as the dynamic viscosity ratio η. From the fact that ρ ≪ η, it follows that ρ ≪ |v l |/|v v | or ρ|v v | ≪ |v l |. For the second equality of Eq. (21a) to be true, the liquid velocity v l · n must be balanced by the interface velocity v I · n, which is represented by Eq. (27a) through Eq. (21b). In Eq. (27b), v I · n has disappeared because from Eq. (27a) it is the same order as v l · n, much smaller than v v · n from the small viscosity ratio η ≪ 1. Therefore, in order for the decoupling of Eq. (21a) into Eqs. (27a) and (27b) to be valid it is essential that ρ ≪ η ≪ 1 be the case.
Solving Eqs. (25d) and (26d) with the boundary conditions (20) and (27g) yields the temperature gradients in both layers Then, the surface temperature T I can be explicitly expressed as Substituting this equation into Eq. (27h) again, we find the expression for the mass flux with The expressions for the vertical velocities w l and ηw v immediately follow from the integration of the continuity equations (25a) and (26a) with the no-slip boundary conditions (20): Integration of Eqs. (25c) and (26c) with Eq. (27d) gives the relation between the liquid and vapor pressure where the liquid pressure gradient ∂ x p l has been eliminated using Eq. (35). Equations (31) and (37) can be combined to eliminate J: Equations (36) and (38) whereĥ andp are infinitesimal quantities. Linearizing the system gives the following growth rate: with Here, positive (negative) values of ω indicate instability growth (decay). The parameters A and k 2 0 represent the effect of lateral vapor pressure fluctuation and that of local saturation temperature variation, respectively. From the dispersion relation (40), one can easily find that the growth rate vanishes in the limit k → 0. Nevertheless, if k 2 0 were not present in Eq. (40), the finite growth rate would remain for k → 0. Therefore, the presence of the saturation temperature variation mitigates long-wave growth rates. Specifically, the saturation temperature is increased (decreased) in higher (lower) vapor pressure regions and thereby the rate of evaporation or condensation is reduced. This effect is prominent for long-wave disturbances. Notice that our model does not admit a quasisteady solution of flat moving interface even if k 0 = 0, because the lateral uniformity leads to ∂ t h = 0 from Eq. (36); this is a direct consequence of the neglect of mass loss or gain. We can consider only the long-wave limit k → 0, where ω = 0 if k 0 = 0.

A. Superheated or supercooled state
To quantify the above results, we consider the water/vapor system at 100 • C and 1 atm.
Using the material properties shown in Table II, we plot the growth rates (40) in Fig. 2, where we set d l = 10 −4 m, ∆T = 0.1 • C, g = 9.8 m/s 2 and α = 1. In the experiment, it is the temperature difference between the plates T 1 − T 2 , not across the liquid layer ∆T that can be controlled. However, in the following we fix ∆T so that the vertical temperature gradient in each phase is constant when we vary the value of the initial thickness ratio d.
From Eq. (17a) ∆T = 0.1 • C corresponds to T 1 −T 2 = 2.8 • C for d = 1 in this system. In the short-wave regime (k ≫ 1), the growth rate is negative because of the effect of the surface tension (S in Eq. (40)), whereas it is reduced in the long-wave regime (k ≪ 1) owing to that of the saturation temperature dependence on the vapor pressure (k 2 0 or Π in Eqs. (40) and (41b)), as mentioned above.
Similarly, the liquid flux is obtained as In terms of q v and q l , the model equations (36) and (37) are rewritten as If we regard h as nearly 1, the vapor and liquid fluxes are Thus, it is shown that the stabilizing vapor flow is responsible for the strong dependence of the instability intensity on the initial thickness ratio and also for that on the viscosity ratio.
It is worthwhile noting that the role of the vapor flow presented here is different from that described in the study of full linear stability analysis of a similar bilayer system, where the penetration of the fluid at the walls is allowed. 8 where If the right hand side of Eq. (46) does not take a real positive value, the cutoff wavenumber k c does not exist and hence the dispersion curve never crosses the line ω = 0. This indicates that the system is linearly stable, because from Eq. (40) the surface tension makes the growth rate negative in the short-wave limit k → ∞. Since k 2 c is negative when a − bk 2 0 < 0, the criterion for the linear stability is expressed as This inequality is equivalent to the condition that the growth rate is negative at infinitesimal wavenumber, which can be found by Taylor  The expression for the fastest growing mode k max can be also analytically obtained from Unlike the results of linear stability analysis of the other film equations, a simple relation between k c and k max cannot be established for our model. Equation (49) seems to be too complicated to find any asymptotic form. However, in the special case k 2 0 = b, it reduces to a simple form This special case is realizable for the water/vapor system considered if we set d l = 2.3 ×10 −5 m for d = 1 and K → ∞.

B. Rayleigh-Taylor instability
We also investigate the Rayleigh-Taylor instability of the system, 6,9,10,13,14,15,18 where gravity acts toward the vapor side. The system of interest is the case that the layers are heated from below or cooled from above, so that the stabilizing effect of evaporation or condensation counteracts the destabilizing one of gravity. In this subsection, we consider the balance between these two effects by changing the signs of gravity and the temperature difference.
Before starting the analysis, we show the difference from the model of evaporating or condensing liquid films with infinitely deep vapor layer. 6 In Ref. 6, the authors made the two assumptions: the much larger gas layer depth than the liquid one and the neglect of the latent heat in the temperature boundary condition at the interface (9) (Eqs. (3) or (4) in Ref. 6). To compare their dispersion relation with ours, we abandon these assumptions and show their dispersion relation derived without imposing them. In our notation, it reads with The assumptions which they made are equivalent to d ≫ 1 and K ≪ 1. Hence, by taking these limit we obtain |A BM | = K|E|, which is consistent with the result of Ref. 6. However, for the water/vapor system considered here K = 1.2 × 10 3 α for d l = 10 −4 m. Therefore their assumption of the negligible latent heat (K ≪ 1) is questionable unless the accommodation coefficient α is very small. Their dispersion relation (51) has the same form as ours (40), except for the presence of k 2 0 . However, their definition of A in Eq. (52) is different from ours (41a) because of the different mechanism: local mass loss or gain at the interface is the main stability mechanism in Ref. 6. In contrast, this effect has been neglected in our system compared to that of the vapor pressure fluctuation. This is also confirmed by taking the ratio between both A: . Note that this ratio is independent of K.
For the water/vapor system, |A BM |/|A| ≃ 0.011 for d = 2. Therefore |A BM | ≪ |A| holds in the system considered here. Here the condition for the neglect of the effect of mass loss or gain, ρ ≪ η, is required, as was mentioned before. However, if d is much larger, the effect of mass loss or gain can be comparable to that of vapor pressure fluctuation, because the latter effect is much more weakened for a thicker vapor layer.
The cutoff wavenumber for this case is given by The second line of this equation suggests the existence of two cutoff wavenumbers when −|a| + |b|k 2 0 < 0. Recall that −|a| + |b|k 2 0 < 0 is the condition that the growth rate is negative at infinitesimal wavenumber. Therefore, the instability starts around k = 0 for the first case of Eq. (53), whereas finite modes between the two cutoff wavenumbers are unstable for the second. The finite critical wavenumber where the instability starts for the second case is k crit = |b| − k 2 0 /2.
In seeking the critical condition for the stability, it is desirable to vary the liquid depth d l independently. In the dispersion relation (40), we have three dimensionless parameters depending on d l , k 2 0 , G and S. Since S is proportional to d l , we choose S as a control parameter and the remaining two parameters are scaled by S to obtain new parameters independent of d l : For simplicity, here we ignore the effect of K by assuming the local thermodynamic equilibrium K → ∞ or K 1 + K + λ/d = 1 in Eqs. (41) because K ≫ 1 for α = 1 as mentioned above.
The stability condition for the Rayleigh-Taylor instability is obtained from the cutoff wavenumber (53) similarly to the previous case and expressed in terms of the above new parameters as where S c = (k 2 * 0 /|G * |) 1/4 . Here the two cases of this condition each correspond to those of the cutoff wavenumber (53). The first line represents the condition that k 2 c is negative, which is essentially identical with the previous stability criterion (48), while the second the one that k 2 c is not real. In Fig. 3 we plot the neutral stability curves in S vs. |E| plane with the remaining parameters fixed as shown in Table III

C. Effect of degree of nonequilibrium on the Rayleigh-Taylor instability
In Fig. 3, we did not consider the effect of nonequilibrium because K ≫ 1 for α = 1 in the system considered. However, the accommodation coefficient α can be much less than unity and thereby K might approach one. Here, we examine its effect on the stability of the system. As was mentioned above, the degree of nonequilibrium K does not enter the stability condition (48) for the thermodynamic unstable case. For the Rayleigh-Taylor unstable case, however, the criterion (55) includes K only in the lower line (S > S c ) when K is finite. Then, it becomes ρd 4 η(4 + 3d)(1 + d) Therefore, the coefficient of the left hand side of this inequality increases as K decreases, which may render the system unstable even if it is stable at large K. In Fig. 5 we show the stability diagram in α vs. |E| plane for the water/vapor system. We have discussed the instability of thin liquid-vapor layers bounded by rigid parallel walls from both below and above. In this system, the interfacial instability is induced by lateral vapor pressure fluctuation, which is in turn attributed to the effect of phase change: the vapor pressure becomes higher at an evaporating portion of the interface and vice versa.
The liquid is driven away from the higher pressure place and pulled up to the lower pressure one. This pressure-induced instability mechanism has not been considered in the past.
In the formulation, the interfacial boundary condition taking into account the pressure dependence of the local saturation temperature, Eq. (16), was imposed. The relation (16) is an extension of the kinetic law and the condition of local thermodynamic equilibrium, conventionally used in the literature. We applied the long-wave approximation to both liquid and vapor layers, assuming that the layers are too thin for thermal convection to occur. The choice of scalings we adopted (24) allows us to decouple the mass flux balance equation (3), resulting in the neglect of mass loss or gain through the interface by evaporation or condensation at the leading order. In the dimensional form, we can formally write the decoupled equations (27a) and (27b) as where the first and second equations correspond to the leading-order and next-order equations in ǫ. The condition for this approximation is ρ ≪ η ≪ 1: ρ ≪ η for the neglect of the effect of mass loss or gain and η ≪ 1 for that of the interface velocity in the second equation of Eq. (57). To our knowledge, this decoupling approximation of the mass flux balance between two phases has never been considered and might have a possibility to be applied to other phase-change problems where ρ ≪ η ≪ 1 holds.
As a result, a set of the equations describing the temporal evolution of the interface of the liquid-vapor layers has been derived. One of the equations (36) is written in the conserved form for the film thickness, as is expected from the decoupling of the mass flux balance mentioned above. Therefore, the total mass of the liquid layer is conserved in our model.
On the other hand, the effect of the lateral vapor pressure gradient induced by phase change is included in Eq. (36) and the vapor pressure is enslaved to the film thickness through Eq. (38). The result of the linear stability analysis of this model shows that the presence of local saturation temperature variation by the vapor pressure mitigates the growth of longwave disturbances. The instability is enhanced for the smaller initial thickness ratio d and larger dynamic viscosity ratio η, which increase the kinetic resistance to the vapor flow. The role of the vapor flow is to mitigate the effect of the lateral vapor pressure variation, which is different from that described in Ref. 8. We also determined the criterion for the linear stability of the system and found that only slight temperature gradients are sufficient to overcome the stabilizing gravitational effect for the water/vapor system.
We also considered the Rayleigh-Taylor instability of the system. Here, the stabilizing vapor pressure effect is balanced with the destabilizing gravitational effect under very small temperature difference between the plates. Again, the thinner vapor layer strengthens the stabilizing effect of lateral vapor pressure fluctuation. However, for this case the instability domain may be widened if the accommodation coefficient α is below a certain critical value.
This value is about 10 −2 for the water/vapor system of d l = 10 −4 m.
In this paper, we addressed only the linear stability of our model. The next step is to proceed to the nonlinear analysis of the model. We shall investigate the behavior of the solution of the equations in the nonlinear regime by means of numerical simulation. Three dimensional computation of the model will reveal the possibility of occurrence of pattern formation as reported in Ref. 6.
The validity of the long-wave approximation, on which our model is based, should be examined. To this aim, the comparison with the full numerical simulation is required.
In particular, we do not know the critical condition for the onset of small-scale cellular convection in liquid-vapor layers, which is neglected within the framework of the long-wave approximation. There should exist a critical temperature gradient or thicknesses of the layers (critical Rayleigh number) for the transition between the conductive and convective states of the temperature fields if the buoyancy effect is taken into consideration. It would be also of interest to make a comparison with the existing full linear stability analyses 8,9 and even weakly nonlinear analysis 27 of bilayer systems. The effect of heat convection should be estimated to clarify the role of the vapor flow, which is different between theirs and ours.
Finally, we make two remarks on our model. First, we set the scalings in ǫ on the dimensionless material parameters as Eq. (24) assuming their values for the water/vapor system at 100 • C and 1 atm as a representative substance. Hence, if the material properties are considerably changed (e.g. near the critical point), we must reset the scalings appropriate for the relevant values, which may lead to different evolution equations. Second, we can incorporate the effect of intermolecular forces as disjoining pressure in the formulation.
This effect will be dominant for layers of thicknesses below 100nm, as in Ref. 16 and 17. However, for this scale the application of continuum theory to the vapor layer would not be valid, because the mean free path of a gas molecule amounts to about 60 or 70nm at atmospheric pressure and becomes much larger at reduced pressure.

Acknowledgments
The author thanks Alexander Oron for numerous suggestions which improve the paper.
The author's visit to his laboratory in Technion-Israel Institute of Technology was financially supported by the Bilateral International Exchange Program (BIEP) of the Global COE "The Next Generation of Physics, Spun from Universality and Emergence" from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. The author also is grateful to Sadayoshi Toh, Tomoaki Kunugi and the participants of the 4th International Marangoni Association Conference (IMA4) held at Noda, Japan in 2008, for helpful discussions.

APPENDIX: EFFECT OF THERMOCAPILLARITY
In the main text, we ignore the thermocapillarity for simplicity. Here we examine its effect on our model. If the thermocapillarity is present, we must add the thermocapillary term in the stress balance at the interface (8): where ∇ s = t · ∇. Accordingly, Eq. (11b) is modified as Here we use the linear approximation for the dependence of the surface tension on the temperature, where σ 0 is the surface tension at the initial equilibrium temperature T sat (p 0 ), corresponding to σ in the main text. The coefficient γ = − dσ dT I is positive for most common substances.
Then, the nondimensionalized tangential stress balance equation (21e) becomes where the Marangoni number is defined by In order to include the thermocapillary effect in the model, the Marangoni number is scaled as under the long-wave approximation. Then, Eq. (27e) reduces to where the tilde over M is omitted. We note here that the temperature dependence of the surface tension (60) also allows the tangential variation of the normal capillary stress in Eq. (58). However, this effect is safely neglected in the framework of the long-wave approximation (see Ref. 1). Using Eq. (64), one of the coefficients in Eq. (33) is modified as while the other does not change. As a result, only one of the equations (36) contains the Marangoni term: Here, the interfacial temperature gradient is calculated from Eq. (30) as ∂ x T I = 1 to d −1 l (independent of ∆T ) and the latter to d −2 l and ∆T . If K decreases, the former decreases and the latter increases. However, the latter does not diverge as does the left hand side of Eq. (56). In both thermodynamic and Rayleigh-Taylor unstable cases, where the vapor pressure and gravity effects counteract each other, the Marangoni effect acts as amplifying the former and diminishing the latter. Finally, we note that we do not know whether the Marangoni effect is negligible in the nonlinear regime, which will be ascertained in the ongoing numerical analysis.