Free Energy Evaluation of Homogeneous Bubble Nucleation based on Stochastic Thermodynamics

1. Introduction

In this research we focus on the energetics of homogeneous bubble nucleation in stretched metastable liquid. Inception of bubble nucleation is a fundamental process in the first order liquid-vapor phase transition in general, and has been investigated long in various fields, such as cavitation [1,2], nucleate boiling [3,4], and metastability limit of liquids [5,6,7,8], but the understanding of microscopic details are still required.

The mechanism of nucleation has been studied based mainly on the so-called classical nucleation theory (CNT) [5,6,9,10,11]. In a simplest form of CNT, the free energy change $\Delta F$ of the system during the growth of a spherical bubble is given by the following equation,

$$\Delta F\left(r\right)=-{\displaystyle \frac{4}{3}}\pi r^3\Delta f_{\mathrm{V}}+4\pi r^2\gamma$$

(1)

where r is the bubble radius, $\Delta f_{\mathrm{V}}$ the difference in free energy density between liquid and vapor phases, and $\gamma$ the surface tension. Although the physics for Eq. (1) is clear as a base of energetics of homogeneous bubble nucleation, discrepancy in the nucleation rate between CNT models and experiments is often extremely large [12], and various improvements in theoretical models have been proposed [6,11,13,14,15].

In this study, we revisit the free energy change in CNT of homogeneous bubble nucleation by utilizing molecular dynamics (MD) simulations combined with stochastic thermodynamics. Several MD simulations on homogeneous bubble nucleation were reported [16,17,18,19,20]. Due to the spatial and temporal scale limits of the simulation method, however, the conditions used in these MD simulations are highly non-equilibrium, closer to the spinodal lines, under which spontaneous bubble generation occurs within the simulation time. Furthermore, kinetics (or the nucleation rate) was the main target in these simulations while energetics is still hard to investigate because the system is in highly non-equilibrium state during this kind of simulations.

Here we want to investigate the free energy during the bubble nucleation. For that purpose, we adopt a basic formula of stochastic thermodynamics, with which we are able to evaluate the free energy from a special type of ensemble average of works to create and expand a bubble.

2. Stochastic Thermodynamics: Theoretical Background

Among several variations of free energy evaluation scheme proposed in the field of stochastic thermodynamics [21], we here adopt the simplest form, the Jarzynski equality. Assuming any path $\mathcal{L}$ in the configurational space from a starting point to a goal, this equality describes a relationship between the free energy $F\left(s\right)$ on a point $s\in \mathcal{L}$ ($s=0$ is the starting point) and the work $W\left(s\right)$ required along the path. If we can assume that the process along the path is always quasistatic, the conventional thermodynamics tells us

$$F\left(s\right)-F\left(0\right)=W\left(s\right)$$

(2)

but in general

$$F\left(s\right)-F\left(0\right)\le W\left(s\right)$$

(3)

which is a representation of the second law of thermodynamics. The difference is of course dissipated, as heat Q in most cases,

$$F\left(s\right)-F\left(0\right)=W\left(s\right)-Q\left(s\right)$$

(4)

In microscale processes where fluctuations are not negligible, W also contains fluctuations, and Eq. (3) should be expressed as

$$F\left(s\right)-F\left(0\right)\le 〈W\left(s\right)〉$$

(5)

where $\langle \rangle$ represents the ensemble average.

In 1997 [22], Jarzynski showed that the following equality holds for “any” processes

$$exp\left(-{\displaystyle \frac{F\left(s\right)-F\left(0\right)}{k_{\mathrm{B}}T}}\right)=〈exp\left({\displaystyle \frac{W\left(s\right)}{k_{\mathrm{B}}T}}\right)〉$$

(6)

where T is the temperature of thermal reservoir, and $k_{\mathrm{B}}$ is the Boltzmann constant. The point in Eq. (6) is that we can evaluate the state quantity $F\left(s\right)$ without the conventional assumption of quasistatic process. This gives a great merit for us to evaluate free energy with MD simulations because a typical time scale of MD simulations is so short that, in many cases, we cannot assume thermal equilibrium however slow the system change is.

