Thermodynamic Bifurcations of Boiling in Solid-State Nanopores

Boiling heat transfer is the basis of many commonly used cooling techniques. In cooling of electronic devices, for example, it is desirable to further miniaturize heat exchangers to achieve higher heat transfer, and thus it is necessary to understand boiling phenomena on shorter spatial and temporal scales. This is especially challenging at the nanometer scale because conventional imaging techniques cannot capture the dynamics of nanobubbles, owing to the Abbe diffraction limit. Here in this research, we utilize the nanopore Joule heating system that enables the generation of nanobubbles and simultaneous diagnosis of their nanosecond resolution dynamics using resistive pulse sensing. When a bias voltage is applied across a silicon nitride nanopore immersed in an aqueous salt solution, Joule heat is generated owing to the flow of ionic current. With increasing voltage, the Joule heating intensifies, and the temperature and entropy production in the pore increase. Our sensing results show that nanopore boiling follows the theory of minimum entropy production and attempts to settle to a minimum dissipative state. This results in two boiling bifurcations corresponding to the transition between different boiling states. These characteristics of nanopore boiling are represented by an"M"-shaped boiling curve, experimentally obtained from the Joule heat variation with the applied voltage. A theoretical framework is proposed to model the thermodynamics of nanopore bubbles and estimate the system dissipation which explains the four arms of the"M"-shaped boiling curve. The present study reveals that the utilization of nanopore boiling as a benchmark platform offers a valuable means for investigating the intricate boiling phenomenon and its correlation with nanoscale bubble dynamics. This would provide fundamental insights into the chaotic transition boiling regime, which is least understood.


I. INTRODUCTION
Boiling is an unsteady nonequilibrium process that plays a vital role in diverse engineering applications, including thermal desalination [1], power generation [2], cooling of electronic devices [3][4][5][6], inkjet printing [7], and spray quenching [8].Much attention has been paid to extending the nucleate boiling regime during which a huge amount of heat flux can be dissipated from a solid surface through microlayer evaporation [9][10][11][12].However, there is a lack of comprehensive understanding of boiling at the nano-to microscale, where much of the process takes place below the optical limit (Abbe diffraction limit) [13,14].
In recent years, both pool boiling experiments [15] and nanoparticle boiling [16][17][18] studies have shown that vapor nanobubbles of dimensions below the roughness scale can exist stably on a hydrophilic boiling surface.Several plausible and self-consistent theoretical explanations, such as contact-line pinning [19], Knudsen effects [20], and dynamic equilibrium [21], have been proposed to explain the stability of nanobubbles created by dissolved gas.However, their applicability to vapor bubble stability in the context of boiling heat transfer is still not known.Vapor bubble dynamics, which dominate boiling heat transfer, differ significantly from gas bubble dynamics, particularly in the mechanisms of interfacial heat and mass transport and the associated time scales [22].Hence, to handle the complexity associated with the nonequilibrium thermodynamics and nonlinear dynamics [23] of nanoscale boiling, we need to start from a thermodynamic stability analysis of vapor nanobubbles.In the case of macroscale pool boiling, Shoji and co-workers [24][25][26] applied Prigogine's theory of self organization [27] to explain the "S"-shaped boiling curve and clarified the entropic origins of the two bifurcation points: the onset of nucleate boiling and the onset of stable film boiling.In a similar vein, we study the thermodynamic bifurcations of nanoscale boiling, based on well-controlled bubble dynamics experiments using a nanopore Joule heating platform [4,5,14,28,29], which allows us to capture vapor nanobubble dynamics at nanosecond resolutions.
Using this platform, we have previously characterized nucleate, transition, and film boiling regimes based on acoustic and resistive pulse sensing signals [14].In the present study, we construct an "M"-shaped boiling curve [Fig.1(e)] summarizing the boiling regimes for a 199 nm diameter solid-state nanopore (Fig. S1 in the Supplemental Material [30]), explaining the two bifurcation points based on the variation of system dissipation.Furthermore, we shed light on the thermodynamics of the film bubble interface, elucidating the role of (i) Knudsen bubble stability and (ii) its out-of-equilibrium behavior on the system dissipation and boiling curve for the nanopore system.
The nanopore was fabricated using focused ion beam etching on a 100 nm thick suspended silicon nitride membrane chip [Fig.1(a)].When bias voltages are applied across the nanopore after its immersion in 3M NaCl solution, the flow of ionic current I results in focused Joule heating in the nanopore liquid [Fig.1(a)], leading to bubble nucleation [5,28,29].The formation of bubbles within the nanopore hinders ionic current flow, which is captured in a high-bandwidth oscilloscope [Figs.1(b) and 1(c)].For detailed information about experimental methods, we refer to our earlier works [5,14].Discrete and periodic current blockage signals represent the nucleation of a homogeneous vapor bubble at the pore center, followed by its thermal growth and ejection, and finally its shrinkage at the pore entrance [Fig.1(b)] [14], following which there is a waiting/heating period before the subsequent bubble nucleation.We term this sequential bubble generation as nanopore nucleate boiling (NB).As the bias voltage is intensified, the nanopore surface becomes blanketed by a torus-shaped vapor film (film boiling, FB), causing a continuous decrease in baseline current [14].In Fig. 1(c), we find that when the voltage is increased to 6.83 V, a step-like decrease in baseline current (corresponding to FB) occurs multiple times during periodic current dips (corresponding to NB).This resem-bles intermittent film boiling, which falls under transition boiling.Following the theory of minimum entropy production, the nanopore system attempts to settle to a minimum dissipative state (superheating, NB, or FB) through fluctuation-driven activation processes such as bubble nucleation and vapor film formation.The tristable dissipation well diagram in Fig. 1(d) shows increasing dissipation wells of superheating and decreasing dissipation wells of NB and FB with bias voltage.This explains the gradual system transition to nucleate and then intermittent film boiling as observed in experiments.
The variation of Joule heat dissipation with bias voltage is summarized in the form of an "M"-shaped boiling curve [Fig.1(e)].In the following sections, after an explanation of the two bifurcations [ONB and OFB in Fig. 1(e)] in the "M"-shaped boiling curve, the superheating and stable film boiling regimes are explained through steady-state continuum simulations of Joule heating, using kinetic theory boundary conditions for interfacial heat and mass transfer (Fig. 2).Furthermore, nucleate boiling and transition boiling, which constitute the unsteady boiling regimes, are described.The dissipation function, which represents the entropy generation of the nanopore system, has been theoretically calculated for all the boiling regimes and shown to be "M"-shaped with respect to bias voltage.

