# Dark solitons under higher order dispersion

Tristram J. Alexander  G. A. Tsolias  Department of Mathematics and
Statistics, University of Massachusetts, Amherst, MA 01003-4515, United States
of America  A. Demirkaya  Mathematics Department, University of Hartford, 200
Bloomfield Ave., West Hartford, CT 06117, United States of America  Robert J.
Decker  Mathematics Department, University of Hartford, 200 Bloomfield Ave.,
West Hartford, CT 06117, United States of America  C. Martijn de Sterke
Institute of Photonics and Optical Science (IPOS), School of Physics, The
University of Sydney, NSW 2006, Australia  P. G. Kevrekidis  Department of
Mathematics and Statistics, University of Massachusetts, Amherst, MA
01003-4515, United States of America

###### Abstract

We show theoretically that dark solitons can exist in the presence of pure
quartic dispersion, and also in the presence of both quadratic and quartic
dispersive effects, displaying a much greater variety of possible solutions
and dynamics than for pure quadratic dispersion. The interplay of the two
dispersion orders may lead to oscillatory non-vanishing tails, which enables
the possibility of bound, potentially stable, multi-soliton states. Dark
soliton-like states which connect to low amplitude oscillations are also shown
to be possible. Dynamical evolution results corroborate the stability picture
obtained, and possible avenues for dark soliton generation are explored.

††journal: ol

Introduction. Following the experimental observation of bright solitons
supported by fourth-order dispersion [1], the theory of bright solitons with
pure-quartic dispersion [2], mixed second- and fourth-order dispersion [3] and
pure dispersion of higher even orders [4] has become well established. The
experimental work was enabled by the development of a fiber laser with
conveniently programmable dispersion [4, 5]. All of these bright solitons
exist when the highest-order of dispersion is negative. Concurrently, this has
inspired substantial mathematical interest in the analysis of such models [6,
7, 8, 9].

To the best of our knowledge, solitons with normal quartic dispersion, i.e.,
dark quartic solitons, have not been investigated. Dark solitons appear as
dark intensity dips on a continuous wave (CW) background [10, 11, 12].
Associated with the intensity minimum is a phase change of $\pi$, and this
“kink" in the phase provides the dark soliton with added (topological)
robustness against noise [13]. Recently, methods using self-induced
modulational instability with normal dispersion have allowed the generation of
dark soliton trains in optical fiber cavities [14], with evidence that dark
solitons are ubiquitous in these conditions [15]. Indicative of their broader
relevance is their emergence in both atomic systems [12, 16], and in water
waves [17, 18].

The effect of dominant higher-order normal dispersion was considered in the
$\phi^{4}$ Klein-Gordon model, in both the pure quartic form [19, 20] and in
the presence of both quadratic and quartic terms [21]. These results for real
fields point to intricate soliton interactions. Dark solitons bear a non-zero
background, unlike bright solitons, and, as we show, the quartic dispersion
leads to oscillations in their non-vanishing tails. This leads to a complex
variety of effects absent in the purely quadratic case, accompanied by rich
instability dynamics.

In this work we consider the fourth-order normal dispersion regime, in the
presence of both normal and anomalous second-order dispersion and cubic Kerr
nonlinearity. In the pure quadratic case there is a single dark soliton family
of stationary solutions, with multiple dark solitons repelling one another. We
find many more possibilities in the mixed quartic-quadratic dispersion case,
opening up new directions for possible dark soliton experiments. We identify
families of (multi-)dark solitons, examine their stability and instability
numerically and dynamically, and conclude with a demonstration of possible
dark soliton generation.

Model & Theoretical Background. We use the generalized nonlinear Schrödinger
equation for the electric field envelope $\tilde{\Psi}$ with quadratic and
quartic dispersion and a Kerr nonlinearity:

$i\frac{\partial\tilde{\Psi}}{\partial\xi}+\frac{\tilde{\beta}_{4}}{24}\frac{\partial^{4}\tilde{\Psi}}{\partial\tau^{4}}-\frac{\tilde{\beta}_{2}}{2}\frac{\partial^{2}\tilde{\Psi}}{\partial\tau^{2}}+\gamma|\tilde{\Psi}|^{2}\tilde{\Psi}=0$
(1)

where $\xi$ is the propagation distance, $\tau$ is the retarded time in the
frame of the pulse and $\gamma$ is the nonlinear coefficient, which we shall
take to be positive. The parameters $\tilde{\beta}_{2}=dv_{g}^{-1}/d\omega$
and $\tilde{\beta}_{4}=d^{3}v_{g}^{-1}/d\omega^{3}$, where $v_{g}$ is the
group velocity and $\omega$ is the pulse carrier frequency, characterize the
quadratic and quartic dispersion strengths respectively. We normalize the
retarded time in units of $t_{0}=1$ ps, $t=\tau/t_{0}$, and similarly the
propagation length in terms of a characteristic propagation length $z_{0}=1$
mm, $z=\xi/z_{0}$. Finally, we rescale $\Psi=\sqrt{z_{0}\gamma}\tilde{\Psi}$.
We restrict ourselves to positive (normal) quartic dispersion, and
consequently fix $\tilde{\beta}_{4}=1$ ps4 mm-1, consistent with experiment
[1], leading to $\beta_{4}=+1$. The resulting normalized model takes the form:

$i\frac{\partial\Psi}{\partial z}+\frac{1}{24}\frac{\partial^{4}\Psi}{\partial
t^{4}}-\frac{\beta_{2}}{2}\frac{\partial^{2}\Psi}{\partial
t^{2}}+|\Psi|^{2}\Psi=0,$ (2)

where the normalized quadratic dispersion parameter is given by
$\beta_{2}=\tilde{\beta}_{2}/(z_{0}t_{0}^{2})$. We look for stationary
solutions $\Psi(t,z)=\psi(t)\exp(i\mu z)$, where $\mu$ is a nonlinearity-
induced phase shift characterizing the stationary solution, leading to the
following equation for the real (without loss of generality) amplitude
$\psi(t)$:

$-\mu\psi+\frac{1}{24}\frac{d^{4}\psi}{dt^{4}}-\frac{\beta_{2}}{2}\frac{d^{2}\psi}{dt^{2}}+\psi^{3}=0.$
(3)

We investigate dark soliton solutions to Eq. (3), connecting to the continuous
wave background $\psi_{cw}\equiv\pm\sqrt{\mu}$ (which implies: $\mu>0$). In
the four dimensional phase space characterizing solutions to Eq. (3) these
dark soliton solutions may either connect the positive CW solution back to
itself (homoclinic solutions, for even numbers of solitons), or connect CW
solutions of different signs (heteroclinic solutions, for odd numbers of
solitons). We apply periodic boundary conditions so we only find homoclinic
solutions. We characterize these using the complementary power, which for a
cavity of length $2L$ is $Q_{c}=\int_{-L}^{L}(\psi_{cw}^{2}-\psi^{2})dt$.