In this work, we targeted simple liquid and conducted a series of MD simulations to observe a “bubble” generation by applying an external force. Then, by adopting Eq. (6), the free energy change during the bubble growth was evaluated. Since Eq. (6) was proposed, several extensions and generalizations have been proposed [21,23,24], such as the Crooks fluctuation theorem [25], which relates the probability distributions of the forward process and the backward one. These extended methods may provide a better tool to investigate the bubble nucleation, which is left for future study.

3. System and Method

We carried out a series of MD simulations to investigate the free energy for a tiny bubble in model liquid. All MD simulations were executed with LAMMPS [26,27]. OVITO [28] was utilized to visualize the atomic configurations.

3.1. Simulation System and Procedure

We adopted a simple monatomic liquid, in which the Lennard-Jones (LJ) 12-6 potential was assumed as the particle-particle interaction,

$$\varphi \left(r\right)=4\epsilon \left[{\left({\displaystyle \frac{\sigma }{r}}\right)}^{12}-{\left({\displaystyle \frac{\sigma }{r}}\right)}^6\right],$$

(7)

where r is the particle-particle distance. $\epsilon$ and $\sigma$ are the parameters for energy and length, respectively. The interaction is truncated at $r=r_c=3.5\sigma$. In the following descriptions, all quantities are expressed with the conventional reduced units, such as $\epsilon$:energy, $\sigma$:length, and m:particle mass. Just for reference, typical values for argon [29] are $\epsilon \simeq 1.66\times {10}^{-21}$ J, $\sigma \simeq 0.34\times {10}^{-9}$ m, and the time unit $\tau \equiv \sigma \sqrt{m/\epsilon }\simeq 2.15\times {10}^{-12}$ s. All MD simulations were executed with time step $0.002\tau$.

With conventional NPT ensemble conditions, where the temperature is controlled with the Nosé-Hoover thermostat [30] and the pressure with the Martyna barostat [31], we performed a series of MD simulations for cavitation under “external force” and evaluated the free energy change based on the Jarzynski formula, eq. (6). Each simulation consists of three steps:

• Liquid in equilibrium at given temperature $T_0$ and pressure $p_0$ is prepared using a standard procedure of MD simulation.
• A spherical force field ${\varphi }_B$ at the center of the system box is applied for all LJ particles, and the system is equilibrated again at $(T_0,p_0)$. We assume an LJ-type field as

  $${\varphi }_B\left(r_o,t=0\right)=4\epsilon \left[{\left({\displaystyle \frac{\sigma }{r_o}}\right)}^{12}-{\left({\displaystyle \frac{\sigma }{r_o}}\right)}^6\right],$$

  (8)

  where $r_o$ is the distance of each particle from the box center. The field is also truncated at $r_o=r_c=3.5\sigma$.
• For the main simulation, the radius of the force field is increased and the resulting “bubble” growth is monitored.

  $${\varphi }_B\left(r_o,t\right)=4\left[{\left({\displaystyle \frac{\sigma }{r_o-R\left(t\right)}}\right)}^{12}-{\left({\displaystyle \frac{\sigma }{r_o-R\left(t\right)}}\right)}^6\right],$$

  (9)

  where $R\left(t\right)$ is a time-dependent variable corresponding to the bubble radius, where the bubble is the region of excluding LJ particles. In this paper we adopted ${\varphi }_B$ with a size linearly increasing with time,

  $$R\left(t\right)=v_{0}t$$

  (10)

  where $v_0$ is a positive constant.

During the MD calculation, atomic data (position ${\overrightarrow{r}}_i$, velocity ${\overrightarrow{v}}_i$, and the exerted external force ${\overrightarrow{f}}_i\equiv -\partial {\varphi }_B/\partial {\overrightarrow{r}}_i$) for each LJ particle i are stored every k steps for later analysis.