A. Boiling bifurcation diagram
In the present study, we apply a ramp voltage pulse [linearly increasing from 4.5 V to 9 V across 30 ms, V app = 4.5 + 150t [s]; Fig. S4(a) in the Supplemental Material [30]] across the nanopore, which was initially in thermal equilibrium at an ambient temperature of 298.15 K.The variation of mean Joule heat generation inside the nanopore during the 30 ms and the boiling modes at discrete time points are studied to explain the global and local aspects of nanoscale boiling.Owing to the very slowly changing bias voltage (nonequilibrium boundary condition), each boiling mode reached during the ramp voltage pulse can be assumed to be in a quasi-steady state.For reference, the thermal relaxation time for the 199 nm nanopore [τ th ∼ R 2 p /D th = O(100 ns)] and the lifetime of nucleate bubbles [t b = 67 ns, Fig. 1(b)] is much shorter than the millisecond-order voltage pulse.Here, R p is the pore radius and D th is the thermal diffusivity of water.
Similar to pool boiling, nanopore boiling is a nonequilibrium phase change process, but the critical difference is the boundary condition that causes nonequilibrium [31].In pool boiling experiments [9], the externally controlled wall temperature is the boundary condition that makes the system nonequilibrium [26], whereas in nanopore boiling, the bias voltage acts as the nonequilibrium boundary condition.Accordingly, we plot the nanopore boiling curve with the applied bias voltage V app on the horizontal axis and the resulting Joule heat dissipation V app × I on the vertical axis [Fig.1(e)].It should be noted here that the creation of Joule heat in a conductor is a significant source of entropy production [32,33] and accordingly an appropriate quantitative measure of the nonequilibrium nature of the nanopore system.Being an open dissipative system, nanopore boiling tries to self-organize and minimize entropy production from Joule heat dissipation.This leads to the two bifurcation points: onset of homogeneous nucleate boiling and onset of stable nanotorus film boiling [Fig.1(e)].At both these points, the slope of Joule heat versus voltage changes from positive to negative, i.e., a decreasing trend of dissipation is achieved, which is a more favorable thermodynamic pathway for the system to progress.