As a starting point, we analyze the nature of the CW, on top of which the dark
solitons are built. We consider purely real perturbations in the form
$\psi=[\psi_{0}+\epsilon\exp(\lambda t)]$, and upon substitution of this
perturbation into Eq. (3) and keeping only linear terms in $\epsilon$ we
obtain the following equation for the perturbation eigenvalue $\lambda$:

$\frac{1}{24}\lambda^{4}-\frac{\beta_{2}}{2}\lambda^{2}-\mu+3\psi_{0}^{2}=0.$
(4)

Setting $\psi_{0}=\psi_{cw}$ and solving for $\lambda$ we find that if
$\mu\leq 3\beta_{2}^{2}/4$ we have four imaginary values for $\lambda$ if
$\beta_{2}<0$, or four real values if $\beta_{2}>0$. The former implies that
the CW is a center, so it is not possible to approach the CW and therefore no
dark solitons connecting to the CW can exist. Conversely, when $\lambda$ is
purely real, the CW is a saddle point and we have the more familiar dark
soliton regime. If instead $\mu>3\beta_{2}^{2}/4$ the $\lambda$ form a complex
quartet, so the CW is a saddle-spiral: any approach to the CW is accompanied
by oscillations. This behaviour is analogous to that observed for the bright
soliton (where instead $\beta_{4}<0$). The dependence of the CW solution on
the 2D parameter plane is shown in Fig. 1(a); the quadratic dependence of
$\mu$ on $\beta_{2}$ (solid curve), is as expected from the model’s scaling
properties.

As we shall see when we look at possible dark soliton solutions, the nature of
the zero solution can also play a role in the form of the stationary states.
Substituting $\psi_{0}=0$ into Eq. (4) the eigenvalues $\lambda$ are an
imaginary pair and a real pair, which makes the zero solution a saddle-center.
It is therefore not possible to approach zero asymptotically when
$\beta_{4}>0$, but it is possible to connect to an oscillatory solution about
zero. This behaviour is distinct from the case of $\beta_{4}<0$ [2, 3, 9]
where the zero solution is either approached asymptotically, or not at all.

While the analysis of the eigenvalues (4) indicates where we can expect
homoclinic (multi-)dark soliton solutions, a necessary condition for stability
of these solutions is that the CW is modulationally stable. To this end, we
examine the growth of linear waves by more generally perturbing the CW
background:

$\displaystyle\Psi(t,z)=\big{[}\psi_{cw}$
$\displaystyle+\epsilon_{1}\exp(ikt)\exp(i\Omega z)$ (5)
$\displaystyle+\epsilon_{2}^{*}\exp(-ikt)\exp(-i\Omega z)\big{]}\exp(i\mu z),$

where the $\epsilon_{i}$ denote small perturbations. Substitution of (5) into
Eq. (2) shows how the perturbation eigenvalue $\Omega$ is related to the
perturbation wavenumber $k$:

$\Omega^{2}=\left(-\mu+\frac{\beta_{2}}{2}k^{2}+\frac{1}{24}k^{4}+2\psi_{cw}^{2}\right)^{2}-\psi_{cw}^{4}.$
(6)

If $\Omega^{2}<0$ then the perturbation undergoes exponential growth, with the
boundary of modulational stability/instability occurring when $\Omega^{2}=0$.
We see that we always have an acoustic branch of modulational instability if
$\beta_{2}<0$ (see Fig. 1(b) for $\beta_{2}=-0.2$). Instability in this region
proceeds with characteristic modulation of the background (Fig. 1(c)) and
exponential growth of the most unstable wavenumber $k_{m}$ (Fig. 1(d)). We see
almost perfect agreement between the growth rate predicted by Eq. (6) and the
numerically observed value. The CW is modulationally stable if $\beta_{2}\geq
0$ and $\mu>0$.

Figure 1: (a) Classification of CW solution with $\beta_{2}$ and $\mu$; (b)
Modulational instability spectrum of CW for $\beta_{2}=-0.2$, with most
unstable wavenumber $k_{m}=1.1$ corresponding to $\Omega=0.77i$; (c)
Instability dynamics of CW for $\beta_{2}=-0.2$; and (d) Dependence of
spectral intensity $|\tilde{\psi}|^{2}$ on time at wavenumber $k_{m}$ for
results in (c), with slope giving $\Omega=0.77i$, agreeing very well with the
prediction. Figure 2: Pure quartic dark soliton stationary solutions
($\beta_{2}=0$, $\mu=5$). (a) A well separated pair of dark solitons; (b) A
mixed solution connecting dark solitons to an oscillation about $0$.

Numerical Findings All propagation results (including the CW dynamics in Fig.
1) are obtained using a fourth-order split-step numerical scheme [22] applied
to Eq. (2). Here, $\Delta t=9.8\times 10^{-3}$ with $\Delta z=7.6\times
10^{-6}$. The robustness of the scheme is monitored by evaluating the
Hamiltonian
$H=\int|d^{2}\psi/dt^{2}|^{2}+\beta_{2}|d\psi/dt|^{2}+|\psi|^{4}dt$ and
verifying that this quantity is conserved during propagation. Stationary
solutions are found by numerically solving Eq. (3) using a conjugate gradient
method [22].

The linear analysis points to new possible features in dark soliton solutions.
Figure 2 shows the two extremes of possible solutions in the pure-quartic
case. Well separated dark solitons can be found (Fig. 2(a)), similar to those
observed in the quadratic case, but with characteristic damped oscillations
approaching the CW, consistent with earlier results for the real case [20].
Fig. 2(b) shows a new possibility enabled by the quartic dispersion, a
connection between the CW and the saddle-centre at the origin.

Figure 3: (a) Bifurcation diagram for the lowest order dark soliton families
with the labels corresponding to solutions shown in the lower panels. (b) and
(c) solutions corresponding to family 0, ultimately connecting to a plane wave
as $\beta_{2}$ decreases; (d),(e),(f) and (g) are solutions from family 1;
(h),(i),(j) and (k) solutions from family 2. Left and right panels correspond
to $\beta_{2}=-1$, $\beta_{2}=1$, respectively. Unstable/stable solutions
shown with dashed/solid lines respectively. In all cases $\mu=5$.

Examining the possible solutions more systematically as a function of a system
parameter, e.g. $\beta_{2}$ in Fig. 3(a), we find that many possible families
exist, each with a characteristic spacing between the dark soliton pairs. The
most closely spaced dark soliton pairs, which we call ‘family $0$’ (Figs. 3(b)
and (c)), connect the CW to the origin, and we find bifurcate from the CW
solution at just less than $\beta_{2}=-2.5$. In contrast, all other families
of increasingly widely spaced dark soliton pairs (e.g. ‘family $1$’ and
‘family $2$’ in Figs. 3(e) and (i) respectively), develop additional
undulations as they proceed to negative $\beta_{2}$ (Figs. 3(d) and (h)),
where they collide (in a saddle-center bifurcation) with an upper branch
associated with large side oscillations (Fig. 3(f) and (j)). These upper
branch solutions have additional dark soliton pairs on either side of the main
pair, i.e., they are a triple pair family, which becomes increasingly evident
as $\beta_{2}$ increases (Fig. 3(g) and (k)). These results are as expected
from the earlier real analysis, with the oscillatory tails enabling isolated
distances at which dark solitons become stationary.