To apply the Jarzynski equality to our system for the free energy evaluation, a large number of samplings are required; for that purpose, we took ensemble average over 100 different atomic configurations at equilibrium as the initial conditions.

The simulation conditions are summarized in Table 1. We chose a relatively low temperature ($T=0.8$) since we have to assume vacuum in the “bubble” for this free energy evaluation with applying an external force field. Four different environmental pressure conditions, $-0.30\le p_0\le -0.15$, are investigated to see how $p_0$ affects the nucleation behavior. In view of the standard phase diagram and the equation of states for LJ fluid [32,33,34], these pressure conditions are extremely mild, i.e., much closer to the binodal line than to the spinodal one, as shown in Figure 1. This is the region difficult for conventional MD simulations to access [16,17,18,35] due to the large activation energy. Since a large size of critical nucleus was expected for the largest pressure condition (Case 4), we prepared a larger system with $108,000$ particles to mitigate the artifact of the finite system size. The growth speed $v_0=0.01$ corresponds to $1.6$ m/s in argon units, being still very fast; the speed dependence of results is shortly discussed later. For the larger system, Case 4, we adopted an even larger speed $v_0=0.02$ to save the computation time.

3.2. Free Energy Evaluation

On completing the MD simulations for each Case, we adopted the Jarzynski equality, Eq. (6), to evaluate the free energy change during the “bubble” growth with the following steps:

Step. 1:

The work during a short period from time t and $t+\Delta t$ is evaluated as

$$w\left(t\right)=\sum _{i=1}{\overrightarrow{f}}_i\left(t\right)·{\overrightarrow{v}}_i\left(t\right)\Delta t$$

(11)

In this simulation we chose $\Delta t=5\mathrm{MD}\mathrm{steps}=0.01$ [$\tau$].

Step. 2:

The accumulated work, which is the sum of the instantaneous work up to time $t=n\Delta t$, is determined.

$$W\left(t\right)=\sum _{n^{\prime }=0}^{n}w(n^{\prime }\Delta t)$$

(12)

Step. 3:

We perform simulations multiple times, with calculating $W\left(t\right)$ for each case. The free energy change $\Delta F\left(t\right)$ from the initial state is evaluated by taking the ensemble average as the Jarzynski equality, Eq. (6), as a function of t.

Step. 4:

Since the bubble “radius” is directly related to time t with Eq. (10), we finally obtain the free energy as a function of bubble radius.

4. Results

4.1. Thermodynamic Properties

We carried out NPT ensemble MD simulations for the main sampling, during which the system volume V changes due to the time-varying external field ${\varphi }_B\left(t\right)$. An example is shown in Figure 2 for Case 1, indicating that the barostat to control the system pressure fails after time $t\sim 800$ [$\tau$]. This is caused by a spontaneous bubble growth, as discussed below.

4.2. Bubble Growth

Figure 3 indicates a series of snapshots (cross-sectional views) for Case 1. The “bubble” generated by the external field ${\varphi }_B$ gradually grows up to $t\simeq 700$; it seems to become unstable at $t\sim 800$, where its shape is distorted and the growth is accelerated, and finally the thin liquid film between the neighboring bubbles due to the periodic boundary conditions is broken.

A quantity relevant to the nucleation analysis is the bubble radius at time t, and we evaluate it in two ways here. One is based on the density profile of LJ particles. An example for Case 1 ($p_0=-0.2$) is shown in Figure 4 (a), where the number density is plotted as a function of the distance r from the box center, i.e., the center of the external field ${\varphi }_B$. The first peak indicates the particles situated on the “bubble surface”, from which we define the bubble radius $r_{Bd}$. As the size R of ${\varphi }_B$ increases with time, Eq. (10), the position of the first peak shifts to larger r accordingly. We expect that $r_{Bd}$ is close to the position of ${\varphi }_B$ minimum, $R_{\min }\left(t\right)\equiv R\left(t\right)+2^{1/6}\sigma$; comparison is made in Figure 4 (b), which shows they are almost identical. In the density profile, it is still possible to see the first peak after the bubble becomes unstable at $t\ge 800$ [$\tau$], but the peak is very small and the density in its vicinity is much less than the bulk liquid one, which suggests that the first peak at this stage is brought by weak “adsorption” on, or entrapment by, the external field.