Superheating
As the voltage is increased from 4.5 V to 6.65 V, the nanopore current increases continuously [Fig.S4(a) in the Supplemental Material [30]].As a result, the generation of Joule heat within the nanopore increases monotonically and constitutes the left arm of the "M"-shaped boiling curve.We use continuum simulations to implement a Joule heating model, which has been described in our earlier papers [5,14] and is restated in Sec.S2 in the Supplemental Material [30].Using the Joule heating model, the steady-state nanopore temperatures at the pore center (T c ) and pore walls (T w ) during superheating [Fig.S4(c) in the Supplemental Material [30]] are obtained.The contour plots in Figs.3(a) and 3(b) show the spatial distribution of nanopore temperature for bias voltages of 6.0 V and 6.65 V, respectively.The theoretically obtained nanopore current and Joule heat generation [yellow lines with markers in Figs.S4(a) and 1(e), respectively] fit well with the experimental results.According to the simulations, we find that T c and T w exceed the saturation temperature of 373.15K at voltages lower than 5.5 V.However, no bubble signals are seen until 6.7125 V, when T c = 562 K is reached.This value is close to the kinetic limit for homogeneous nucleation, 575 K [34].The enhanced range of metastability of the nanopores and the prevention of heterogeneous nucleation at low voltages were shown by Golovchenko and co-workers [28,29], and were later explained by Paul et al. [5] through a nucleation theory model based on ripening competition between homogeneous and heterogeneous bubble clusters.This unique characteristic is due to a large temperature difference between T c and T w , which originates from the focused nature of Joule heating.
We also calculate the total system dissipation during superheating, Ψ = V ψ dV , where the local dissipation in the liquid volume can be written as ψ = T d i s/dt = T ṡgen [33].Here, T is the liquid temperature and d i s/dt is the entropy generation per unit volume, which can be expressed as Here, J u = k w ∇T is the heat flux in the liquid and J = σE is the ionic current, where k w is the thermal conductivity and σ is the electrical conductivity, and E = −∇ϕ is the electric field generated in the nanopore liquid due to the applied bias voltage ∆ϕ = V app .The first and second terms on the right-hand side of Eq. ( 1) represent the entropy generation due to heat conduction and Joule heating, respectively.Here, we neglect the viscous dissipation term τ : (∇u ) is assumed, which is consistent with the nanopore electroosmotic velocity [35,36].Using Eq. ( 1), we calculate the system dissipation Ψ, which increases monotonically with voltage during superheating [Fig.3(g)].The breakdown of dissipation sources is shown in Fig. S4(c) in the Supplemental Material [30], where Joule heating in the liquid is the dominant source of dissipation, while heat conduction in the silicon nitride membrane makes the least contribution.The superheating regime is followed by periodic nucleate boiling and intermittent film boiling, during which the time-averaged dissipation decreases with voltage (discussed further in the following sections), causing the formation of a stable dissipative structure.