We examine the stability of the solutions shown in Fig. 3 by perturbing the
stationary solution and numerically computing the linearization eigenvalues
[22]. A full presentation of the linear stability analysis is beyond the scope
of this work, however we find results consistent with the real case [20]. We
find that family $0$ and family $2$ are stable for $\beta_{2}\geq 0$ (Fig.
3(c)), so stability alternates with spacing. This is topologically justified
from the underlying alternating minima and maxima in the relevant effective
landscape (c.f., e.g., Fig. 5 of [21] in the real-field case). All solutions
are unstable for $\beta_{2}<0$ due to the modulational instability of the CW
background. The continuous spectrum, as is typical for NLS dark solitons,
spans the entire imaginary axis [16]. Stable and unstable families of
solutions are shown respectively as solid and dashed lines in Fig. 3.

Figure 4: Numerical propagation of dark soliton solutions from different
families at $\beta_{2}=0$. (a) Family 1 is unstable to fission; (b) Family 2
is stable for $\beta_{2}>0$, due to inversion of the stability one step
further apart; (c) Family 2 upper branch, all upper branch composite solutions
are unstable; (d) Solution shown in Fig. 2(b), a very weak oscillatory
instability is present.

Examples of possible stability and instability dynamics are shown in Fig. 4
for $\beta_{2}=0$. Figs. 4(a) and (b) confirm respectively the instability of
family $1$ and the stability of family $2$. While all upper branch solutions
are unstable, the instability dynamics show long-lived non-stationary soliton
dynamics which, as we shall see, appear to play a significant role in dark
soliton generation processes. The instability dynamics shown in Fig. 4(c) are
one example, with more closely bound pairs appearing robust until disturbed by
the instability products of the inner pair. The generalisation of family $0$
to multiple oscillations about the zero solution, corresponding to Fig. 2(d),
appears remarkably robust, as seen in Fig. 4(d), but stability analysis
reveals a weak oscillatory instability.

Given the complexities of controlling dispersion (e.g., using a WaveShaper) in
the bright soliton case[5, 4], the careful pulse control of early dark soliton
generation experiments [11] appears to be practically challenging. We
therefore restrict considerations to an, arguably, prototypical initial
condition of a CW state with intensity notches, either of Gaussian form or
domain wall-like form using a hyperbolic tangent. With a Gaussian initial
condition we consider the dynamics in three different regimes: the pure
quartic case ($\beta_{2}=0$, Fig. 5(a)), positive $\beta_{2}$
($\beta_{2}=0.5$, Fig. 5(b)), and negative $\beta_{2}$ ($\beta_{2}=-0.2$, Fig.
5(c)). In the pure-quartic case a dark soliton pair appears, and exhibits a
periodic cycle of collisions and separations, before escaping the local
potential induced through the oscillatory tails, whereas for the positive
$\beta_{2}$ example, the dark solitons repel each other (which we note is the
only possible dynamics in the pure positive quadratic dispersion case). For
the negative $\beta_{2}$ case the dark soliton pair is more tightly bound, but
sits on a modulationally unstable CW background. A domain wall-like initial
condition produces multiple dark soliton pairs, and even what appears to be a
travelling complex of dark solitons, Fig. 5(d), similar to the instability
dynamics observed earlier. We see that in the presence of instability (soliton
or CW), solitons emerge, although they may be moving solitons. The nature of
these moving solitons is an interesting direction for further research,
particularly in light of recent results indicating the non-Galilean invariance
of solitons in the presence of quartic dispersion [23].

Figure 5: Dark soliton generation from intensity notches in a CW background.
(a,b,c) Initial condition $\psi(t,0)=\sqrt{\mu}(1-\exp(-t^{2}/2))$, (d)
initial condition $\psi(t,0)=\sqrt{\mu}(1+(\tanh(t-5)-\tanh(t+5))/2)$, all
with $\mu=5$; (a) Pure quartic dark soliton bound state with $\beta_{2}=0$;
(b) Repulsive dark soliton pair ($\beta_{2}=0.5$); (c) Dark soliton pair when
the CW is modulationally unstable ($\beta_{2}=-0.2$); (d) Dark soliton complex
at $\beta_{2}=0$.

Typical experimental values for associated optical settings would involve a
nonlinearity coefficient of $\gamma=4.07$ W-1mm-1 [1]. In this case, taking
the temporal units to be $1\leavevmode\nobreak\ {\rm ps}$ and the propagation
units to be $1\leavevmode\nobreak\ {\rm mm}$, the power of the CW considered
throughout is $1.2\leavevmode\nobreak\ {\rm W}$, the pulse length is
$40\leavevmode\nobreak\ {\rm ps}$, and the propagation distance is
$100\leavevmode\nobreak\ {\rm mm}$. These parameters can be adjusted via the
scaling properties of the model.

Conclusions & Outlook. This work provides the first study of higher-order
dispersion dark solitons in an optical context, focusing on soliton families
and their properties in pure quartic and quadratic/quartic dispersion
settings. A first example of a mixed state dark soliton is given, in which the
CW is connected to oscillations about the zero solution. The complex
bifurcation picture uncovered for dark soliton pairs requires further study,
including examining the properties of the waveforms as the modulationally
unstable regime is approached. The challenge of generating these states
indicates a wealth of phenomena that could be further explored, including
soliton complexes and gray solitons. Our simulations indicate that even for
the simple unstable states, it appears possible to observe long-lived
interacting states in the dynamics. We find no stable stationary states in the
presence of positive $\beta_{4}$ and negative $\beta_{2}$, raising the
question of the nature of the field in this region of parameter space.