The change of the total volume of the system during the simulation $V\left(t\right)$, an example of which is shown in Figure 2, provides another way to define the bubble radius. We assume that the density of “bulk” liquid is constant, which is apparent in Figure 4 (a). Then the change of $V\left(t\right)$ is brought only by the bubble growth, and we define and evaluate the bubble radius $r_{Bv}$ with spherical bubble assumption.

The obtained $r_{Bv}$ is compared with $r_{Bd}$ in Figure 5. The fluctuations in box size bring the large fluctuations in $r_{Bv}$ especially at the initial stage of $t\le 300$ [$\tau$], but the agreement of $r_{Bd}$ and $r_{Bv}$ is reasonable for $t\le 700$ [$\tau$], as seen in Figure 5 (b). Based on these data, we adopt a rather arbitrarily determined criterion in the following data analyses that the spontaneous bubble growth starts when the difference $r_{Bv}-r_{Bd}$ exceeds $1\sigma$; in this time region, the application of the Jarzynski’s scheme fails because the “work” may not be properly evaluated with Eq. (11).

4.3. Work and Free Energy

Example data of stepwise work $w\left(t\right)$ and cumulative one $W\left(t\right)$ are shown in Figure 6; although fluctuations are large in $w\left(t\right)$, their running sum $W\left(t\right)$ seems to have a reasonable shape as expected in typical nucleation processes. By utilizing the relation between $r_{Bd}$ and t as in Figure 4 (b), the cumulative work W is plotted as a function of bubble radius $r_{Bd}$ in Figure 7, which illustrates the work required to generate a bubble of size $r_{Bd}$. As expected in CNT, $W\left(r_{Bd}\right)$ increases up to some “critical size”; bubbles larger than this size spontaneously grows, leading to negative W.

To evaluate the free energy change, we carried out 100 MD simulations starting from different initial configurations, each of which gives a different $W\left(r_{Bd}\right)$ curve. All curves are superimposed in Figure 8, indicating large fluctuations of $\sim 50$ [$ϵ$]. Ensemble average based on Eq. (6) gives a relatively smooth curve of free energy difference

$$\Delta F\left(r_{Bd}\right)\equiv F\left(r_{Bd}\left(t\right)\right)-F\left(r_{Bd}(t=0)\right)$$

(13)

Note that the applied external field ${\varphi }_B$, Eqs. (9) and (10), has a finite size even at $t=0$, thus $\Delta F\left(r_{Bd}\right)$ does not start from the origin but from $r_{Bd}\simeq 2^{1/6}$ [$\sigma$]. Figure 8 clearly indicates the existence of a critical nucleus at $\simeq 6.5$ [$\sigma$] with activation energy $190\pm 10$ [$ϵ$]. It is apparent that the Jarzynski equality gives the $\Delta F$ curve very close to the data of minimum W.

The results of $\Delta F(r_{Bd}$) are summarized for all four cases of different pressure conditions in Figure 9 (a). We successfully obtained a smooth curve for each condition; as expected, the critical bubble size and the energy barrier reduces as $p_0$ becomes more negative. The critical size and the activation energy are roughly estimated from the data in Figure 9 (b); the dependence on the surrounding pressure $p_0$ seems reasonable.

The evaluated activation energy is very large because we investigated systems under mild $p_0$ conditions. For example, the obtained energy is about 400 [$ϵ$] at $p_0=-0.15$, leading to an extremely small Boltzmann factor $exp\left(-400/0.8\right)\sim {10}^{-216}$, which is not accessible with conventional molecular simulation techniques.

5. Discussion

5.1. Comparison with CNT

Here we only make comparison with the simplest form of CNT. The obtained $\Delta F\left(r_{Bd}\right)$ data are fitted to a simple equation by the least squares method,