Stable film boiling
Following the second bifurcation at 8 V, stable film boiling (SFB) commences.During this regime, a torusshaped vapor film blankets the cylindrical pore surface and exists stably despite the steep temperature fields and complex heat flow within the nanopore.As the bubble occludes the pore volume, there are catastrophic decreases in the nanopore current [Fig.S4(a) in the Supplemental Material [30]] and Joule heat generation [Fig.1(e)].In our previous work [14], we found that the magnitude of the current drop was proportional to the volume of the pinned torus bubble.This relationship was found to be valid for 340 nm, 420 nm, and 460 nm diameter pores.These trends strongly suggested the presence of a pinned torus bubble on the pore surface.In Fig. S4(a) in the Supplemental Material [30], we find that as the voltage is increased, the nanopore current decreases, which can be explained by the bulging out of the vapor film pinned at the nanopore wall.In other words, the contact angle of the bubble [θ in Fig. 2(a)] increases with voltage.As the bubble occupies a larger fraction of the pore volume, both experimental [purple trace in Fig. 1(e)] and simulation [black line in Fig. 1(e)] studies of vapor film bulging out show that Joule heat dissipation is reduced.We study the exceptional stability of the torus bubble through a continuum model of heat and mass transfer and explain the bulging-out effect of the vapor film with increasing bias voltage.Details of the simulation results for stable film bubbles are given in Sec.S3 in the Supplemental Material [30].
Nonuniform and concentrated Joule heat generation [Fig.S8(c) in the Supplemental Material [30]] causes nonuniform thermal gradients on the torus bubble surface.The shape of the torus bubble is determined under the assumption that it is pinned to the surface of a cylindrical pore of thickness L = 100 nm with uniform mean curvature K [Fig.2(a)].A further description can be found in our previous paper [14].Under these constraints, the contact angle θ is the only free geometric parameter that controls the bubble shape.The larger the angle θ, the more the bubble expands [Figs.3(d)-3(f)], blocking much of the pore volume and limiting current flow and Joule heat generation.Now, the tip-to-base height of the bubble satisfies h max ≈ 3λ, where λ is the thermal mean free path of vapor molecules, which satisfies λ = 1.922λ . λ HS is the mean free path per hard sphere model.Here, m is the mass of one vapor molecule, ρ v is the vapor density, and d is the collision diameter, assumed to be 2.7 Å.For the Knudsen bubble, the volume can be divided into liquid-vapor (LV) and solid-vapor (SV) kinetic interfaces and the vapor bulk regions.Each interface has a height of one mean free path (h LV ∼ h SV ∼ λ), and the vapor bulk region has a length of less than one mean free path [h bulk < λ in Fig. 2(b)].Figure 3(h) shows the height variations of the equilibrium film bubble at varying voltages, where we find the average bubble height of the SFB to be h avg ≈ 1.16λ.Hence, we consider temperature drops across the LV and SV interfaces (net thickness 2λ) to account for the net temperature drop across the bubble.
Considering the heat transfer at the SV interface to be governed by molecular reflections [Fig.2(c)], we express the heat flux as where R g is the universal gas constant, and k m is the thermal conductivity of the silicon nitride membrane, taken as 3.2 W m −1 K −1 [28].q bulk is the incoming heat flux from the vapor bulk region to the SV interface.ρ v , c v , and U v are the vapor density, specific heat at constant volume, and internal energy obtained from the phase diagram of water according to the IAWPS formulation [38] [Fig.S7].As shown in Fig. 2  the LV interface can be written as Here, the subscript "sat" denotes the vapor phase evaporating at the liquid side of the LV interface, which is assumed to be saturated at T LV .The latent heat of evaporation or condensation (h fg ṁ) associated with the net emission or trapping of molecules at the liquid side of the LV interface [39] consumes a significant amount of the incoming heat flux from the bulk liquid [q in = k w n • ∇T in Figs.2(d) and 3(f)].The remaining heat (q int = q evap − q cond + q ref ) is transferred through the molecular fluxes of evaporation, condensation, and reflection at the LV interface [J evap , J cond , and J r,i as shown in Fig. 2(d)].The molecular velocities in these fluxes are defined based on Maxwellian distributions [40][41][42].Accordingly, . The net evaporation/condensation mass flux across the LV interface can then be written as where α e and α c are the evaporation and condensation accommodation coefficients, respectively.For the results shown in Fig. 3, α e = 0.9 is uniformly applied on the LV interface.Applying the local balance of fluxes on the LV interface (⟨J out ⟩ − ⟨J coll ⟩ = ⟨J evap ⟩ − ⟨J cond ⟩), α c can be written as We also apply the Laplace equation on the LV surface along with closure conditions of no heat and mass accumulation to solve for the steady-state bubble size and temperature: Here, P v and P w are the vapor and liquid pressures, respectively.S LV and S SV are the liquid-vapor and solidvapor interfaces respectively as shown in Fig. S2 (inset) in SI [30].γ is the surface tension of water defined at T v .We solve the governing equations for Joule heating of nanopores described in Sec.S2 in the Supplemental Material [30] under the above-described heat and mass flux boundary conditions on a finite volume mesh [Figs.S2] using the multiphysics software, arb [43].This yields the Figure 3(f) shows that the bulging-out effect results in locally high Joule heating (high ṡgen ) near the bubble tip and locally high heat inflow (high q in ).Figures S5(a A stable nanobubble with the simultaneous influx and outflux of dissolved gas on the bubble interface was termed dynamic equilibrium or controlled nonequilib-rium by Brenner and Lohse [21].In that paper, the authors constructed a model based on volume conservation and local gas over-saturation near the solid surface to explain the unusual stability of gas nanobubbles from dissolution.Later, Liu et al. [44] performed continuum simulations, showing that simultaneous influx and outflux of dissolved hydrogen on the gas-liquid interface can also create a dynamic equilibrium, explaining both the stability of electrochemically generated nanobubbles and the corresponding steady-state current drop on the nanoelectrode.On a similar note, we show here that heat and mass fluxes can also reorganize on the liquid-vapor interface of a bubble, lending it stability.We define the term, dynamic thermal equilibrium for this mode of thermal nanobubble stability. It should be noted here that the solutions for the dynamic equilibrium size θ eqm for different values of α e at the same bias voltage are very close [blue and turquoise traces in Fig. 3(h)].However, as shown in Figs.S5(f) and S9(f) in the Supplemental Material [30], small values of α e cause negative and unrealistic values of α c in order to satisfy the heat flow from the liquid.Hence, only for high values of evaporation coefficients a stable LV interface can exist.Furthermore, compared to α e = 0.5, when α e = 0.9, we find that α c is less non-uniform and varies over a shorter range of 0.82-0.98 on the bubble surface (Fig. S5f and h).
Here, we should note that a similar dynamic thermal equilibrium can be caused by the self-organization of other interfacial parameters.For example, α c can remain uniform on the LV interface, while α e self-organizes to conduct the non-uniform heat fluxes coming from the liquid.To test this hypothesis, we solved for the bubble size under the condition of α c = 0.9 and obtained a spatially varying α e ranging from 0.51 to 0.94 (Fig. S6c).When uniform α c = 0.9 is applied on the LV interface, the equilibrium size of θ eqm = 117°is obtained at 8.04 V (Fig. S6a).On the other hand, when uniform α e = 0.9 is applied on the LV interface, equilibrium size of θ eqm = 117°is obtained at 8.02 V (Fig. 2d).Also, the vapor temperature is 415.91 K when α e = 0.9 (Fig. S5b) and 416.33 K (Fig. S6b) when α c = 0.9.This indicates that we converge to a similar dynamic thermal equilibrium bubble size and temperature from both assumptions.It should be noted that apart from accommodation coefficient variations, variations in vapor molecule velocity distributions on the bubble surface can also accommodate the non-uniform heat fluxes creating a dynamic thermal equilibrium.Our hypothesis of self-organization of interfacial parameters creating a dynamic equilibrium suggests the existence of dissipative structures [27,45].The spontaneous evolution of far-from-equilibrium systems to organized states is commonly observed in biology and chemistry.In the case of nanoscale bubbles, the exact molecular mechanism may be very complex, and large-scale nonequilibrium molecular dynamics (NEMD) simulations [13] might prove useful in resolving the whole picture.
Despite the large heat and mass transport at the bubble interface, dissipation at the bubble interface is negligible compared with the dissipation due to Joule heating and heat conduction in the liquid [Fig.S13 in the Supplemental Material [30]].Thus, although the interfacial flux intensifies as the voltage increases to maintain dynamic thermal equilibrium, the increase in dissipation at the bubble interface is much smaller than the decrease in dissipation due to the suppression of Joule heat generation by the expansion of the bubble.Hence, as the voltage is increased, the increasing nonuniform temperature distribution at the bubble interface is accommodated through a dynamic equilibrium mechanism, which locally stabilizes the bubble interface.At the same time, the expansion of the bubble reduces dissipation throughout the system, resulting in thermodynamically stable film boiling.

Nucleate boiling
When the linearly increasing voltage reaches 6.7125 V at 14.75 ms, the temperature at the pore center reaches 562 K, and the first homogeneous bubble blockage signal (67 ns in Fig. S3 in the Supplemental Material [30]) is seen, which marks the onset of nucleate boiling.During this regime, the homogeneous bubble nucleating at the pore center undergoes inertio-thermal bubble growth [5,22].We simulated the spherical bubble growth-collapse cycle through a one-dimensional moving boundary model [4,46] discussed in our earlier works.A spherical bubble seed is inserted at the pore center, and the hotspot steady-state temperature distribution along the pore axis at 6.75 V is used as the initial radial liquid temperature distribution around the bubble surface.The vapor temperature and density inside the bubble are assumed to follow the saturation line, and the bubble center is fixed at the pore center.Although only radial motion is considered, good agreement is found in the estimated bubble lifetime of 98 ns compared with the experimental blockage duration of 67 ns.Approximation of radial growth dynamics in the inertio-thermal regime [47] has proven quite useful in simplifying cavitation bubble dynamics [48] and plasmonic bubble dynamics [49].According to our model [4], the bubble radius R reaches a maximum value of 677 nm before shrinking back to the liquid phase [Fig.4(a)].As the bubble radius was larger than the pore radius for almost the entire lifetime, we assume total cutoff of Joule heat generation inside the nanopore.This causes a drastic reduction in dissipation post-nucleation [Figs.4(b) and 4(c)].Furthermore, the growth cycle causes the pre-nucleation hotspot [Fig.4(a)] temperature distribution in the liquid to even-out, leading to a gradual reduction in dissipation during the bubble lifetime [Fig.4(c)].As the Joule heat generation is disrupted during the bubble lifetime, the sharp temperature distribution in the silicon nitride membrane also relaxes, which was simulated using a onedimensional conduction model assuming zero heat flux on the pore surface [5].This rough assumption allows us to simplify the complicated conjugate heat transfer problem during nucleate bubble growth and collapse.However, as the dissipation in the silicon nitride membrane has an overall minor contribution to the total system dissipation [Fig.S4c], said assumption should not have a major impact on the period-averaged system dissipation values during nucleate boiling regime (blue marker in Fig. 2g).
After the bubble growth-collapse cycle, the nanopore Joule heating re-commences, and this is simulated using a transient Joule heating model [5] implemented on an axisymmetric finite volume mesh [Fig.S2].We find that the system dissipation increases [Figs.4(b) and 4(c)] as Joule heating is switched on and the thermal hotspot at the pore center is regenerated after a waiting period [Fig.1(b)], leading to subsequent bubble nucleation.From experiments, we found that the blockage duration remained almost the same, while the waiting times between bubble nucleations decreased steadily as the bias voltage was increased [Fig.4(d)].Accordingly, the periodaveraged system dissipation during the nucleate boiling cycle also decreased with increasing bias voltage (blue trace in Fig. S12 in the Supplemental Material [30]).

Transition boiling
As the pore surface temperatures increase with increasing voltage, nucleation of heterogeneous bubbles becomes likely, and these can coalesce to form a thin vapor film blanketing the pore surface, leading to a decrease in the baseline current.At 6.83 V, we find intermittent film bubbles [the FB zone in Fig. 1(c)] separating periodic bubble and nucleate bubble blockage signals, which signifies transition boiling.We simulate the film bubble shape and temperature using the Joule heating model as described in Sec.S2 in the Supplemental Material [30].However, as the film bubble is thinner than the bulged-out bubble in stable film boiling (SFB), we neglect the SV interfacial temperature drop and apply the boundary condition T SV = T v instead of Eq. ( 2).Accordingly, a steadystate current of 31 µA is obtained for the dynamic equilibrium bubble size [Fig.3(c)], which slightly underestimates the experimental nanopore FB current at 6.83 V [Fig.1(c)].As the average bubble height satisfies h avg > λ during unstable film boiling (UFB) [Fig.3(h)], we estimate the temperature drop in the vapor bulk, ∆T bulk , using a one-dimensional thermomass model that considers the diffusio-ballistic (non-Fourier) heat transport [50][51][52][53][54]. Details of the simulation results for unstable film bubbles are given in Sec.S4 in the Supplemental Material [30].However, when the vapor film has h ave = 1λ, ∆T bulk increases as the bubble expands with increasing voltage [Fig.S9(e) in the Supplemental Material [30]], revealing the limitations of this model.In other words, if the thin vapor film stays in the range h ave < λ, it will not reach mechanical equilibrium and is prone to collapse.This was actually observed in the experiments [Fig.1(c)], where the FB disappears after a few tens of microseconds, allowing nucleate boiling to resume.We estimated the average system dissipation during intermittent film boiling.Details are given in Sec.S5 of the Supplemental Material [30].The average system dissipation (green diamond markers in Fig. S12) was obtained by a weighted average of the nucleate boiling dissipation [Figs.4(b) and 4(c)] and the unstable film bubble dissipation [Figs.S13(c) and S13(d)], based on their relative probabilities of occurrence P NB and P FB as read from the experimental current signals.
As can be seen in Fig. 3(g), following the first bifurcation, the dissipation starts to decrease with the onset of nucleate boiling, and further decreases during intermittent film boiling (IFB), which leads overall to a greater reduction in dissipation, albeit one that is not fully stable [Fig.1(d)].
In the boiling curve [Fig.1(e)], we find that the Joule heat dissipation rises from 7.1 V to 8 V before the second bifurcation event.The mean current I m [the purple trace in Fig. 4(e)] flattens in this voltage range, indicating minor bubble expansion.On the other hand, the current spectrogram [Fig.4(f)] reveals the appearance of a distinct frequency band that rises from 5 MHz to 13 MHz in this voltage range.However, as shown in Fig. 4(e), in the second ramp voltage pulse (P2), unlike the first ramp voltage pulse (P1), the mean current I m decreases monotonically with voltage up to 8.1 V.The spectrogram also shows no frequency band up to 8.1 V at P2 (Fig. S11 in the Supplemental Material [30]).Therefore, we can infer that the mean bubble position is out of equilibrium (less bulged position), where the bubble is in pinned volumetric self-oscillation [14,[55][56][57].In comparison, during nucleate boiling (NB) there is a power spectrum distributed over a wide frequency range, originating from non-sinusoidal current fluctuations comprising unequal blockage durations and waiting times.Accordingly, we can discern the nucleate boiling regime comprising of transient bubble nucleations from the frequency signature as well.
Although estimating entropy generation during UFB is quite cumbersome [58], if we assume that the mean bubble position is determined from dynamic thermal equilibrium while holding I m constant, the system dissipation increases with bias voltage [the black line in Fig. 4(g)].This rise in dissipation and Joule heating due to incomplete film bubble expansion constitutes the second rising arm of the "M"-shaped boiling curve [Figs.1(e) and 3(g)].
On the other hand, assuming that the mean bubble position is determined from mechanical equilibrium, the system dissipation decreases with increasing bias voltage [the pink line in Fig. 4(g)].
For the P2 pulse, the current rises sharply at 8.1 V (jump 1), from which a small frequency band appears and continues until 8.3 V (Fig. S11), and the current drops ubbles in a thermal steady state.Darker shades of blue and green correspond to larger film bubbles.The pink lines shows the variations of current and dissipation when the bubble size also satisfies mechanical equilibrium.The black line indicates the rise in dissipation when the mean bubble is out of mechanical equilibrium and only in a thermal steady state, albeit satisfying constant mean current.
sharply, marking the second bifurcation event.The mean current agrees well with the stable film boiling (SFB) simulation results [the red trace in Fig. 4(e)], and the oscillation spectrum disappears for both P1 [Fig.4(f)] and P2 (Fig. S11).In the UFB region, the bubble swelling increases ∆T bulk , and the heat transport becomes less ballistic, destabilizing the bubble [20].However, it is unclear why the bubble sometimes starts oscillating at a mean position that is not in equilibrium.Furthermore, for the P2 pulse, a shift was observed between the extrap-olated UFB mean current and the SFB current (jump 2).This suggests that a membrane or film bubble with h avg in the range from λ to 2λ (at UFB) follows a different heat transport and stabilization mechanism than one with h avg = 2λ (at SFB), causing this quantized behavior.In summary, the second bifurcation follows a stochastic jump process [59] which can be potentially useful as a switch in nanofluidic computing applications [60][61][62].However, further studies using large-scale nonequilibrium molecular dynamics (NEMD) simulations [13] may help explain the anomalous behavior of the Knudsen membrane bubble.

III. CONCLUSIONS
Nanopore Joule heating serves as an excellent platform to detect single-bubble dynamics at nanosecond resolutions.Based on an "M"-shaped boiling curve for a 199 nm nanopore, superheating, homogeneous nucleate boiling, transition boiling, and stable film boiling regimes are classified.Two bifurcation points, namely, the onset of nucleate boiling and the onset of stable film boiling, are observed, at which the system finds a dissipation reduction mechanism and self-organizes into a more stable boiling structure.During nucleate boiling, the periodicity of bubble nucleation increases gradually with increasing bias voltage, leading to dissipation reduction.At the second bifurcation point, the self-oscillating and entropygenerating film bubble expands catastrophically into a bigger film bubble, leading to a stepwise reduction in system dissipation.These bifurcations are explained by theoretical calculations of the variation in dissipation function using a continuum model.The model yields thermal properties of the film bubble as well as nanopore currents that are in good agreement with experimental results.We show that in addition to thermodynamic effects such as dynamic thermal equilibrium and ballistic heat transfer, confinement effects like contact-line pinning must be considered to explain stable film boiling in experiments.Furthermore, our model suggests that a self-organization of thermal accommodation coefficients at the liquid-vapor interface can account for the transport of high heat and mass fluxes through nanoscale bubbles, creating a dynamic thermal equilibrium, which causes nanobubble stability and unique boiling structures in extreme thermal environments.However, a detailed understanding of the molecular dynamics within such Knudsen bubbles is still a challenge, for which further investigation would be required.

FIG. 1 .
FIG. 1. Boiling bifurcations in a 199 nm diameter nanopore.(a) Schematic of nanopore Joule heating setup.Boiling signals are captured simultaneously using a passive (more noise) and an active (less noise, owing to higher capacitance) oscilloscope probe.(b) Homogeneous nucleate boiling comprises discrete bubble events, resulting in periodic bubble blockage signals.The plot shows filtered current from the active probe.Tc and Tw denote the pore center and pore surface temperatures, respectively.(c) Transition boiling shows intermittent film bubbles disrupting nucleate boiling.The plot shows filtered current from the passive probe.(d) A phenomenological plot of the system dissipation function versus system coordinates (superheating, nucleate, and film boiling) explains the first bifurcation point.(e) An "M"-shaped nanopore boiling curve showing the distinct zones of subcooled liquid, metastable liquid, nucleate boiling, transition boiling, and stable film boiling.The two bifurcation points at the onset of nucleate boiling (ONB) and the onset of stable film boiling (OFB) are indicated.

FIG. 2 .
FIG. 2. Schematic of the model of heat transfer in the Knudsen vapor film bubble.(a) Nanopore torus bubble consisting of the liquid-vapor interface, the solid-vapor interface, and the bulk-vapor region between these two interfaces.(b) Ballistic heat flux within the bulk vapor region.(c) Explanation of heat transfer at the solid-vapor interface, where heat transfer takes place through molecular reflection.(d) Heat transfer at the liquid-vapor interface, with evaporation, condensation, and reflection fluxes.

FIG. 3 .
FIG. 3. Contour plots showing distributions of temperature (left) and entropy generation inside a nanopore during superheating at (a) 6.00 V and (b) 6.65 V, intermittent film boiling at (c) 6.83 V, and steady-state film boiling at (d) 8.02 V, (e) 8.59 V, and (f) 9.04 V.The dotted line in (f) shows the smaller bubble shape at 8.02 V.The vectors in each subplot denote the heat flux.(g) System dissipation function calculated through continuum simulations.(h) Variations in height, contact angle, and mean free path of vapor molecules within the equilibrium film bubble for changing bias voltage.
)-S5(f) in the Supplemental Material[30]  show the variations of vapor pressure, temperature, density, net heat transfer, interfacial temperature drops, and condensation coefficients of the film bubble for different equilibrium sizes θ eqm .As Joule heating is greater near the pore cen-ter in the axial direction and greatest at the bubble tip in the radial direction [Fig.S8(c) in the Supplemental Material[30]], the highest thermal gradient [Fig.S8(d)] occurs at the bubble tip, resulting in a higher heat and mass inflow into the bubble at the tip.On the other hand, near the bubble base, there is condensation together with partial heat outflux [Fig.3(d)].The remaining heat input from the liquid is conducted to the surface of the silicon nitride membrane.The fact that the liquid temperature at the LV interface is nonuniform means that there are places where heat and mass flow in and out of the LV interface simultaneously.This is captured mathematically [Eq.(5)] by the spatial variation of the condensation coefficient [Fig.S5(h) in the Supplemental Material[30]].

FIG. 4 .
FIG. 4. (a)-(d) Nucleate bubble dynamics.(a) Simulated radial temperature distribution variation during the 98 ns bubble growth-collapse cycle at 6.75 V. (b) and (c) Calculated dissipation functions according to simulations during periodic nucleate boiling at 6.75 V and 6.83 V, respectively.(d) Variations of blockage durations and waiting time with voltage as seen in experiments.The error bars denote the standard deviations.(e)-(g) Unstable film bubble characteristics.(e) Mean current variation with bias voltage near the second bifurcation point.(f) Current spectrogram of ramp voltage pulse, P1.(g) Variations of simulated nanopore mean current (−△−) and system dissipation (− ⋄ −) with bias voltage for θ = 90°, 94°, 98°, and 102°bubbles in a thermal steady state.Darker shades of blue and green correspond to larger film bubbles.The pink lines shows the variations of current and dissipation when the bubble size also satisfies mechanical equilibrium.The black line indicates the rise in dissipation when the mean bubble is out of mechanical equilibrium and only in a thermal steady state, albeit satisfying constant mean current.