## References

  * [1] A. Blanco-Redondo, C. M. de Sterke, J. Sipe, T. F. Krauss, B. J. Eggleton, and C. Husko, Nature Communications 7, 10427 (2016).
  * [2] K. K. K. Tam, T. J. Alexander, A. Blanco-Redondo, and C. M. de Sterke, Opt. Lett. 44, 3306 (2019).
  * [3] K. K. K. Tam, T. J. Alexander, A. Blanco-Redondo, and C. M. de Sterke, Phys. Rev. A 101, 043822 (2020).
  * [4] A. F. J. Runge, Y. L. Qiang, T. J. Alexander, M. Z. Rafat, D. D. Hudson, A. Blanco-Redondo, and C. M. de Sterke, Phys. Rev. Research 3, 013166 (2021).
  * [5] A. F. J. Runge, D. D. Hudson, K. K. K. Tam, C. M. de Sterke, and A. Blanco-Redondo, Nature Photonics 14, 492 (2020).
  * [6] A. Demirkaya and M. Stanislavova, Discrete Contin. Dyn. Syst. B 24, 197 (2019).
  * [7] I. Posukhovskyi and A. Stefanov, Discrete Contin. Dyn. Syst. A 40, 4131 (2020).
  * [8] R. Parker and A. Aceves, Physica D: Nonlinear Phenomena 422, 132890 (2021).
  * [9] R. I. Bandara, A. Giraldo, N. G. R. Broderick, and B. Krauskopf, Phys. Rev. A 103, 063514 (2021).
  * [10] P. Emplit, J. Hamaide, F. Reynaud, C. Froehly, and A. Barthelemy, Optics Communications 62, 374 (1987).
  * [11] A. M. Weiner, J. P. Heritage, R. J. Hawkins, R. N. Thurston, E. M. Kirschner, D. E. Leaird, and W. J. Tomlinson, Phys. Rev. Lett. 61, 2445 (1988).
  * [12] D. J. Frantzeskakis, J. Phys. A-Math. Theor. 43, 213001 (2010).
  * [13] W. Zhao and E. Bourkoff, Opt. Lett. 14, 703 (1989).
  * [14] Y. F. Song, J. Guo, L. M. Zhao, D. Y. Shen, and D. Y. Tang, Opt. Lett. 39, 3484 (2014).
  * [15] D. Y. Tang, L. Li, Y. F. Song, L. M. Zhao, H. Zhang, and D. Y. Shen, Phys. Rev. A 88, 013849 (2013).
  * [16] P. Kevrekidis, D. Frantzeskakis, and R. Carretero-González, _The Defocusing Nonlinear Schrodinger Equation_ (Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015).
  * [17] A. Chabchoub, O. Kimmoun, H. Branger, N. Hoffmann, D. Proment, M. Onorato, and N. Akhmediev, Phys. Rev. Lett. 110, 124101 (2013).
  * [18] F. Tsitoura, U. Gietz, A. Chabchoub, and N. Hoffmann, Phys. Rev. Lett. 120, 224102 (2018).
  * [19] R. J. Decker, A. Demirkaya, N. S. Manton, and P. G. Kevrekidis, Journal of Physics A: Mathematical and Theoretical 53, 375702 (2020).
  * [20] R. J. Decker, A. Demirkaya, P. Kevrekidis, D. Iglesias, J. Severino, and Y. Shavit, Communications in Nonlinear Science and Numerical Simulation 97, 105747 (2021).
  * [21] G. A. Tsolias, R. J. Decker, A. Demirkaya, T. J. Alexander, and P. G. Kevrekidis, Journal of Physics A: Mathematical and Theoretical 54, 225701 (2021).
  * [22] J. Yang, _Nonlinear Waves in Integrable and Nonintegrable Systems_ , vol. 16 of _Mathematical modeling and computation_ (SIAM, 2010).
  * [23] J. Widjaja, E. Kobakhidze, T. R. Cartwright, J. P. Lourdesamy, A. F. J. Runge, T. J. Alexander, and C. M. de Sterke, Phys. Rev. A 104, 043526 (2021).


# Thermodynamic Bifurcations of Boiling in Solid-State Nanopores

Soumyadeep Paul  Corresponding author<EMAIL_ADDRESS>Department of Mechanical Engineering, The University of Tokyo, 7-3-1, Hongo,
Bunkyo-ku, Tokyo 113-8656, Japan  Yusuke Ito  Department of Mechanical
Engineering, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-8656,
Japan  Wei-Lun Hsu  Department of Mechanical Engineering, The University of
Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan  Hirofumi Daiguji
Corresponding author<EMAIL_ADDRESS>Department of Mechanical
Engineering, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-8656,
Japan

###### Abstract

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 much-needed fundamental insights into the chaotic transition
boiling regime, which is least understood.

Nanopore, Ionic Joule heating, Resistive pulse sensing, Boiling bifurcation,
System dissipation

††preprint: APS/123-QED

##  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 [28, 29, 5, 4, 14], which allows us to
capture vapor nanobubble dynamics at nanosecond resolutions.

Figure 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. $T_{\mathrm{c}}$ and
$T_{\mathrm{w}}$ 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.

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. LABEL:fig:a1 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 [29, 28, 5]. 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 resembles 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.

##  II Results and Discussion