$$\Delta F\left(r\right)=a{r}^3+b{r}^2+c$$

(14)

where a, b, and c are the fitting parameters. In a simplest CNT as shown in Eq. (1), a and b are related to the free energy density difference between the bulk liquid and bulk vapor, $\Delta f_{\mathrm{V}}\equiv f_{\mathrm{liq}}-f_{\mathrm{vap}}$, and the surface tension $\gamma$, respectively. As Figure 9 (a) indicates, the obtained free energy curves are well approximated by Eq. (14). Plotted in Figure 10 are the $p_0$ dependence of extracted values,

$$\Delta f_{\mathrm{V}}={\displaystyle \frac{a}{4\pi /3}}$$

(15)

and

$$\gamma ={\displaystyle \frac{b}{4\pi }}$$

(16)

Since the inside of the generated “bubble” is vacuum in this scheme, $\Delta f_{\mathrm{V}}$ should correspond to the free energy density of expanded liquid, and the results seem reasonable that $f_{\mathrm{V}}$ of liquid increases for larger metastability, or negatively larger $p_0$.

The pressure dependence of surface tension is worthy to be noted. As far as we have noticed, no data for surface tension (interfacial tension between liquid and vapor) in non-equilibrium states are reported, although a number of studies exist for its curvature dependence for nano-bubbles [36,37,38,39]. For equilibrium LJ fluid at $T=0.8$. $\gamma \simeq 0.8-0.9$ [$ϵ/{\sigma }^2$] is often referred [33]. Figure 10 suggests that $\gamma$, or the surface excess free energy, is an increasing function of the “degree of non-equilibrium,” which seems reasonable; the evaluation of its accuracy is left for future investigation.

5.2. Sampling

As the main part of free energy evaluation, we used 100 samples for each Case. To see if 100 samples are sufficient or not, we compared the evaluated $\Delta F$ in Figure 11 for Case 1, where the results using the first 20, the first 50, and all 100 samples are plotted. The difference in the averaged values seems negligibly small, and we may conclude that 100 samples are sufficient. However, closer look at the distribution of work may give a different view.

Plotted in Figure 12 are a histogram of instantaneous work $w\left(t\right)$ at different time. Since data with smaller w has larger weight in Jarzynski’s average, Eq. (6), we need sufficient number of samples in the region of smaller w. However we have only four samples at $t=300$ and two at $t=600$ which are smaller than the obtained $\Delta F$, suggesting that 100 samples may not be sufficient. Thus, further increase in the number of samples may significantly lower the free energy.

5.3. Growth Speed

Finally, we briefly discuss how the bubble growth speed $v_0$ affects the evaluation of $\Delta F$. In preliminary investigations, five cases with different $v_0$ ($0.005\le v_0\le 0.04$) are compared, as in Figure 13, where $\Delta F$ is evaluated with Eq. (6) with ten samples for each case. The free energy increases for larger expanding speeds probably because the liquid structure is not sufficiently relaxed for faster bubble expansion, but the results for $v_0=0.01$ and $0.005$ look similar, and thus we have chosen $v_0=0.01$ for the main simulations.

6. Conclusions

Using a standard molecular dynamics (MD) simulation technique, we investigated the free energy change during homogeneous bubble nucleation based on a stochastic thermodynamics formulation, namely, the Jarzynski’s equality, for a simple Lennard-Jones fluid system with an expanding external force field. It is shown that, when compared with the direct application of conventional MD simulations for various types of nucleation processes, we can evaluate the free energy under much milder (i.e., closer to the binodal lines) conditions. Yet assumption of thermal equilibrium is not required for free energy evaluation, in contrast to the often-adopted molecular simulations for thermal properties.

It should be noted, however, sufficient number of samplings is essential, and may require much computational resources. Also exist technical restrictions for this method, that the inside of the bubble should be vacuum due to the external force field, and the bubble shape should be spherical. Therefore application of this method is currently limited only to bubble nucleation at low temperature conditions.