###  II.1 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_{\mathrm{app}}=4.5+150t~{}[\mathrm{s}]$; Fig.
LABEL:fig:a2(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
[$\tau_{\mathrm{th}}\sim R_{\mathrm{p}}^{2}/D_{\mathrm{th}}=\mathcal{O}$(100
ns)] and the lifetime of nucleate bubbles [$t_{\mathrm{b}}=67$ ns, Fig. 1(b)]
is much shorter than the millisecond-order voltage pulse. Here,
$R_{\mathrm{p}}$ is the pore radius and $D_{\mathrm{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_{\mathrm{app}}$ on the horizontal axis
and the resulting Joule heat dissipation $V_{\mathrm{app}}\times 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.

###  II.2 Steady states

####  II.2.1 Superheating

As the voltage is increased from 4.5 V to 6.65 V, the nanopore current
increases continuously [Fig. LABEL:fig:a2(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_{\mathrm{c}}$) and pore walls
($T_{\mathrm{w}}$) during superheating [Fig. LABEL:fig:a2(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. LABEL:fig:a2(a)
and 1(e), respectively] fit well with the experimental results. According to
the simulations, we find that $T_{\mathrm{c}}$ and $T_{\mathrm{w}}$ exceed the
saturation temperature of 373.15 K at voltages lower than 5.5 V. However, no
bubble signals are seen until 6.7125 V, when $T_{\mathrm{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 [29, 28], 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_{\mathrm{c}}$ and
$T_{\mathrm{w}}$, which originates from the focused nature of Joule heating.

We also calculate the total system dissipation during superheating,
$\Psi=\int_{V}{\psi}\,dV$, where the local dissipation in the liquid volume
can be written as $\psi=Td_{i}s/dt=T\dot{s}_{\mathrm{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

$\frac{d_{i}s}{dt}=\mathbf{J_{\mathrm{u}}}\cdot\nabla\left(\frac{1}{T}\right)+\frac{1}{T}\mathbf{J}\cdot\mathbf{E}.$
(1)

Here, $\mathbf{J_{\mathrm{u}}}=k_{\mathrm{w}}\nabla T$ is the heat flux in the
liquid and $\mathbf{J}=\sigma\mathbf{E}$ is the ionic current, where
$k_{\mathrm{w}}$ is the thermal conductivity and $\sigma$ is the electrical
conductivity, and $\mathbf{E}=-\nabla\phi$ is the electric field generated in
the nanopore liquid due to the applied bias voltage
$\Delta\phi=V_{\mathrm{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
$\bm{\tau}\,\mathord{:}\,(\nabla\mathbf{u})\approx\mu({\partial u}/{\partial
x})^{2}\approx\mu({u}/{R_{\mathrm{p}}})^{2}=\mathcal{O}(${10}^{11}\text{\text{~{}}}\mathrm{W}\mathrm{/}\mathrm{m}^{3}$)$,
which is much smaller than the dissipation originating from Joule heating,
$\sigma{|\mathbf{E}|}^{2}=\mathcal{O}(${10}^{16}\text{\text{~{}}}\mathrm{W}\mathrm{/}\mathrm{m}^{3}$)$
[see Fig. LABEL:fig:a35(a) in the Supplemental Material [30]], and heat
conduction, $k_{\mathrm{w}}{\nabla
T}^{2}/T=\mathcal{O}(${10}^{15}\text{\text{~{}}}\mathrm{W}\mathrm{/}\mathrm{m}^{3}$)$
[see Fig. LABEL:fig:a35(b) in the Supplemental Material [30]]. Here,
$u=\mathcal{O}($1\text{\text{~{}}}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$)$ is
assumed, which is consistent with the nanopore electroosmotic velocity [35,
36]. Using Eq. (1), we calculate the system dissipation $\Psi$, which
increases monotonically with voltage during superheating [Fig. 3(g)]. The
breakdown of dissipation sources is shown in Fig. LABEL:fig:a2(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.

Figure 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.

####  II.2.2 Stable film boiling

Following the second bifurcation at 8 V, stable film boiling (SFB) commences.
During this regime, a torus-shaped 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. LABEL:fig:a2(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. LABEL:fig:a2(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 [$\theta$ 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. LABEL:fig:a35(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 $\theta$ is the only free geometric parameter that controls the bubble
shape. The larger the angle $\theta$, 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_{\mathrm{max}}\approx
3\lambda$, where $\lambda$ is the thermal mean free path of vapor molecules,
which satisfies
$\lambda=1.922\lambda_{\mathrm{HS}}=1.922m/(\sqrt{2}\,\pi\rho_{\mathrm{v}}d^{2})$
[37]. $\lambda_{\mathrm{HS}}$ is the mean free path per hard sphere model.
Here, $m$ is the mass of one vapor molecule, $\rho_{\mathrm{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_{\mathrm{LV}}\sim
h_{\mathrm{SV}}\sim\lambda$), and the vapor bulk region has a length of less
than one mean free path [$h_{\mathrm{bulk}}<\lambda$ 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_{\mathrm{avg}}\approx 1.16\lambda$. Hence, we consider temperature drops
across the LV and SV interfaces (net thickness $2\lambda$) to account for the
net temperature drop across the bubble.

Figure 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.

Considering the heat transfer at the SV interface to be governed by molecular
reflections [Fig. 2(c)], we express the heat flux as

$q_{\mathrm{m}}=-k_{\mathrm{m}}\frac{\partial T}{\partial
n}=q_{\mathrm{bulk}}=\rho_{\mathrm{v}}\sqrt{\frac{R_{\mathrm{g}}T_{\mathrm{v}}}{2\pi}}\,c_{\mathrm{v}}(T_{\mathrm{v}}-T_{\mathrm{SV}}),$
(2)

where $R_{\mathrm{g}}$ is the universal gas constant, and $k_{\mathrm{m}}$ is
the thermal conductivity of the silicon nitride membrane, taken as
$3.2\text{\text{~{}}}\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$
[28]. $q_{\mathrm{bulk}}$ is the incoming heat flux from the vapor bulk region
to the SV interface. $\rho_{\mathrm{v}}$, $c_{\mathrm{v}}$, and
$U_{\mathrm{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. LABEL:fig:a33]. As shown in Fig. 2(d), the heat
flux balance at the LV interface can be written as

$\displaystyle q_{\mathrm{in}}={}$ $\displaystyle
q_{\mathrm{evap}}-q_{\mathrm{cond}}+q_{\mathrm{ref}}+h_{\mathrm{fg}}\dot{m}$
$\displaystyle={}$
$\displaystyle\alpha_{\mathrm{e}}\rho_{\mathrm{v,sat}}U_{\mathrm{v,sat}}\sqrt{\frac{R_{\mathrm{g}}T_{\mathrm{LV}}}{2\pi}}-\alpha_{\mathrm{c}}\rho_{\mathrm{v}}U_{\mathrm{v}}\sqrt{\frac{R_{\mathrm{g}}T_{\mathrm{v}}}{2\pi}}$
$\displaystyle+(1-\alpha_{\mathrm{e}})\rho_{\mathrm{v,sat}}\sqrt{\frac{R_{\mathrm{g}}T_{\mathrm{LV}}}{2\pi}}\,c_{\mathrm{v}}(T_{\mathrm{LV}}-T_{\mathrm{v}})+h_{\mathrm{fg}}\dot{m}.$
(3)

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_{\mathrm{LV}}$. The latent heat of evaporation or condensation
($h_{\mathrm{fg}}\dot{m}$) 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_{\mathrm{in}}=k_{\mathrm{w}}\hat{n}\cdot\nabla T$ in Figs. 2(d) and 3(f)].
The remaining heat
($q_{\mathrm{int}}=q_{\mathrm{evap}}-q_{\mathrm{cond}}+q_{\mathrm{ref}}$) is
transferred through the molecular fluxes of evaporation, condensation, and
reflection at the LV interface [$J_{\mathrm{evap}}$, $J_{\mathrm{cond}}$, and
$J_{\mathrm{r,i}}$ as shown in Fig. 2(d)]. The molecular velocities in these
fluxes are defined based on Maxwellian distributions [40, 41, 42].
Accordingly, $\langle
J_{\mathrm{out}}\rangle=\rho_{\mathrm{v,sat}}\sqrt{{R_{\mathrm{g}}T_{\mathrm{LV}}}/{2\pi}}$
and $\langle
J_{\mathrm{coll}}\rangle=\rho_{\mathrm{v}}\sqrt{{R_{\mathrm{g}}T_{\mathrm{v}}}/{2\pi}}$,
alongwith $\langle J_{\mathrm{evap}}\rangle=\alpha_{\mathrm{e}}\langle
J_{\mathrm{out}}\rangle$ and $\langle
J_{\mathrm{cond}}\rangle=\alpha_{\mathrm{c}}\langle J_{\mathrm{coll}}\rangle$
[39]. The net evaporation/condensation mass flux across the LV interface can
then be written as

$\dot{m}=\alpha_{\mathrm{e}}\rho_{\mathrm{v,sat}}\sqrt{\frac{R_{\mathrm{g}}T_{\mathrm{LV}}}{2\pi}}-\alpha_{\mathrm{c}}\rho_{\mathrm{v}}\sqrt{\frac{R_{\mathrm{g}}T_{\mathrm{v}}}{2\pi}}\,,$
(4)

where $\alpha_{\mathrm{e}}$ and $\alpha_{\mathrm{c}}$ are the evaporation and
condensation accommodation coefficients, respectively. For the results shown
in Fig. 3, $\alpha_{\mathrm{e}}=0.9$ is uniformly applied on the LV interface.
Applying the local balance of fluxes on the LV interface ($\langle
J_{\mathrm{out}}\rangle-\langle J_{\mathrm{coll}}\rangle=\langle
J_{\mathrm{evap}}\rangle-\langle J_{\mathrm{cond}}\rangle$),
$\alpha_{\mathrm{c}}$ can be written as

$\alpha_{\mathrm{c}}=1-\frac{(1-\alpha_{\mathrm{e}})\rho_{\mathrm{v,sat}}}{\rho_{\mathrm{v}}}\sqrt{\frac{T_{\mathrm{LV}}}{T_{\mathrm{v}}}}\,.$
(5)

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:

$\displaystyle P_{\mathrm{v}}$ $\displaystyle=K\gamma+P_{\mathrm{w}},$ (6)
$\displaystyle\iint_{S_{\mathrm{LV}}}\dot{m}dS$ $\displaystyle=0,$
$\displaystyle\iint_{S_{\mathrm{LV}}}q_{\mathrm{in}}\,dS$
$\displaystyle=\iint_{S_{\mathrm{SV}}}q_{\mathrm{m}}\,dS.$

Here, $P_{\mathrm{v}}$ and $P_{\mathrm{w}}$ are the vapor and liquid
pressures, respectively. $S_{\mathrm{LV}}$ and $S_{\mathrm{SV}}$ are the
liquid–vapor and solid-vapor interfaces respectively as shown in Fig.
LABEL:fig:a22 (inset) in SI [30]. $\gamma$ is the surface tension of water
defined at $T_{\mathrm{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. LABEL:fig:a22] using the multiphysics software, arb [43].
This yields the temperature distribution of nanopores for varying bias voltage
[Figs. 3(d)–3(f)]. We find that the film bubble expands
[$\theta_{\mathrm{eqm}}$ increases from $117\degree$ to $125\degree$ as shown
in Figs. 3d–3(f) and 3(h)] to satisfy the closure equations (6) for increasing
voltage. The right panels in Figs. 3(d)–3(f) show the net entropy generation
in the bulk liquid and solid regions, which is a combination of Joule heating
[Fig. LABEL:fig:a35(c) in the Supplemental Material [30]] and thermal gradient
[Fig. LABEL:fig:a35(d)] sources according to Eq. (1). Figure 3(f) shows that
the bulging-out effect results in locally high Joule heating (high
$\dot{s}_{\mathrm{gen}}$) near the bubble tip and locally high heat inflow
(high $q_{\mathrm{in}}$). Figures LABEL:fig:a3(a)–LABEL:fig:a3(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
$\theta_{\mathrm{eqm}}$. As Joule heating is greater near the pore center in
the axial direction and greatest at the bubble tip in the radial direction
[Fig. LABEL:fig:a35(c) in the Supplemental Material [30]], the highest thermal
gradient [Fig. LABEL:fig:a35(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. LABEL:fig:a3(h) in the Supplemental Material [30]].

A stable nanobubble with the simultaneous influx and outflux of dissolved gas
on the bubble interface was termed dynamic equilibrium or controlled
nonequilibrium 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
$\theta_{\mathrm{eqm}}$ for different values of $\alpha_{\mathrm{e}}$ at the
same bias voltage are very close [blue and turquoise traces in Fig. 3(h)].
However, as shown in Figs. LABEL:fig:a3(f) and LABEL:fig:a38(f) in the
Supplemental Material [30], small values of $\alpha_{\mathrm{e}}$ cause
negative and unrealistic values of $\alpha_{\mathrm{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
$\alpha_{\mathrm{e}}=0.5$, when $\alpha_{\mathrm{e}}=0.9$, we find that
$\alpha_{\mathrm{c}}$ is less non-uniform and varies over a shorter range of
0.82-0.98 on the bubble surface (Fig. LABEL:fig:a3f 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,
$\alpha_{\mathrm{c}}$ can remain uniform on the LV interface, while
$\alpha_{\mathrm{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 $\alpha_{\mathrm{c}}=0.9$ and obtained a spatially
varying $\alpha_{\mathrm{e}}$ ranging from 0.51 to 0.94 (Fig. LABEL:fig:a32c).
When uniform $\alpha_{\mathrm{c}}=0.9$ is applied on the LV interface, the
equilibrium size of $\theta_{\text{eqm}}=117\degree$ is obtained at 8.04 V
(Fig. LABEL:fig:a32a). On the other hand, when uniform
$\alpha_{\mathrm{e}}=0.9$ is applied on the LV interface, equilibrium size of
$\theta_{\mathrm{eqm}}=117\degree$ is obtained at 8.02 V (Fig. 2d). Also, the
vapor temperature is 415.91 K when $\alpha_{\mathrm{e}}=0.9$ (Fig.
LABEL:fig:a3b) and 416.33 K (Fig. LABEL:fig:a32b) when
$\alpha_{\mathrm{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. LABEL:fig:a7 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.

###  II.3 Unsteady states

####  II.3.1 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. LABEL:fig:a21 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 [22, 5]. We simulated the spherical bubble growth–collapse cycle
through a one-dimensional moving boundary model [46, 4] 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 one-dimensional 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. LABEL:fig:a2c], 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. LABEL:fig:a22]. 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
period-averaged system dissipation during the nucleate boiling cycle also
decreased with increasing bias voltage (blue trace in Fig. LABEL:fig:a6 in the
Supplemental Material [30]).

Figure 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 ($-\triangle-$) and
system dissipation ($-\diamond-$) with bias voltage for $\theta=90\degree$,
$94\degree$, $98\degree$, and $102\degree$ 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.

####  II.3.2 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_{\mathrm{SV}}=T_{\mathrm{v}}$ instead of Eq. (2). Accordingly, a steady-
state current of $31\text{\text{~{}}}\mathrm{\SIUnitSymbolMicro 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_{\mathrm{avg}}>\lambda$ during unstable
film boiling (UFB) [Fig. 3(h)], we estimate the temperature drop in the vapor
bulk, $\Delta T_{\mathrm{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_{\mathrm{ave}}=1\lambda$, $\Delta T_{\mathrm{bulk}}$ increases as
the bubble expands with increasing voltage [Fig. LABEL:fig:a38(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_{\mathrm{ave}}<\lambda$,
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. LABEL:fig:a6) was obtained
by a weighted average of the nucleate boiling dissipation [Figs. 4(b) and
4(c)] and the unstable film bubble dissipation [Figs. LABEL:fig:a7(c) and
LABEL:fig:a7(d)], based on their relative probabilities of occurrence
$P_{\mathrm{NB}}$ and $P_{\mathrm{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_{\mathrm{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_{\mathrm{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. LABEL:fig:a5 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 [55, 14,
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_{\mathrm{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. LABEL:fig:a5),
and the current drops 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. LABEL:fig:a5). In the UFB region, the bubble swelling
increases $\Delta T_{\mathrm{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
extrapolated UFB mean current and the SFB current (jump 2). This suggests that
a membrane or film bubble with $h_{\mathrm{avg}}$ in the range from $\lambda$
to $2\lambda$ (at UFB) follows a different heat transport and stabilization
mechanism than one with $h_{\mathrm{avg}}=2\lambda$ (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 entropy-generating 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.

###### Acknowledgements.

The authors would like to gratefully acknowledge Juan G Santiago at Stanford
University for his valuable comments and discussions on the manuscript. This
work was supported by the Japan Society for the Promotion of Science (JSPS)
KAKENHI Grant Nos. 20H02081 and 20J22422. Part of this work was conducted at
the Advanced Characterization Nanotechnology Platform of the University of
Tokyo, which was supported by the Nanotechnology Platform of the Ministry of
Education, Culture, Sports, Science, and Technology (MEXT), Japan, Grant No.
JPMXP09A21UT0089.

## References

  * Liu and Mauter [2022] T. Liu and M. S. Mauter, Heat Transfer Innovations and their Application in Thermal Desalination Processes, Joule 6, 1199 (2022).
  * Chen  _et al._ [2009] R. Chen, M. C. Lu, V. Srinivasan, Z. Wang, H. H. Cho, and A. Majumdar, Nanowires for Enhanced Boiling Heat Transfer, Nano Lett. 9, 548 (2009).
  * Dhillon  _et al._ [2015] N. S. Dhillon, J. Buongiorno, and K. K. Varanasi, Critical Heat Flux Maxima during Boiling Crisis on Textured Surfaces, Nat. Commun. 6, 8247 (2015).
  * Paul  _et al._ [2021] S. Paul, W.-L. Hsu, M. Magnini, L. R. Mason, Y. Ito, Y.-L. Ho, O. K. Matar, and H. Daiguji, Analysis and Control of Vapor Bubble Growth inside Solid-State Nanopores, J. Therm. Sci. Technol. 16, JTST0007 (2021).
  * Paul  _et al._ [2020] S. Paul, W.-L. Hsu, M. Magnini, L. R. Mason, Y.-L. Ho, O. K. Matar, and H. Daiguji, Single-Bubble Dynamics in Nanopores: Transition between Homogeneous and Heterogeneous Nucleation, Phys. Rev. Res. 2, 043400 (2020).
  * Thome and Dupont [2007] J. R. Thome and V. Dupont, Bubble Generator (2007), US Patent 7,261,144.
  * Asai [1989] A. Asai, Application of the Nucleation Theory to the Design of Bubble Jet Printers, Jpn. J. Appl. Phys. 28, 909 (1989).
  * Jiang  _et al._ [2022] M. Jiang, Y. Wang, F. Liu, H. Du, Y. Li, H. Zhang, S. To, S. Wang, C. Pan, J. Yu, _et al._ , Inhibiting the Leidenfrost Effect above 1,000°C for Sustained Thermal Cooling, Nature 601, 568 (2022).
  * Dhir [1998] V. K. Dhir, Boiling Heat Transfer, Annu. Rev. Fluid Mech. 30, 365 (1998).
  * Zhang  _et al._ [2023] L. Zhang, C. Wang, G. Su, A. Kossolapov, G. Matana Aguiar, J. H. Seong, F. Chavagnat, B. Phillips, M. M. Rahman, and M. Bucci, A Unifying Criterion of the Boiling Crisis, Nat. Commun. 14, 2321 (2023).
  * Bureš and Sato [2021] L. Bureš and Y. Sato, On the Modelling of the Transition between Contact-Line and Microlayer Evaporation Regimes in Nucleate Boiling, J. Fluid Mech. 916, A53 (2021).
  * Thome [2004] J. R. Thome, Boiling in Microchannels: A Review of Experiment and Theory, Int. J. Heat Fluid Flow 25, 128 (2004).
  * Lavino  _et al._ [2021] A. D. Lavino, E. Smith, M. Magnini, and O. K. Matar, Surface Topography Effects on Pool Boiling via Non-Equilibrium Molecular Dynamics Simulations, Langmuir 37, 5731 (2021).
  * Paul  _et al._ [2022] S. Paul, W.-L. Hsu, Y. Ito, and H. Daiguji, Boiling in Nanopores through Localized Joule Heating: Transition between Nucleate and Film Boiling, Phys. Rev. Res. 4, 043110 (2022).
  * Popov  _et al._ [2018] E. Popov, L. He, E. Dominguez-Ontiveros, and Y. Melnichenko, Detection of Vapor Nanobubbles by Small Angle Neutron Scattering (SANS), Appl. Phys. Lett. 112, 153704 (2018).
  * Lombard  _et al._ [2016] J. Lombard, T. Biben, and S. Merabia, Ballistic Heat transport in Laser Generated Nano-Bubbles, Nanoscale 8, 14870 (2016).
  * Maheshwari  _et al._ [2018] S. Maheshwari, M. van der Hoef, A. Prosperetti, and D. Lohse, Dynamics of Formation of a Vapor Nanobubble around a Heated Nanoparticle, J. Phys. Chem. C 122, 20571 (2018).
  * Hołyst  _et al._ [2017] R. Hołyst, M. Litniewski, and D. Jakubczyk, Evaporation of Liquid Droplets of Nano- and Micro-Meter Size as a Function of Molecular Mass and Intermolecular Interactions: Experiments and Molecular Dynamics Simulations, Soft Matter 13, 5858 (2017).
  * Lohse and Zhang [2015] D. Lohse and X. Zhang, Surface Nanobubbles and Nanodroplets, Rev. Mod. Phys. 87, 981 (2015).
  * Seddon  _et al._ [2011] J. R. T. Seddon, H. J. W. Zandvliet, and D. Lohse, Knudsen Gas Provides Nanobubble Stability, Phys. Rev. Lett. 107, 116101 (2011).
  * Brenner and Lohse [2008] M. P. Brenner and D. Lohse, Dynamic Equilibrium Mechanism for Surface Nanobubble Stabilization, Phys. Rev. Lett. 101, 214505 (2008).
  * Prosperetti [2016] A. Prosperetti, Vapor Bubbles, Annu. Rev. Fluid Mech. 49, 221 (2016).
  * Nicolis and Prigogine [1989] G. Nicolis and I. Prigogine,  _Exploring Complexity: An Introduction_ (W. H. Freeman, 1989).
  * Shoji [2004] M. Shoji, Studies of Boiling Chaos: A Review, Int. J. Heat Mass Transf. 47, 1105 (2004).
  * Chai and Shoji [2001] L. Chai and M. Shoji, Boiling Curves—Bifurcation and Catastrophe, Int. J. Heat Mass Transfer 44, 4175 (2001).
  * Chai and Shoji [2002] L. H. Chai and M. Shoji, Self-Organization and Self-Similarity in Boiling Systems, J. Heat Transf. 124, 507 (2002).
  * Nicolis and Prigogine [1977] G. Nicolis and I. Prigogine,  _Self-Organization in Nonequilibrium Systems: From Dissipative Structures to Order Through Fluctuations_, A Wiley-Interscience publication (Wiley, 1977).
  * Levine  _et al._ [2016] E. V. Levine, M. M. Burns, and J. A. Golovchenko, Nanoscale Dynamics of Joule Heating and Bubble Nucleation in a Solid-State Nanopore, Phys. Rev. E 93, 13124 (2016).
  * Nagashima  _et al._ [2014] G. Nagashima, E. V. Levine, D. P. Hoogerheide, M. M. Burns, and J. A. Golovchenko, Superheating and Homogeneous Single Bubble Nucleation in a Solid-State Nanopore, Phys. Rev. Lett. 113, 24506 (2014).
  * [30] See Supplemental Material at [URL will be inserted by publisher] for further information.
  * Prigogine [1978] I. Prigogine, Time, Structure, and Fluctuations, Science 201, 777 (1978).
  * Onsager [1931] L. Onsager, Reciprocal Relations in Irreversible Processes. I., Phys. Rev. 37, 405 (1931).
  * Demirel [2007] Y. Demirel,  _Nonequilibrium Thermodynamics: Transport and Rate Processes in Physical, Chemical and Biological Systems_ (Elsevier, 2007).
  * Avedisian [1985] C. T. Avedisian, The Homogeneous Nucleation Limits of Liquids, J. Phys. Chem. Ref. Data 14, 695 (1985).
  * Di Muccio  _et al._ [2022] G. Di Muccio, B. Morozzo della Rocca, and M. Chinappi, Geometrically Induced Selectivity and Unidirectional Electroosmosis in Uncharged Nanopores, ACS Nano 16, 8716 (2022).
  * Wang  _et al._ [2020] Z. Wang, W.-L. Hsu, S. Tsuchiya, S. Paul, A. Alizadeh, and H. Daiguji, Joule Heating Effects on Transport-Induced-Charge Phenomena in an Ultrathin Nanopore, Micromachines 11, 1041 (2020).
  * Dongari  _et al._ [2011] N. Dongari, Y. Zhang, and J. M. Reese, Modeling of Knudsen Layer Effects in Micro/Nanoscale Gas Flows, J. Fluids Eng. 133 (2011).
  * Wagner and Pruß [2002] W. Wagner and A. Pruß, The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use, J. Phys. Chem. Ref. Data 31, 387 (2002).
  * Cao  _et al._ [2011] B. Y. Cao, J. F. Xie, and S. S. Sazhin, Molecular Dynamics Study on Evaporation and Condensation of _n_ -Dodecane at Liquid–Vapor Phase Equilibria, J. Chem. Phys. 134, 10.1063/1.3579457 (2011).
  * Aursand and Ytrehus [2019] E. Aursand and T. Ytrehus, Comparison of Kinetic Theory Evaporation Models for Liquid Thin-Films, Int. J. Multiph. Flow 116, 67 (2019).
  * Kobayashi  _et al._ [2018] K. Kobayashi, T. Nagayama, M. Watanabe, H. Fujii, and M. Kon, Molecular Gas Dynamics Analysis on Condensation Coefficient of Vapour during Gas–Vapour Bubble Collapse, J. Fluid Mech. 856, 1045 (2018).
  * Kobayashi  _et al._ [2016] K. Kobayashi, K. Hori, M. Kon, K. Sasaki, and M. Watanabe, Molecular Dynamics Study on Evaporation and Reflection of Monatomic Molecules to Construct Kinetic Boundary Condition in Vapor–Liquid Equilibria, Heat Mass Transf. 52, 1851 (2016).
  * Harvie [2012] D. J. E. Harvie, An Implicit Finite Volume Method for Arbitrary Transport Equations, ANZIAM J. 52, 1126 (2012).
  * Liu  _et al._ [2017] Y. Liu, M. A. Edwards, S. R. German, Q. Chen, and H. S. White, The Dynamic Steady State of an Electrochemically Generated Nanobubble, Langmuir 33, 1845 (2017).
  * Kondepudi [2008] D. Kondepudi,  _Introduction to Modern Thermodynamics_ (Wiley, 2008).
  * Robinson and Judd [2004] A. J. Robinson and R. L. Judd, The Dynamics of Spherical Bubble Growth, Int. J. Heat Mass Transf. 47, 5101 (2004).
  * Sullivan  _et al._ [2022] P. Sullivan, D. Dockar, M. K. Borg, R. Enright, and R. Pillai, Inertio-Thermal Vapour Bubble Growth, J. Fluid Mech. 948, A55 (2022).
  * Zhong  _et al._ [2020] X. Zhong, J. Eshraghi, P. Vlachos, S. Dabiri, and A. M. Ardekani, A Model for a Laser-Induced Cavitation Bubble, Int. J. Multiph. Flow 132, 103433 (2020).
  * Zhang and Prosperetti [2021] Y. Zhang and A. Prosperetti, Dynamics, Heat and Mass Transfer of a Plasmonic Bubble on a Solid Surface, Int. J. Heat Mass Transf. 167, 120814 (2021).
  * Jou  _et al._ [2005] D. Jou, J. Casas-Vazquez, G. Lebon, and M. Grmela, A Phenomenological Scaling Approach for Heat Transport in Nano-Systems, Appl. Math. Lett. 18, 963 (2005).
  * Joshi and Majumdar [1993] A. A. Joshi and A. Majumdar, Transient Ballistic and Diffusive Phonon Heat Transport in Thin Films, J. Appl. Phys. 74, 31 (1993).
  * Dong  _et al._ [2012] Y. Dong, B. Y. Cao, and Z. Y. Guo, General Expression for Entropy Production in Transport Processes Based on the Thermomass Model, Phys. Rev. E 85, 1 (2012).
  * Sellitto and Cimmelli [2012] A. Sellitto and V. A. Cimmelli, A Continuum Approach to Thermomass Theory, J. Heat Transf. 134, 112402 (2012).
  * Guo and Wang [2015] Y. Guo and M. Wang, Phonon Hydrodynamics and Its Applications in Nanoscale Heat Transport, Phys. Rep. 595, 1 (2015).
  * Jenkins [2013] A. Jenkins, Self-Oscillation, Phys. Rep. 525, 167 (2013).
  * Nguyen  _et al._ [2018] D. M. Nguyen, L. Hu, J. Miao, and C.-D. Ohl, Oscillate Boiling from Electrical Microheaters, Phys. Rev. Appl. 10, 044064 (2018).
  * Li  _et al._ [2017] F. Li, S. R. Gonzalez-Avila, D. M. Nguyen, and C.-D. Ohl, Oscillate Boiling from Microheaters, Phys. Rev. Fluids 2, 014007 (2017).
  * Otsubo  _et al._ [2022] S. Otsubo, S. K. Manikandan, T. Sagawa, and S. Krishnamurthy, Estimating Time-Dependent Entropy Production from Non-Equilibrium Trajectories, Commun. Phys. 5, 11 (2022).
  * Van Kampen [1992] N. G. Van Kampen,  _Stochastic processes in physics and chemistry_ , Vol. 1 (Elsevier, 1992).
  * Noy and Darling [2023] A. Noy and S. B. Darling, Nanofluidic Computing Makes a Splash, Science 379, 143 (2023).
  * Marr and Munakata [2007] D. W. Marr and T. Munakata, Micro/nanofluidic Computing, Communications of the ACM 50, 64 (2007).
  * Karnik  _et al._ [2005] R. Karnik, R. Fan, M. Yue, D. Li, P. Yang, and A. Majumdar, Electrostatic Control of Ions and Molecules in Nanofluidic Transistors, Nano Lett. 5, 943 (2005).