# Investigation of the particle-particle interaction effects in the cosmic
Zevatron based on cyclotron auto-resonance by particle-in-cell simulations

Yousef I. Salamin  Corresponding author<EMAIL_ADDRESS>Department of
Physics, American University of Sharjah, POB 26666, Sharjah, United Arab
Emirates  Qian Zhao<EMAIL_ADDRESS>School of Physics, Xi’an Jiaotong
University, Xi’an 710049, China  Ting Sun<EMAIL_ADDRESS>School
of Physics, Xi’an Jiaotong University, Xi’an 710049, China

###### Abstract

Cyclotron autoresonance acceleration has been recently advanced as a potential
mechanism for accelerating nuclei to ZeV energies (1 ZeV = $10^{21}$ eV). All
results have been based on single- and many-particle calculations employing
analytic solutions to the relativistic equations of motion in the combined
magnetic and radiation fields, excluding effects related to the particle-
particle interactions. Here, results from many-particle calculations and
Particle-In-Cell (PIC) simulations, are presented which lend support to the
single-particle investigations. Each single-particle result is found to lie
well within one standard deviation about the ensemble average obtained from
the corresponding many-particle calculation. The PIC simulations show that,
even for number densities far exceeding those employed in the non-interacting
case, the energy gain drops markedly due to the particle-particle
interactions, over the first $\sim 8~{}mm$ of the acceleration length.
Together with the substantial attenuation, this finding supports the
conclusion that the particle-particle interaction effects can be negligibly
small over acceleration lengths of typically many kilometers.

††preprint: APS/123-QED

##  I Introduction

The binary neutron-star merger, source of the gravitational waves detected
recently [1], was followed by the emission of ultraintense radiation, with
frequencies covering a substantial part of the known electromagnetic spectrum
[2, 3, 4, 5]. Beamed gamma-rays, in the form of a gamma ray burst (GRB),
x-rays, and visible light were emitted. Among other things, these emissions
are a clear indicator of stellar nucleosynthesis and the presence of atoms [6,
7, 8, 9, 10].

Investigation of the emitted radiation can, in principle, be a source of
valuable information about the merger and subsequent evolution of the newly
formed entity. On the other hand, interaction of charged particles with the
beamed radiation, especially in the added presence of superstrong magnetic
fields associated with the merging entities, can drastically influence the
subsequent evolution of such particles. The question thus arises as to whether
atomic nuclei can be accelerated to ZeV energies [11, 12, 13, 14, 15, 16] and
ejected as ultra-high-energy cosmic-rays (UHECR) as a result.

This work is part of efforts dedicated to answering this question [17].
Detection of such particles is quite rare. Only 72 events, with energies
exceeding 57 EeV (1 EeV = $10^{18}$ eV) were detected by the Telescope Array
experiment [18] between 2008 and 2013.

The mechanism of cyclotron autoresonance acceleration (CARA) has recently been
advanced [19] as a possible explanation for the ZeV energies of UHECR
particles. Calculations have demonstrated ZeV energy gains by the nuclei of
hydrogen, helium and iron, due to interaction with ultraintense radiation and
superstrong uniform magnetic fields. The radiation-reaction effects were shown
to be important in CARA but not to lower the energy gain substantially from
the ZeV level.

The investigations in [19, 20] were general in nature and aimed at theoretical
proof-of-principle demonstration of CARA [21, 22, 23, 24, 25, 26, 27, 28, 29]
in an astrophysical context. They did not make specific reference to any known
astrophysical environment where the resonance conditions (and ultraintense
radiation and superstrong magnetic fields) may be found. These conditions may
exist during the brief merger time of two compact objects, over a small area
around the polar cap of one such object, during a magnetar-powered supernova
explosion, among other possibilities [30, 31, 32]. Away from the polar caps,
topology of the steady-state magnetic field of a compact object can be much
more complex than uniform and its lines can be severely curved. The requisite
radiation-field intensity for CARA to work deserves some discussion, too. This
is offered at the very end of Sec. II below.

As such, CARA can be put forth as a mechanism for cosmic-ray acceleration,
alternative to or complementing the widely discussed models based on, for
example, shock waves, magnetic reconnection and unipolar induction [13]. The
existing models describe acceleration to energies close to the EeV level
inside a potential cosmic-ray source, where a plasma background plays a
central role. To reach the ZeV energy levels, it seems plausible to assume
that a particle is first pre-accelerated inside the source by the shock wave
mechanism, for example, and subsequently receives a big energy boost from CARA
outside the source. This assumption will be made throughout this work.

The recent investigations employing CARA [19] have also been single-particle.
Many-particle calculations have been carried out in [20] but did not discuss
the particle-particle interactions. In the present work, many-particle
simulations are carried out to lend support to the single-particle
calculations and to take into account the inter-particle corrections. Thus,
the main working equations of CARA need to be amended slightly. The entities
to be accelerated will be assumed to be initially picked randomly from an
ensemble of $N$ particles. Shape and size of the ensemble will be decided
plausibly, and $N$ will be chosen so that the number of particles per unit
volume is kept way below solid density.

The aim of this work is two-fold: (a) to support the findings of the single-
particle calculations in [19] with many-particle simulations, and (b) to
strengthen the case for CARA by performing simulations which employ a more
realistic set of astrophysical parameters than has been used in [19]. Included
in the latter aim is also presenting, for the first time, results for
acceleration by CARA of an ensemble of nickel nuclei.

Figure 1: CARA: a many-particle schematic diagram.

In Section II, the CARA working equations will be revisited in order to
incorporate the set of initial conditions appropriate for an ensemble of
particles. Dynamics of the ensemble of particles will be investigated, based
on the revised equations, in Section III, employing a parameter set (and for
nuclei) the same as in [19]. In section IV, similar simulations will be
performed for: (a) iron and nickel nuclei, the latter not covered by the
single-particle calculations in [19], (b) a more realistic parameter set, and
(c) smaller ensembles, to ensure that the particle-particle interactions may
be considered negligible. Results from Particle-In-Cell (PIC) simulations,
whose purpose is to shed some light on the particle-particle interaction
effects, will be presented in Sec. V. A brief discussion of our results will
be conducted, and some concluding remarks will be given, in Section VI.

Figure 2: Log-Log plot of the exit kinetic energy of the nuclide Fe+26 with
the radiation-field intensity. For acceleration by CARA using light of
wavelength $\lambda=1~{}\mu$m, the initial injection kinetic energy is
$K_{0}=150$ MeV. For the GBR of wavelength $\lambda=5\times 10^{-11}$ m,
$K_{0}=20$ TeV. In both plots interaction is with 5 radiation-field phase-
cycles. The horizontal dotted lines represent exit kinetic energies of 1 EeV
and 1 ZeV, respectively.

##  II The equations

Figure 1 is a schematic diagram showing the initial ensemble of $N$ identical
particles, each of mass $M$ and charge $Q$, moving along the $z$-axis of a
Cartesian coordinate system. Their positions are uniformly distributed inside
a cylinder of radius $R$ and height $H$. Their initial speeds are derived from
a normal distribution of their injection kinetic energies, of mean
$\bar{K}_{0}$ and standard deviation $\Delta K_{0}$. This makes the number
density $n_{d}=N/(\pi R^{2}H)$. The schematic diagram also shows a uniform
magnetic field of strength $B_{s}$, oriented along $+z$, and a radiation wave
propagating along the same direction. Only the size and shape of the initial
ensemble, and the mean and spread of the initial kinetic energies, will be
fixed. In all many-particle calculations in this work, the initial ensemble
size will be determined by the choices $R=5$ m and $H=\lambda$, the wavelength
of the radiation field employed.

The electromagnetic fields will be modeled by [19]

$\displaystyle\bm{E}$ $\displaystyle=$
$\displaystyle\hat{\bm{i}}E_{0}\sin\eta,$ (1) $\displaystyle\bm{B}$
$\displaystyle=$
$\displaystyle\hat{\bm{j}}\frac{E_{0}}{c}\sin\eta+\hat{\bm{k}}B_{s}.$ (2)

In these equations, $E_{0}$ is the constant amplitude of the plane-wave
radiation, $c$ is the speed of light in vacuum, $\eta=\omega t-kz$,
$k=\omega/c$, and $\hat{\bm{i}},\hat{\bm{j}}$ and $\hat{\bm{k}}$ are unit
vectors in the positive $x$-, $y$\- and $z$-directions, respectively. Recall,
at this point, the resonance condition that characterizes CARA. It ties the
particle and radiation and magnetic field parameters by $r=1$, where

$r=\frac{\omega_{c}}{\omega}\sqrt{\frac{1+\beta_{0}}{1-\beta_{0}}};\quad\omega_{c}=\frac{QB_{s}}{M},$
(3)

in which $\omega$ is the radiation frequency and $\omega_{c}$ is the cyclotron
frequency of the particle around the lines of the magnetic field. Resonance
occurs when the cyclotron frequency matches the Doppler-shifted frequency of
the radiation field, which the particle senses in its own rest frame. With
$Q,M$, and $\omega$ fixed, the resonance condition is essentially a
relationship between $B_{s}$ and $\beta_{0}$, the randomly-selected initial
speed.

Table 1:  Parameters, exit kinetic energies, and resonance magnetic fields for $N=100$ protons accelerated by CARA. Intensity (W/m2) |  $\lambda$ (m) |  $\bar{K}_{e}\pm\Delta K_{e}$ (ZeV) |  $\bar{B}_{s}\pm\Delta B_{s}$ (MT)
---|---|---|---
$10^{38}$ | $10^{-6}$ | $0.437\pm 0.002$ | $11.25\pm 0.03$
$10^{43}$ | $5\times 10^{-11}$ | $2.663\pm 0.028$ | $9.22\pm 0.09$

We proceed now to amend the solutions to the equations of motion of a single
particle in the presence of the electromagnetic fields, by properly
incorporating the above-mentioned initial conditions. The obtained equations
will be used in the next section to carry out the promised many-particle
simulations. On-resonance solutions to the relativistic Newton-Lorentz
equations of motion follow essentially the same steps as in [19]. With the
initial conditions on position expressed as $\eta_{0}=-kz_{0}$, one finally
obtains

$\displaystyle x(\eta)$ $\displaystyle=$ $\displaystyle
x_{0}+\frac{ca_{0}}{2\omega}\gamma_{0}(1+\beta_{0})\left[(\sin\eta-\sin\eta_{0})-(\eta\cos\eta-\eta_{0}\cos\eta_{0})\right],$
(4) $\displaystyle y(\eta)$ $\displaystyle=$ $\displaystyle
y_{0}+\frac{ca_{0}}{2\omega}\gamma_{0}(1+\beta_{0})\left[(\eta\sin\eta-\eta_{0}\sin\eta_{0})+2(\cos\eta-\cos\eta_{0})\right],$
(5) $\displaystyle z(\eta)$ $\displaystyle=$ $\displaystyle
z_{0}+\frac{c}{\omega}\left(\frac{1+\beta_{0}}{1-\beta_{0}}\right)\left\\{\left(\frac{\beta_{0}}{1+\beta_{0}}\right)(\eta-\eta_{0})+\left(\frac{a_{0}^{2}}{24}\right)(\eta-\eta_{0})^{2}(\eta+2\eta_{0})\right.$
(6) $\displaystyle\left.+\frac{a_{0}^{2}}{16}\left[(\eta-2\eta_{0})\cos
2\eta_{0}+\eta\cos 2\eta+(2\eta_{0}^{2}-2\eta\eta_{0}-1)\sin 2\eta_{0}-\sin
2\eta\right]\right\\},$ $\displaystyle\gamma(\eta)$ $\displaystyle=$
$\displaystyle\gamma_{0}+\frac{a_{0}^{2}}{8}\gamma_{0}(1+\beta_{0})\left[(\eta^{2}-\eta_{0}^{2})+(\sin^{2}\eta-\sin^{2}\eta_{0})-(\eta\sin
2\eta-\eta_{0}\sin 2\eta_{0})\right].$ (7)

In these equations, $\gamma_{0}=(1-\beta_{0}^{2})^{-1/2}$ and $a_{0}\equiv
QE_{0}/(M\omega c)$. Note that $a_{0}^{2}$ may be thought of as a
dimensionless radiation-field intensity parameter, whereas the radiation field
intensity in W/m2 is $I=c\epsilon_{0}E_{0}^{2}/2$, where $\epsilon_{0}$ is the
permittivity of free space. Equations (4)-(6) give a parametric representation
of the particle’s trajectory. Equation (7) is the particle’s Lorentz factor
(its energy scaled by $Mc^{2}$).

For the special case of initial position at the origin of coordinates,
evolution of the kinetic energy of a particle with $\eta$ may be written as
[19]

$\displaystyle K(\eta)$ $\displaystyle=$ $\displaystyle
K_{0}+\left[\frac{Q^{2}}{16\pi^{2}\epsilon_{0}Mc^{3}}\right]\gamma_{0}(1+\beta_{0})(I\lambda^{2})$
(8) $\displaystyle\hskip 56.9055pt\times\left[\eta^{2}+\sin^{2}\eta-\eta\sin
2\eta\right],$

where $K_{0}=(\gamma_{0}-1)Mc^{2}$ is the initial injection kinetic energy.
Figure 2 shows log-log plots of the exit kinetic energies against $I$, the
radiation-field intensity, at the end of interaction with 5 phase cycles of
near-visible light and a GRB. The main assumption here is that the particle is
pre-accelerated to kinetic energies of 150 MeV (visible) and 20 TeV (GRB). For
these initial conditions, the resonance magnetic field strengths are 38.9316
MT and 1.09132 GT, respectively. Note that the flat parts of the $K_{e}$ vs.
$I$ curves reflect those initial injection energies. The particle’s kinetic
energy begins to increase substantially after some threshold intensity has
roughly been passed ($\sim 10^{25}$ W/m2, for near-visible light, and $\sim
10^{35}$ W/m2, for the GRB). The energy range of 1 EeV to 1 ZeV is bounded by
the two horizontal dotted lines in Fig. 2. The figure clearly shows that, for
the chosen parameters, the minimum radiation-field intensities required to
reach the EeV to ZeV kinetic energy levels are $\sim 10^{35}$ W/m2 (near-
visible) and $\sim 10^{42}$ W/m2 (GRB).

Figure 3: Proton acceleration by CARA employing near-visible light. (a)
Initial ensemble of $N=100$ protons inside a disk of radius 5 m and thickness
$\lambda=1~{}\mu$m (number density $n_{d}\simeq 1.27324\times 10^{6}$ m-3).
Initial ensemble kinetic energy: normal distribution of mean $\bar{K}_{i}=150$
MeV and standard deviation $\Delta K_{0}=1.5$ MeV. (b) Actual trajectories of
the ensemble members during interaction with the radiation and magnetic
fields. (c) Distribution of the ensemble particles at the end of an
interaction time equivalent to 5 phase cycles ($\Delta\eta=10\pi$) of the
radiation field. (d) Kinetic energy evolution with the excursion distance for
all of the particles in the ensemble. The radiation field intensity is
$I=10^{38}$ W/m2 and the resonance magnetic field strength sensed by the
ensemble members and calculated on the basis of Eq. (3) is $B_{s}=11.2529\pm
0.0354$ MT.

##  III The many-particle calculations

The particle ensemble dynamics in the combined radiation plus uniform magnetic
fields will now be discussed using the equations in terms of the radiation
field phase. Without loss of generality, the examples will focus on the nuclei
H+1, He+2 and Fe+26, as in [19]. Cosmic rays are close to 90% protons, H+1,
the simplest atomic nucleus. Alpha particles, He+2, account for about 10%, and
the rest are heavier nuclei. Fe+26 is one of the most stable nuclei in nature.
Furthermore, recent measurements by the Pierre Auger Observatory [33] suggest
that most UHECR particles are nuclei heavier than the proton.

Figure 4: Acceleration by CARA of protons employing the fields of a GRB. (a)
Initial ensemble of $N=100$ protons inside a disk of radius 5 m and thickness
$\lambda=5\times 10^{-11}$ m (number density $n_{d}\simeq 2.54648\times
10^{10}$ m-3). Initial ensemble kinetic energy: normal distribution of mean
$\bar{K}_{0}=20$ TeV and standard deviation $\Delta K_{0}=0.2$ TeV. (b) Actual
trajectories of the ensemble members during interaction with the radiation and
magnetic fields. (c) Distribution of the ensemble particles at the end of an
interaction time equivalent to 5 phase cycles ($\Delta\eta=10\pi$) of the
radiation field. (d) Kinetic energy evolution with the excursion distance for
all of the particles in the ensemble. The radiation field intensity is
$I=10^{43}$ W/m2 and the resonance magnetic field strength sensed by the
ensemble members and calculated on the basis of Eq. (3) is $B_{s}=9.23921\pm
0.09534$ MT.

As shown schematically in Fig. 1, members of the initial ensemble are assumed
to have been pre-accelerated to relativistic velocities along the common
directions of $\bm{B}_{s}$ and $\bm{k}$, the latter being the wavevector of
the radiation field, by shock waves or any other means [34]. For simplicity,
it will be assumed that the wavefront of a radiation wave catches up with a
particle at $t=0$ when the latter is at the spatial coordinates
($x_{0},y_{0},z_{0}$) and has a speed $\beta_{0}$, derived from the
corresponding initial normal distribution of kinetic energies of mean
$\bar{K}_{0}$ and standard deviation $\Delta K_{0}$. In all cases considered,
the interaction time will be equivalent to five radiation-field phase cycles,
$\Delta\eta=\eta_{e}-\eta_{0}=10\pi$, with $e$ standing for exit.

As examples, we first investigate the dynamics of $N=100$ protons accelerated
by CARA, without reference to any astrophysical environment, known to a good
degree of certainty, where the conditions for acceleration may be met [35, 36,
37, 38, 39, 40, 41, 42, 43]. Figure 3 displays the results for acceleration
using near-visible light of intensity $I=10^{38}$ W/m2 and wavelength
$\lambda=1~{}\mu$m. Speeds of the particles of the initial ensemble are
derived from a normal distribution of kinetic energies ($\bar{K}_{0}=150$ MeV,
and $\Delta K_{0}=1.5$ MeV). Figure 3(a) shows the initial ensemble, a uniform
distribution of initial positions ($x_{0},y_{0},z_{0}$). Interactions are
assumed to commence at $t=0$ (or, equivalently, at $\eta_{0}=-kz_{0}$) causing
the particles to follow the trajectories shown in Fig. 3(b). Figure 3(c) shows
the positions through which the particles pass at the end of the interaction
time. In other words, the initial ensemble in (a) evolves to the final spatial
distribution of particles shown in (c) as a result of the acceleration
process. Assuming that each particle’s own initial conditions launch it into
cyclotron autoresonance, not necessarily exactly, this will result in
tremendous energy gain. In Fig. 3(d) the exit kinetic energy of each particle
of the ensemble is shown as a function of its axial excursion along the
$z$-direction. Exit (end-of-interaction) results for this example are
displayed in the first row of Table 1. Note, in particular, that the magnetic
field strength shown in the last column is given as a mean $\pm$ some spread.
This is due to the fact that once a value for $\beta_{0}$ has been picked at
random, a value for $B_{s}$ will be dictated by the resonance condition, Eq.
(3). Nevertheless, the spread in those values does not seem to disturb
resonance appreciably and the particles end up attaining ZeV kinetic energies.

Table 2: Exit statistics of $10^{4}$ particles accelerated by CARA, employing ultraintense near-visible light.  Nucleus |  $\bar{K}_{e}\pm\Delta K_{e}$ (ZeV) |  $\bar{x}_{e}\pm\Delta x_{e}$ (m) |  $\bar{y}_{e}\pm\Delta y_{e}$ (m) |  $\bar{z}_{e}\pm\Delta z_{e}$ (km) |  $\bar{B}_{s}\pm\Delta B_{s}$ (MT)
---|---|---|---|---|---
H+1 | $0.437\pm 0.002$ | $-0.22\pm 2.48$ | $-0.001\pm 2.516$ | $1352\pm 12$ | $11.25\pm 0.03$
He+2 | $0.334\pm 0.001$ | $-0.05\pm 2.49$ | $-0.029\pm 2.487$ | $197.5\pm 1.4$ | $29.44\pm 0.04$
Fe+26 | $3.281\pm 0.013$ | $-0.06\pm 2.49$ | $0.011\pm 2.503$ | $113.0\pm 0.7$ | $38.93\pm 0.01$

The second illustrative example also involves acceleration of $N=100$ protons,
albeit employing the fields of a GRB of intensity $I=10^{43}$ W/m2 and
wavelength $\lambda=5\times 10^{-11}$ m. In this case, the initial ensemble
kinetic energy (normal) distribution has mean $\bar{K}_{0}=20$ TeV and spread
$\Delta K_{0}=0.2$ TeV [44, 45]. Figure 4 displays results of simulations for
this example similar to those of Fig. 3. Numerical values of the exit
dynamical quantities pertaining to this example are displayed in the second
row of Table 1. The exit kinetic energies of the protons from interaction with
the GRB are substantially larger than from interaction with the lower-
intensity near-visible radiation, as expected. The resonance magnetic field in
this case is lower than in the case of interaction with the lower-frequency
near-visible light, as Eq. (3) predicts.

Further results from simulations performed essentially along the same lines,
albeit involving a much bigger ensemble, two more nuclei, and exhibiting more
exit numerical values, will be presented next. Only the exit numerical values
will be displayed in tabular format, for the nuclei H+1, He+2, and Fe+26.
Table 2 shows results for the acceleration of $N=10^{4}$ particles by CARA,
employing near-visible light of intensity $I=10^{38}$ W/cm2 and wavelength
$\lambda=1~{}\mu$m. Fairly good estimates of the exit mean and spread of the
kinetic energy and spatial coordinates, as well as the resonance magnetic
field strength, may be read from the tabulated results. For example, from the
last row for iron, one concludes that an ensemble of $N=10^{4}$ nuclei
uniformly distributed initially inside a cylinder of radius $R=5$ m and height
$H=\lambda=1~{}\mu$m, evolves into roughly a cylinder bounded by a box of
dimensions $2\Delta x_{e}\sim 5$ m, $2\Delta y_{e}\sim 5$ m, and $2\Delta
z_{e}\sim 1.4$ km.

Table 3: Exit statistics of $10^{4}$ particles accelerated by CARA, employing parameters of a GRB. Nucleus |  $\bar{K}_{e}\pm\Delta K_{e}$ (ZeV) |  $\bar{x}_{e}\pm\Delta x_{e}$ (m) |  $\bar{y}_{e}\pm\Delta y_{e}$ (m) |  $\bar{z}_{e}\pm\Delta z_{e}$ (km) |  $\bar{B}_{s}\pm\Delta B_{s}$ (MT)
---|---|---|---|---|---
H+1 | $2.664\pm 0.029$ | $-0.030\pm 2.507$ | $0.003\pm 2.486$ | $10064\pm 213$ | $9.22\pm 0.09$
He+2 | $0.675\pm 0.007$ | $-0.021\pm 2.501$ | $0.044\pm 2.500$ | $161.6\pm 3.4$ | $72.8\pm 0.7$
Fe+26 | $0.585\pm 0.006$ | $0.024\pm 2.497$ | $0.014\pm 2.509$ | $0.72\pm 0.02$ | $1091\pm 11$

Table 3 is similar to Table 2, but using the fields of a GRB of intensity
$I=10^{43}$ W/m2 and wavelength $\lambda=5\times 10^{-11}~{}$m. The results
presented in Tables 2 and 3 follow different patterns, as they correspond to
two widely differing sets of parameters. Their injection energies, radiation
field intensities, and radiation wavelengths differ by about 5 orders of
magnitude. This leads to different resonance magnetic field strengths. In
Table 3, $K_{e}$ decreases with increasing mass, due to its dependence on the
ratio $Z^{2}/A$, where $Z$ is the atomic number and $A$ is the mass number of
the particle [19]. In Table 2 $K_{e}$ follows the opposite trend, due to
increasing resonance magnetic field values, in jumps of one order of
magnitude. The effect on $K_{e}$ in Table 3 due to $B_{s}$ is not great
because the $B_{s}$ values there are comparable.

##  IV A more realistic parameter space

It may be argued that the parameters employed in our calculations thus far
have been unrealistic. The assumption has been made that a big portion of the
energy output from the source, like a binary-star merger, is radiant and
beamed through a small circle, which leads to the GRB intensities exceeding
$10^{42}$ W/m2, for example, that have been employed in [19]. In the
scientific literature of relevance, the assumption is often made that energy
is radiated isotropically, not in a beam. The intensity calculated based on
this assumption must, therefore, be many orders of magnitude smaller than
$10^{42}$ W/m2.

Another assumption made in our many-particle calculations has been that the
particle-particle interactions are negligible. This has been justified by the
fact that the number densities employed are very small compared to those in a
typical solid, where such interactions can not be ignored. This issue will be
revisited in Sec. V.

Figure 5: Kinetic energy evolution with the excursion distance of an ensemble
of nuclei of Fe+26 (blue) and Ni+28 (red) employing the fields of infrared
light of wavelength $\lambda=0.12$ mm and intensity $I=6\times 10^{32}$ W/m2.
The initial ensemble has $N=10$ nuclei inside a disk of radius 5 m and
thickness $H=\lambda$ (number density $n_{d}\simeq 1061$ m-3). Initial
ensemble kinetic energy: normal distribution of mean $\bar{K}_{0}=1$ GeV and
standard deviation $\Delta K_{0}=10$ MeV. The average resonance magnetic field
strengths sensed by the ensemble members and calculated on the basis of Eq.
(3) are $B_{s}\simeq 287747\pm 202$ T (iron) and $298794\pm 161$ T (nickel).
Interaction time is equivalent to 5 phase cycles ($\Delta\eta=10\pi$).

Our aim thus far has been to lend support to the single-particle results [19]
by performing simulations for non-interacting many particles. Not only do the
many-particle simulations agree, in general, with the single-particle
calculations, but they do not seem to depend on the size of the ensemble
employed.

These points will be addressed together by employing a more realistic
parameter set. Key departure from the old parameter set is employing infrared
radiation of wavelength 0.12 mm and intensity $I=6\times 10^{32}$ W/m2. An
ensemble of only 10 nuclei is employed, making the particle density
$n_{d}\simeq 1061$ m-3 and, therefore, strengthening the argument in support
of the particle-particle interactions being negligibly small. Simulations have
been performed for Fe+26 and Ni+28, two of the most stable nuclei in nature
[46].

Figure 5 shows evolution of the kinetic energies of the particles of ensembles
of nuclei of iron and nickel during interaction with 5 phase cycles of the
radiation field. All other ensemble parameters are given in the figure
caption. The exit kinetic energies are as follows: $K_{e}\simeq 0.3197\pm
0.0014$ ZeV, for iron, and $K_{e}\simeq 0.3317\pm 0.0013$ ZeV, for nickel.
Note that the difference between the two sets of results is quite small. This
is due to the fact that the masses and charge-to-mass ratios of the two nuclei
are quite close. Also, the resonance magnetic field values are comparable.

##  V one-dimensional Particle-in-cell simulations

To approximately asses the effects associated with the particle-particle
interactions in CARA, we consider injection initially of a high-density ion
slab into the interaction region with the radiation and magnetic fields under
near-resonance conditions. The above investigations have shown that the
acceleration length (axial excursion length of the accelerated particles) can
be extremely large, essentially many many kilometers. For that reason alone,
PIC simulations cannot be employed. As will be seen below, dephasing and
spatial spreading will quickly result in a substantial drop in the particle
density. Thus, the particle-particle interaction effects on the exit kinetic
energies of the accelerated particles are expected to be negligibly small.

Despite the above, one-dimensional PIC simulations have been performed using
EPOCH [47] over a small part of the axial excursion length. In each
simulation, a 160 $\mu$m simulation segment, covering the range
$0<z<160~{}\mu$m, moves along the $z-$axis at the speed of light. The segment
is subdivided into $160\times 512$ cells with 64 particles per cell. The
interaction scenario involves He+2 ions initialized into a 1-$\mu$m slab with
a density of $n_{i0}=10^{18}$ cm-3, energy $K_{0}=150\pm 1.5$ MeV, and a
longitudinal static magnetic field $\vec{B}_{0}=B_{0}\hat{z}$, with
$B_{0}=29.65$ MT obtained from the resonance condition. The plane
electromagnetic (PEM) wave is initialized with a wavelength of
$\lambda_{0}=1~{}\mu$m, intensity of $I_{0}=10^{30}$ Wm-2, and a
longitudinally super-Gaussian profile
$f(t)=\exp{[-(t-\tau_{0})^{6}/\tau_{l}^{6}]}$, where $\tau_{0}=55~{}T_{0}$ is
pulse center, and $\tau_{l}=110~{}T_{0}$ is pulse width (with
$T_{0}=\lambda_{0}/c$). In the interaction scenario, the ion slab is first
injected into the region $0<z<1~{}\mu$m of the simulation segment. Then the
PEM is launched behind the slab and overtakes it.

Figure 6: (a)-(c) Snapshots, at propagation time $t=5000~{}T_{0}$, for: (a)
Distribution of the particle number $N$ in phase space (color bar in
logarithmic scale), (b) Distribution of particle energy in the transverse
momentum space, and (c) The local (sensed by the ions in the slab) electric
field component, $e_{x}$, scaled by $E_{0}$ of the PEM wave. Also shown in (c)
is the particle number density $n_{i}$, normalized by the critical density
$n_{c}=1.116\times 10^{21}$ cm-3. All quantities displayed are for
$\alpha-$particles (He+2 ions) in interaction with a PEM wave of 1 $\mu$m
wavelength and a magnetic field of strength $B_{0}=29.65$ MT. (d) Evolution of
the energy gains with axial excursion distance. Here, $\bar{K}$ is the average
kinetic energy of the phase-locked particles inside the particle slab,
$K_{non}$ is corresponding quantity in the non-interacting case based on Eq.
(7), $K_{max}$ is the maximum kinetic energy of the dephased particles.

During the propagation, because of the particle-particle Coulomb interactions
and interaction of the particles with the current-induced magnetic field, in
addition to the energy spread of the initial ensemble, many ions get disturbed
continuously and violate the resonance condition, i.e., they phase-slip or get
dephased. The dephased ions are thus knocked out of the slab and slip back
through the flat top and down the ramp of the PEM pulse. Consequently, they
form the belt distribution in phase space shown in Fig. 6(a). As may be
inferred from the phase-space trajectory oscillations, slipping back through
the flat top causes the dephased ions to gain higher energies, while going
down the PEM ramp results in energy loss. From the energy distribution in the
transverse momentum space shown in 6(b) one can see that the rise in energy
gain is associated with increasing transverse momentum. Note also that the
maximum transverse momentum is a constant, determined by the specific set of
values taken by the parameters $I_{0}$, $B_{0}$ and $K_{0}$.

Propagation of the ion slab is illustrated in Fig. 6(c). It is shown that the
front-edge phase of the PEM catches up with the slab and co-propagates with
it. Density of the ion slab is increased by an order-of-magnitude compared to
the initial density due to the compression caused by the focusing current-
induced magnetic field. On the other hand, the slab is shortened due to the
loss of dephased ions.

The gain curves, showing the maximum kinetic energy $K_{max}$ of the dephased
ions and the averaged kinetic energy $\bar{K}$ of the phase-locked ions within
the slab, are presented in Fig. 6(d), together with the average energy gain of
the non-interacting particles $K_{non}$. Note that $\bar{K}$ of the ion slab
behaves roughly the same as $K_{non}$ obtained from the many- and single-
particle calculations based on Eq. (7) with constant $a_{0}=1.5$,
corresponding to $a_{x}\equiv Qe_{x}/(M\omega c)\simeq 1.5$, where $e_{x}$ is
the $x-$component of the local electric field sensed by the ions in the slab.
The average kinetic energy $\bar{K}$ is in general lower than $K_{non}$ due to
the effect of up-ramp of the PEM pulse. It seems from Fig. 6(d) that, after
propagation for 8 mm, the particle-particle interaction effects lower the
energy gain from about 6.3 to 5 TeV (or by roughly 20%) and the particle
density by more than an order-of-magnitude.

##  VI Concluding remarks

This work is part of efforts to lend support to conclusions arrived at
recently regarding the acceleration to ZeV energies, of protons and other bare
atomic nuclei, by cyclotron autoresonance, in astrophysical environments such
as the merger of a binary neutron-star system or a magnetar-powered supernova
explosion. Those conclusions were based on single-particle calculations [19],
as well as non-interacting many-particle simulations [20] in which the number
density is kept well below that of a typical solid. The current study has
advanced the calculations to the level of interacting many particles,
employing a particle number density close to that of a solid. The process is
assumed to take place outside a compact object, or equivalent, so that the
particle-particle interactions may be ignored. That the particle-particle
interaction effects can be ignored has been supported by PIC simulations
employing, initially, a thin high-density ion slab. The effects can be
strongest between the ions initially. However, interaction with the combined
radiation and magnetic fields, as well as the inter-particle Coulomb forces
and interaction with the ion-generated currents, change that quickly. Many
ions violate the resonance condition, get dephased and slip back and to the
sides. This causes the (accelerated) ion density to fall sharply and the
particle-particle separations to grow axially and render the particle-particle
interactions negligibly small. However, these conclusions cannot be easily
generalized to apply to the ions over many kilometers.

Results of the many-particle simulations strongly agree with those of the
single-particle calculations, with the recognition that the former is more
statistically significant. Let $X$ stand for a physical quantity pertaining to
the particles accelerated by CARA. In [19] and in the present study
$X\in\\{K,x,y,z,B_{s}\\}$. Denote by $X_{e}$ and $X^{\prime}_{e}$ the exit
values of $X$ obtained from the single-particle calculations and the many-
particle simulations, respectively. In all cases considered, and for all
quantities investigated, $X_{e}$ is found to lie within the range
$\bar{X^{\prime}}_{e}\pm\Delta X^{\prime}_{e}$, where $\bar{X^{\prime}}_{e}$
and $\Delta X^{\prime}_{e}$ represent the mean and standard deviation of
$X^{\prime}_{e}$, respectively. Finally, inspection of the numbers (for the
protons, in particular) displayed in Table 1, on the one hand, and 2 and 3, on
the other, reveals that the end results do not depend strongly on the size of
the ensemble.

###### Acknowledgements.

This work is supported in part by the National Natural Science Foundation of
China (Grant number 12105217). The work of YIS has been funded by an Alexander
von Humboldt Fellowship (Wiederaufnahme) and a Faculty Research Grant (FRG
number AS1811) from the American University of Sharjah.

## References

  * Abbott  _et al._ [2017] B. P. Abbott _et al._ , Phys. Rev. Lett. 119, 161101 (2017).
  * Drout  _et al._ [2017] M. R. Drout _et al._ , Science 358, 1570 (2017).
  * Hallinan and others. [2017] G. Hallinan and others., Science 358, 1579 (2017).
  * Kilpatrick  _et al._ [2017] C. D. Kilpatrick _et al._ , Science 358, 1583 (2017).
  * Wu and MacFadyen [2019] Y. Wu and A. MacFadyen, The Astrophysical Journal 880, L23 (2019).
  * Pian  _et al._ [2017] E. Pian _et al._ , Nature 551, 67 (2017).
  * Arcavi  _et al._ [2017] I. Arcavi _et al._ , Nature 551, 64 (2017).
  * Smartt  _et al._ [2017] S. J. Smartt _et al._ , Nature 551, 75 (2017).
  * Kajino  _et al._ [2019] T. Kajino, W. Aoki, A. Balantekin, R. Diehl, M. Famiano, and G. Mathews, Progress in Particle and Nuclear Physics 107, 109 (2019).
  * Wang  _et al._ [2020] X. Wang _et al._ , The Astrophysical Journal 893, 92 (2020).
  * Nagano and Watson [2000] M. Nagano and A. A. Watson, Rev. Mod. Phys. 72, 689 (2000).
  * Harari [2014] D. Harari, Physics of the Dark Universe 4, 23 (2014).
  * Aloisio [2017] R. Aloisio, Progress of Theoretical and Experimental Physics 2017 (2017).
  * Anchordoqui [2018] L. A. Anchordoqui, Phys. Rev. D 97, 063010 (2018).
  * Aab  _et al._ [2020a] A. Aab _et al._ (The Pierre Auger Collaboration), Phys. Rev. Lett. 125, 121106 (2020a).
  * Aab  _et al._ [2020b] A. Aab _et al._ (The Pierre Auger Collaboration), Phys. Rev. D 102, 062005 (2020b).
  * Bell [2013] A. R. Bell, Astroparticle Physics 43, 56–70 (2013).
  * Abbasi  _et al._ [2014] R. U. Abbasi _et al._ , The Astrophysical Journal Letters 790, L21 (2014).
  * Salamin  _et al._ [2021] Y. I. Salamin, M. Wen, and C. H. Keitel, The Astrophysical Journal 907, 24 (2021).
  * Salamin [2021] Y. I. Salamin, Physics Letters A 397, 127275 (2021).
  * Kolomenskii and Lebedev [1963] A. A. Kolomenskii and A. N. Lebedev, Soviet Physics Doklady 7, 745 (1963).
  * Hafizi  _et al._ [1994] B. Hafizi, P. Sprangle, and J. L. Hirshfield, Phys. Rev. E 50, 3077 (1994).
  * Wang and Hirshfield [1995] C. Wang and J. L. Hirshfield, Phys. Rev. E 51, 2456 (1995).
  * LaPointe  _et al._ [1996] M. A. LaPointe, R. B. Yoder, C. Wang, A. K. Ganguly, and J. L. Hirshfield, Phys. Rev. Lett. 76, 2718 (1996).
  * Salamin and Faisal [1998] Y. I. Salamin and F. H. M. Faisal, Phys. Rev. A58, 3221 (1998).
  * Faisal and Salamin [1999] F. H. M. Faisal and Y. I. Salamin, Phys. Rev. A 60, 2505 (1999).
  * Salamin  _et al._ [2000] Y. I. Salamin, F. H. M. Faisal, and C. H. Keitel, Phys. Rev. A 62, 053809 (2000).
  * Galow  _et al._ [2013] B. J. Galow, J.-X. Li, Y. I. Salamin, Z. Harman, and C. H. Keitel, Phys. Rev. ST Accel. Beams 16, 081302 (2013).
  * Salamin  _et al._ [2015] Y. I. Salamin, J.-X. Li, B. J. Galow, and C. H. Keitel, Optics Express 23, 17560 (2015).
  * Vink [2008] J. Vink, Advances in Space Research 41, 503 (2008).
  * Kaspi and Beloborodov [2017] V. M. Kaspi and A. M. Beloborodov, Annual Review of Astronomy and Astrophysics 55, 261 (2017).
  * Soker and Gilkis [2017] N. Soker and A. Gilkis, The Astrophysical Journal 851, 95 (2017).
  * Aab  _et al._ [2017] A. Aab _et al._ (Pierre Auger Collaboration), Phys. Rev. D 96, 122003 (2017).
  * Schlickeiser and Oppotsch [2017] R. Schlickeiser and J. Oppotsch, The Astrophysical Journal 850, 160 (2017).
  * Hillas [1984] A. M. Hillas, Ann. Rev. Astron. Astrophys. 22, 425 (1984).
  * Drury  _et al._ [1994] L. O. Drury, F. Aharonian, and H. J. Völk, Astronomy & Astrophysics 287, 959 (1994).
  * Chen  _et al._ [2002] P. Chen, T. Tajima, and Y. Takahashi, Phys. Rev. Lett. 89, 161101 (2002).
  * Aharonian  _et al._ [2004] F. A. Aharonian _et al._ , Nature 432, 75 (2004).
  * Honda and Honda [2004] M. Honda and Y. S. Honda, The Astrophysical Journal 617, L37 (2004).
  * Drury [2009] L. O. Drury, Plasma Physics and Controlled Fusion 51, 124005 (2009).
  * Osmanov  _et al._ [2014] Z. Osmanov, S. Mahajan, G. Machabeli, and N. Chkheidze, Monthly Notices of the Royal Astronomical Society 445, 4155 (2014).
  * Fang and Murase [2018] K. Fang and K. Murase (Pierre Auger), Nature Physics 14, 396–398 (2018).
  * Liu  _et al._ [2017] R.-Y. Liu, F. M. Rieger, and F. A. Aharonian, The Astrophysical Journal 842, 39 (2017).
  * Kann  _et al._ [2019] D. A. Kann _et al._ , Astronomy & Astrophysics 624, A143 (2019).
  * Lipari and Vernetto [2020] P. Lipari and S. Vernetto, Astroparticle Physics 120, 102441 (2020).
  * Adriani  _et al._ [2021] O. Adriani _et al._ (CALET Collaboration), Phys. Rev. Lett. 126, 241101 (2021).
  * Arber  _et al._ [2015] T. D. Arber _et al._ , Plasma Physics and Controlled Fusion 57, 113001 (2015).


# Thin film Rupture from the Atomic Scale

Muhammad Rizwanur Rahman1<EMAIL_ADDRESS>Li Shen1 James P. Ewen1
Benjamin Collard1,2 D. M. Heyes1 Daniele Dini1 E. R. Smith3 1Department of
Mechanical Engineering, Imperial College London, South Kensington Campus,
London SW7 2AZ, United Kingdom, 2Department of Materials Science, Imperial
College London, South Kensington Campus, London SW7 2AZ, United Kingdom,
3Department of Mechanical and Aerospace Engineering, Brunel University London,
Uxbridge UB8 3PH, United Kingdom

###### Abstract

The retraction of thin films, as described by the Taylor-Culick (TC) theory,
is subject to widespread debate, particularly for films at the nanoscale. We
use non-equilibrium molecular dynamics simulations to explore the validity of
the assumptions used in continuum models, by tracking the evolution of holes
in a film. By deriving a new mathematical form for the surface shape and
considering a locally varying surface tension at the front of the retracting
film, we reconcile the original theory with our simulation data to recover a
corrected TC speed valid at the nanoscale.

Thin films, rupture, Taylor-Culick speed, film retraction, local surface
tension.

††preprint: APS/123-QED

Introduction: Thin film rupture is of premier interest in numerous physical
and chemical transport processes (Thiele  _et al._ , 2001; Masson and Green,
2002; De Gennes  _et al._ , 2004; de Gennes, 2008; Bird  _et al._ , 2010;
Bico, 2015; Petit  _et al._ , 2015; Poulain  _et al._ , 2018), ranging from
disease transmission (Bourouiba, 2021) to planetary scale environmental and
oceanic science (Wu, 1981; Veron, 2015). As such, the film rupture process has
been investigated for over a century (Marangoni  _et al._ , 1872; Ranz, 1959;
Chatzigiannakis and Vermant, 2020); and remains an active area of debate.
Taylor (Taylor, 1959) and Culick (Culick, 1960) independently derived an
equation for film retraction by considering that the growth of a nucleated
hole of radius, $R$, is driven purely by the balance between the surface
force, $F_{\gamma}$, and the inertia of the liquid mass, $m$, collected in a
rim surrounding the hole, i.e., $F_{\gamma}=\frac{d}{dt}(m\dot{R})$. Applying
the initial condition at (time) $t=0$, $R=0$ for finite $\dot{R}$, one
obtains, upon integration, the Taylor-Culick (TC) speed,

$U_{\mathrm{TC}}=\sqrt{{\phi\gamma}/{\rho h_{0}}}$ (1)

where $h_{0}$ is the thickness of the unperturbed part of the film with
density $\rho$ and surface tension $\gamma$. In Taylor (Taylor, 1959) and
Culick (Culick, 1960), $\phi=2$, but later work use smaller values for thinner
films (McEntee and Mysels, 1969). This surprisingly simple, yet elegant
equation captures the overall retraction dynamics, and corrects Dupre’s
inappropriate energy balance assumption (Dupré, 1869). However, its accuracy
is debatable, and it markedly fails for films thinner than $100\,\mbox{nm}$
(McEntee and Mysels, 1969; Evers  _et al._ , 1996, 1997; Planchette  _et al._
, 2019; Pierson  _et al._ , 2020) \- this region constitutes the key focus of
this letter.

Taylor’s resolution is based on a few key assumptions, (i) when rupture
propagates, $h_{0}$ remains uniform except at the rim, and (ii) all liquid-
mass of the expanding hole is collected inside the rim, i.e., $m=\rho\pi
R^{2}h_{0}$. In a more formal derivation, Culick made the assumptions of
uniform thickness and constant surface tension everywhere. Whilst the
assumptions leading to Eq. (1) are not incorrect at the global scale, their
local deviations become increasingly important as thickness decreases down to
nanometers.

Significant deviations from $U_{\mathrm{TC}}$ have been documented, especially
in the limit of very thin films (McEntee and Mysels, 1969; Planchette  _et
al._ , 2019), with a range of continuum based models attempting to explain
these differences by consideration of the non-uniformity of film thickness
(Keller, 1983), surface elasticity (Petit  _et al._ , 2015), surface tension
gradients (De Corato  _et al._ , 2022) and viscous effects (Pierson  _et al._
, 2020; Savva and Bush, 2009). With continuum models, any additional physics
must be included by adding another conservation or balance law. In this
letter, molecular dynamics (MD) simulation is applied to fully understand the
retraction process from the atomic scale. Any rupture event must start at the
smallest scale, so MD is an ideal technique to fully understand the origins of
this process. It provides the full picture of a fluid down to the
electrostatic atomic interactions; and the interface dynamics, surface tension
and viscosity emerge as outputs from MD simulations. Hence, many fluid
properties have been successfully studied by MD, a previous letter in this
journal, Koplik and Banavar (2000) used MD to show the origins of dewetting,
alongside the studies of film breakup (Koplik and Banavar, 1993; Zhao  _et
al._ , 2019), boiling and nucleation (Diemand  _et al._ , 2014) and a range of
cases (Kadau  _et al._ , 2010; Todd and Daivis, 2017). Kono  _et al._ (Kono
_et al._ , 2014) studied a two-dimensional film breakup with MD, and observed
temperature rise due to rapid rupture which would cause a decrease in surface
tension.

Figure 1: (a) Computational domain considered for this study. A hole is
nucleated at the center of the film in the $yz$ plane. A cylindrical needle
schematically shows the hole-poking process. Zoomed view of the hole: atoms in
red denote the atoms initially in the hole at $t=0$ that are displaced; faint
schematic shows an idealization of the rim surrounding the hole which borders
the unperturbed film. (b) Radial average of a retracting film over time:
$R(t)$ is measured from the center of the hole to the tip of the rim. (c)
Simplified schematic of panel (b), film thickness, $h_{0}$, is measured at the
unperturbed part of the film.

In this work, we model the full three-dimensional film breakdown and
systematically outline all possible contributions from the atomic scale. Using
novel local pressure measurement techniques, interface tracking and dynamic
time-evolving mappings, local surface tension is calculated. This uncovers the
full picture of surface force and interface shape, putting the thermal effects
in context and exposing other key factors contributing to the dynamics. The TC
model does not take into account any surface-species and their contributions
to the surface properties of the film (Tammaro  _et al._ , 2021; Manikantan
and Squires, 2020). By construction, the film investigated in this work is
pure in nature, and hence, closest to the TC model. Present work thereby tests
this model by removing external factors and shows that only when certain
atomic-scale corrections are made, TC equation can successfully predict the
retraction process.

Simulation setup: Films investigated in this study were modeled using a
Lennard-Jones (LJ) fluid simulated in the extensively validated and verified
Flowmol MD code (Smith, 2013). The initial simulation domain was a cubic box
of dimensions $L_{x}=76.19$ and $L_{y}=L_{z}=609.56$ in LJ units. These
dimensions correspond to $\sim 207$ nm wide films with $h_{0}<5$ nm. Rest of
this letter reports all quantities in LJ units. The middle $20\%$ of the
simulation box in the $x$ direction was designated as liquid and setup with a
density, $\rho\approx 0.7$ and the remaining as vapor with a lower density,
$\rho=0.01$. Schematic in Fig. 1 (a) shows the system modelled in this study
where all atoms ($\sim 3.54$ million) are identical, forming a central liquid
film coexisting with surrounding vapor at equilibrium. After equilibrating the
film at $T=0.78\pm 0.03$, a hole of initial radius, $R_{0}$ was induced in the
otherwise stable film, and $R$ was measured over time. The inset in Fig. 1 (a)
(also SI Fig. 1 (mrr, )) shows the liquid atoms surrounding an expanding hole.
Displaced atoms (in red) from the hole collect around the hole to form a rim.
The radial-averaged film extracted from MD data, in panel (b), shows its
temporal evolution - with simplified schematic in panel (c).

The hole creation technique may invoke slight differences to the initial
conditions, but any difference rapidly ceases (Frankel and Mysels, 1969). We
confirm this by applying two separate techniques to create holes, by (i)
deleting atoms from a circular region of radius, $R_{0}$, or (ii) applying a
small radially decaying force for a short time that models poking a film.
Henceforth, these two types will be referred to as cut- and poked-holes,
respectively. The growth of spontaneously nucleated holes (by allowing holes
to form naturally due to thermal fluctuations) were found to give similar
dynamics, albeit spontaneous holes take considerably longer time to nucleate.
Each reported case in this letter was averaged over at least three separate
ensembles. The short-time growth, being more susceptible to thermal noise, was
averaged over ten separate ensembles, further details in  S1 (mrr, ).

Surface deformations arising from thermal motion increases the likelihood of
spontaneous rupture (Vrij, 1966; Vrij and Overbeek, 1968; Ivanov  _et al._ ,
1970). Therefore, to find a stable, and computationally tractable case, films
with a range of thicknesses were tested without applying any external
perturbation. A thickness, $h_{0}=14$ was selected as films thinner than this
were observed to break spontaneously before reaching a constant retraction
rate. Other computationally affordable thicknesses were also examined (S1.1
(mrr, )).

Figure 2: (a) Short-term nonlinear retraction of film for cut- and poked-holes
with initial thickness is $h_{0}\sim 14$, time is non-dimensionalized by the
inertial time scale, i.e., $t^{*}=t/\sqrt{\rho h_{0}^{3}/2\gamma}$. Faint and
solid symbols denote, respectively, data from individual simulation, and data
averaged over 10 simulations. Solid line is an exponential fit to the data of
the form: $R(t)-R_{0}=e^{\alpha t+\beta}$. (b) Retraction velocity (blue:
poked-, red: cut-holes) normalized by the TC speed, $U^{*}=U/U_{\mathrm{TC}}$
with $2\%$ standard deviation. Faint symbols denote cases where the classical
$U_{\mathrm{TC}}$ is used, whereas, solid symbols denote cases where
$U_{\mathrm{TC}}$ is corrected considering the local surface tension
($U^{*}_{\gamma}$), and the existence of a transition regime
($U^{*}_{\gamma,\alpha}$). Pie-chart in the inset shows various corrections
arising from the change of $\gamma$ due to temperature, $\epsilon_{\gamma,T}$
and curvature, $\epsilon_{\gamma,\kappa}$, and due to the transition regime,
$\epsilon_{\alpha}$. In panels a and b, the non-linear growth region is
shaded. (c) Schematic in the top-panel shows an idealized retracting film. The
mid-panel sheds molecular insights where a transition region is present
between the rim and the unperturbed film. Green arrows show the directions
normal ($n$) to the local film surface. Color map shows radial-averaged
instantaneous temperature field with velocity vectors (black arrows) overlaid.
Bottom panel shows the local to global (time-averaged) surface tension ratio,
($\gamma^{*}=\gamma_{r}/\gamma$) against distance, $r$ from the tip of the
retracting film.

Results and discussions: We chart the hole-growth as predicted by MD and
compare to experiments, CFD and TC theory, starting from initial hole
formation, through a short-time exponential regime before eventually reaching
a constant speed which can be shown to match TC when we account for local
changes in shape and in surface tension.

Initial hole formation: For systems at macroscopic length scale, growth of
holes with an initial radius, $R_{0}<h_{0}$ is not energetically favoured,
i.e., an open hole closes (Taylor and Michael, 1973). Hence, earlier
experiments or numerical investigations (McEntee and Mysels, 1969; Pandit and
Davidson, 1990; Savva and Bush, 2009; Sanjay  _et al._ , 2022) on sheet
retraction considered an initial hole dimension satisfying, $R_{0}\gg h_{0}$.
Additionally, continuum assumptions require long-wave approximations to
describe the rupture process which is valid for $h_{0}$ much larger than the
atomic scale (Vaynblat  _et al._ , 2001). However, in the limit of a few
nanometers thickness (essentially Newtonian black films), the disjoining force
acts as a destabilizer opposing the stabilizing effect of surface tension, and
therefore, becomes responsible for the puncture of the film (Neél and
Villermaux, 2018; Villermaux, 2020). At this scale; MD is the only appropriate
tool to scrutinize the hole opening process. It was observed that $R_{0}\gg
h_{0}$ no longer remains a requirement for holes to expand ( S1.2 (mrr, )).
This agrees with the seminal study of Koplik and Banavar (2000), where the
authors reported that for any dry patch larger than a few molecular diameters,
the separated molecules suffer from mutual interactions, and the hole always
grows.

Short-term retraction: At the onset of retraction, an exponential growth for a
short period of time is observed, which converges to a linear terminal speed
at later times. Fig. 2 (a) portrays this short-term hole expansion for cut-
and poked-holes. In agreement with the observations in (Debrégeas  _et al._ ,
1995; Masson and Green, 2002), it is evident that the early-time hole growth
shows an exponential behavior, i.e., $R(t)-R_{0}=\mathrm{exp}\,(ct+d)$, where
$c$ and $d$ are fitting parameters. While Debrégeas  _et al._ (Debrégeas  _et
al._ , 1995) ascribed this short-term growth to the viscoelastic effects, it
was confirmed to be a generic feature of circular hole retraction (Savva and
Bush, 2009), and was also accredited to the initial shape of the puncture
(Roth  _et al._ , 2005; Savva and Bush, 2009).

Taylor-Culick retraction: Beyond the short-term exponential regime, the
retraction velocity approaches a plateau - the constant TC speed (Savva and
Bush, 2009). This is shown by the faint symbols in Fig. 2 (b) for cut- (red)
and poked-holes (blue). Retraction speed, $U$ is measured over $50$
consecutive time-steps, and is non-dimensionalized by $U_{\mathrm{TC}}$ from
Eq. (1) with $\phi=2$, i.e., $U^{*}=U/U_{\mathrm{TC}}$. Evident from Fig. 2
(b), $U_{\mathrm{TC}}$ is appreciably faster than the measured speed
($U^{*}\sim 0.8$), that is MD predicts a speed which is only 80% of the
Taylor-Culick prediction. This is not unexpected, and documented in the
literature (McEntee and Mysels, 1969; Müller  _et al._ , 2007; Trittel  _et
al._ , 2013; Petit  _et al._ , 2015; Tammaro  _et al._ , 2021). Notably,
McEntee and Mysels (1969) studied the retraction of soap-films with a wide
range of thicknesses. In their experiments, the retraction speeds for films of
thickness smaller than $100\,\mbox{nm}$ were found to be substantially slower
than $U_{\mathrm{TC}}$, requiring $\phi\sim 0.7$. $U_{\mathrm{TC}}$ is rather
closer to our measured speed requiring $\phi\sim 1.35\pm 0.05$ for $h_{0}\sim
4.75$ nm. However, having the effects of surfactants and visco-elasticity
(Vrij and Overbeek, 1968) completely eliminated in the present study (by
considering pure film), why the TC theory still fails to estimate the
retraction rate is explored now.

A number of factors, i.e., temperature variation due to internal heat
generation, curvature effects and viscous contribution (Pierson  _et al._ ,
2020) might play a critical role in slowing down the retraction. These are
often characterized by the Ohnesorge number,
$\mathrm{Oh}\,=\mu/\sqrt{\rho\sigma h_{0}}$, which is $\sim 0.8$ for the
present case.

Surface tension - temperature effects: Adjusting the value of $\phi$ marries
the theoretical estimate with the measured speed, but it does not necessarily
explain the underlying physics responsible for the slower retraction. This
disagreement was attributed to the disturbance caused by the presence of a
precursor ahead of the rim arising from surface tension gradient due to
surface shrinkage and surfactant concentration (McEntee and Mysels, 1969;
Mysels and Stikeleeather, 1971). Frankel and Mysels (1969) inferred that the
slow growth-rate should be due to a large variation of the local surface
tension at the film-front. Although such local fluctuation is not unusual
(Moody and Attard, 2001; Wen  _et al._ , 2021), capturing these with continuum
models is difficult, and the only MD study is the quasi-2D case of Kono  _et
al._ (2014) who reported considerable temperature rise ($>5K$) in the liquid
surrounding the hole for $h_{0}\sim 5$ nm.

In contrast to 2D systems, fully 3D simulations are computationally more
expensive, but avoid the well-known deficiencies of 2D fluids, such as
divergent transport coefficients Irving and Kirkwood (1950). One expects the
behavior of an interface to be dominated by the number of interacting
neighboring molecules - which is much smaller in 2D. Therefore, a 2D interface
would be overly susceptible to rupture Koplik and Banavar (1993). Even with
periodic boundary conditions, the strong system-size dependency of the
diffusion coefficient cannot also be ignored (Yeh and Hummer, 2004). Present
study therefore minimizes these effects by simulating a fully three-
dimensional (3D) system which is $\sim 8$ times larger (in terms of particle
number) than the 2D system of (Kono  _et al._ , 2014). A warmer region
surrounding the expanding hole is observed, and the temperature increase is
rather less pronounced, $2.9$K vs. $>5$K,  see S2.2 (mrr, ). Such an increase
in temperature constitutes around $5-7\%$ drop in $\gamma$. This is, of
course, non-negligible; but is not sufficient to explain the slower
retraction.

Surface tension - curvature effects: The retraction process is mainly driven
by the tangential component of the pressure tensor, while the normal component
is responsible for the radial growth of the rim. Resolving the pressure tensor
spatially along the film profile in the local tangential and normal directions
( S2.2 (mrr, )) reveals a noticeable depression of the tangential pressure.
From these constituent components, $\gamma$ can be locally calculated using
the Kirkwood-Buff definition (Kirkwood and Buff, 1949), i.e.,
$\gamma_{r}=\frac{1}{2}\int_{-\infty}^{\infty}\left[P_{xx}(r,x)-P_{rr}(r,x)\right]dx$,
where $r$ denotes the radial distance from the film tip, subscripts $xx$ and
$rr$ denote direct pressure along the x and r axis, respectively. A full
rotation and mapping of the pressure tensor using a radially-varying normal
vector does not give a substantially different surface tension (see  S2.2
(mrr, )). The tension acts as a ‘pulling force’ in Taylor-Culick formulation,
purely in the direction of retraction, seen from the Kirkwood-Buff with
$P_{rr}$ the main contribution to $\gamma$. The bottom-panel in Fig. 2 (c)
shows the ratio of the local to global surface tension with increasing $r$.
Remarkably, this gives evidence of a substantial drop in $\gamma_{r}$ ($\sim
15\%$) at the tip of the retracting film ($r=0$), and $\gamma_{r}\to\gamma$
far from the rim. While the temperature rise explains a $\sim 5-7\%$ drop
($\epsilon_{\gamma,T}$) in $\gamma$, the remaining can be attributed to
surface curvature and shrinkage, $\epsilon_{\gamma,\kappa}$, so
$\epsilon_{\gamma}=\epsilon_{\gamma,T}+\epsilon_{\gamma,\kappa}$. Accounting
for this local variation, an improvement
$\epsilon_{\gamma}=1-\sqrt{\gamma_{r=0}/\gamma}$ is obtained; and for the
present case, using the surface tension averaged over the film front in Eq.
(1) results in agreement with the measured speed within only a $\sim 8\%$
deviation, shown by the dashed line denoted by $U^{*}_{\gamma}$ in Fig. 2 (b).
These corrections, however, will diminish for thicker films as the variations
of temperature (Kono  _et al._ , 2014) and curvature become less pronounced (
S2.2 (mrr, )).

Transition regime: The top panel in Fig. 2 (c) shows an idealized picture of
film retraction upon which the TC formalism is based on. This simplified
description requires the liquid from the hole to be collected in a spherical
rim bordering on the unperturbed film. Nature will rectify such abrupt
discontinuity through the effects of viscosity (Rayleigh, 1880). Contrasting
Taylor’s assumption of a sudden transformation, present study rather agrees
with experimentally (Lhuissier and Villermaux, 2009a; Petit  _et al._ , 2015)
and numerically (Koplik and Banavar, 2000) observed transition, shown in the
mid-panel of Fig. 2 (c). Previous studies reported an ‘aureole-shaped’ regime
for films of liquid with surfactants, and attributed this to surface tension
gradient, and surface elasticity (De Corato  _et al._ , 2022; Lhuissier and
Villermaux, 2009b; Petit  _et al._ , 2015). Existence of a similar ‘extended
rim’ in the present study - where the film is made of pure liquid - indicates
a gradient of the surface force, but caused by curvature effects (geometry)
and viscosity (Savva and Bush, 2009).

The classical formulation discounts any possibility of a transition regime
between the rim and the undisturbed film. Our results confirm that the rim is
followed by a transition regime (Fig 1 b). This can be explored through a mass
balance argument by assuming an elliptic blob, of aspect ratio $\alpha=b/a$,
connected by a smooth transition region, of length $L$, to the neck of width
$h_{0}$. The length, $L$ is given by,

$L={\pi h_{0}}(\beta^{2}-1)/{2\alpha}$ (2)

where, $\beta=2b/h_{0}$ (see S2.3 (mrr, ) for full derivation). Transition
must take place over the finite length, $L$ (Frankel and Mysels, 1969) \-
which we confirm by comparing $L$ to our simulations (SI Fig. 6(b) (mrr, )),
thus disproving, for thinner films, the assumption of an unperturbed film
bordering on the rim. As Eq. (2) suggests, the transition length becomes more
prominent with thinner films, and disappears when $b\sim h_{0}/2$
corresponding to viscous film retraction (Brenner and Gueyffier, 1999; Savva
and Bush, 2009).

The MD data indicates that the normal pressure difference of the surface is
minimal at the point where the rim meets the transition region. In order to
satisfy the Young-Laplace equation, it follows that the absolute value of the
curvature must also be minimized. Under these assumptions, the atomistic form
of Eq. (1) can be expressed as

$U=\left(1+\alpha^{2}\right)^{-1/4}U_{\mathrm{TC}}.$ (3)

For simplicity, assuming a spherical rim, i.e., $\alpha=1$ yields a surface-
shape correction factor, $\epsilon_{\alpha}=2^{-1/4}\approx 0.84$. When
$\epsilon_{\alpha}$ is accounted for, in addition to $\epsilon_{\gamma}$, the
corrected TC model ($U^{*}_{\gamma,\alpha}$) shows close agreement with the MD
results. This is shown by the solid symbols in Fig. 2 (b), and as evidenced by
the figure, $U^{*}_{\gamma,\alpha}\sim 1$.

There remain a few other aspects which plausibly affect the dynamics of film
retraction. Any mutual dependency of the two correction factors developed in
this study cannot be ruled out. The growth of the blob is non-negligible as
compared to $\dot{R}$. In situations like this, Pierson  _et al._ (2020) and
others (Savva and Bush, 2009; De Corato  _et al._ , 2022) argued that the
velocity within the blob is non-uniform; the velocity vectors overlaid on the
temperature field in Fig. 2 (c) confirm this hypothesis. The effect of this
velocity would contribute to the departure of the blob shape from being
circular and therefore, is assumed to be included in the surface shape
correction, $\epsilon_{\alpha}$. The effect of evaporation, however, in
agreement with previous study (Kono  _et al._ , 2014), was found to be
negligible.

Conclusion: This study examines the Taylor-Culick theory for film retraction
at atomic scale using large-scale MD simulations. Similar growth rate of holes
made by two separate methods, either poking or cutting, confirms that the
dynamics is independent of hole creation mechanism. This then allows a
systematic study, where we repeatably show a short term exponential growth
period followed by a constant velocity regime over an ensemble of cases - both
behaviours observed in the experimental literature. A Lennard-Jones model is
used to remove viscoelastic effects so the various driving contributions to
the growth dynamics can be exposed. The theoretical conjecture of a constant
terminal retraction rate is accurate at this scale, however, we show that
estimating the ‘true’ retraction speed requires corrections for (i) the local
variation of surface tension arising from surface curvature and temperature
rise, and, (ii) the existence of a transition regime between the rim and the
unperturbed film - both of these factors were ignored in the classical
formulation. Surface tension, viscosity and interface shapes are outputs of MD
simulations and we observe how they contribute to the dynamics, and quantify
their impact. While their contributions are marginal for macroscopic films,
these are the key reasons for the failure of the classical theory as film
thickness is reduced below $100$ nm.

This letter provides compelling evidence that the pioneering theory of Taylor
and Culick, which in its original form fails at the atomic scale, can be
corrected to accurately predict the film retraction speed.

###### Acknowledgements.

M.R. R. thanks Shell and the Beit Trust for PhD funding through a Beit
Fellowship for Scientific Research. L.S. thanks the Engineering and Physical
Sciences Research Council (EPSRC) for a Postdoctoral Fellowship
(EP/V005073/1). J.P.E. was supported by the Royal Academy of Engineering
(RAEng) through their Research Fellowships scheme. B.C. was supported by Shell
and the EPSRC via an iCASE PhD studentship (EP/T517690/1). D.D. acknowledges a
Shell/RAEng Research Chair in Complex Engineering Interfaces and the EPSRC for
an Established Career Fellowship (EP/N025954/1).

## References

  * Thiele  _et al._ (2001) U. Thiele, M. G. Velarde, and K. Neuffer, Dewetting: Film rupture by nucleation in the spinodal regime, Physical Review Letters 87, 016104 (2001).
  * Masson and Green (2002) J.-L. Masson and P. F. Green, Hole formation in thin polymer films: A two-stage process, Physical Review Letters 88, 205504 (2002).
  * De Gennes  _et al._ (2004) P.-G. De Gennes, F. Brochard-Wyart, D. Quéré, _et al._ ,  _Capillarity and wetting phenomena: drops, bubbles, pearls, waves_ , Vol. 315 (Springer, 2004).
  * de Gennes (2008) P. G. de Gennes, Bubbles, foams and other fragile objects, in  _A Century of Ideas: Perspectives from Leading Scientists of the 20th Century_, edited by B. G. Sidharth (Springer Netherlands, Dordrecht, 2008) pp. 77–85.
  * Bird  _et al._ (2010) J. C. Bird, R. De Ruiter, L. Courbin, and H. A. Stone, Daughter bubble cascades produced by folding of ruptured thin films, Nature 465, 759 (2010).
  * Bico (2015) J. Bico, Cracks in bursting soap films, Journal of Fluid Mechanics 778, 1 (2015).
  * Petit  _et al._ (2015) P. Petit, M. Le Merrer, and A.-L. Biance, Holes and cracks in rigid foam films, Journal of Fluid Mechanics 774 (2015).
  * Poulain  _et al._ (2018) S. Poulain, E. Villermaux, and L. Bourouiba, Ageing and burst of surface bubbles, Journal of Fluid Mechanics 851, 636 (2018).
  * Bourouiba (2021) L. Bourouiba, The fluid dynamics of disease transmission, Annual Review of Fluid Mechanics 53, 473 (2021).
  * Wu (1981) J. Wu, Evidence of sea spray produced by bursting bubbles, Science 212, 324 (1981).
  * Veron (2015) F. Veron, Ocean spray, Annual Review of Fluid Mechanics 47, 507 (2015).
  * Marangoni  _et al._ (1872) C. Marangoni, P. Stefanelli, and R. Liceo, Monografia sulle bolle liquide, Il Nuovo Cimento (1869-1876) 7, 301 (1872).
  * Ranz (1959) W. E. Ranz, Some experiments on the dynamics of liquid films, Journal of Applied Physics 30, 1950 (1959).
  * Chatzigiannakis and Vermant (2020) E. Chatzigiannakis and J. Vermant, Breakup of thin liquid films: From stochastic to deterministic, Physical Review Letters 125, 158001 (2020).
  * Taylor (1959) G. I. Taylor, The dynamics of thin sheets of fluid. $iii$. disintegration of fluid sheets, Proceedings of the Royal Society of London, Series A 253, 313 (1959).
  * Culick (1960) F. E. Culick, Comments on a ruptured soap film, Journal of Applied Physics 31, 1128 (1960).
  * Dupré (1869) A. Dupré,  _Théorie mécanique de la chaleur_ (Gauthier-Villars, 1869).
  * McEntee and Mysels (1969) W. R. McEntee and K. J. Mysels, Bursting of soap films. I. An experimental study, The Journal of Physical Chemistry 73, 3018 (1969).
  * Evers  _et al._ (1996) L. J. Evers, S. Y. Shulepov, and G. Frens, Rupture of thin liquid films from Newtonian and viscoelastic liquids Bursting behaviour of Newton-black films, Faraday Discussions 104, 335 (1996).
  * Evers  _et al._ (1997) L. J. Evers, S. Y. Shulepov, and G. Frens, Bursting Dynamics of Thin Free Liquid Films from Newtonian and Viscoelastic Solutions, Physical Review Letters 79, 4850 (1997).
  * Planchette  _et al._ (2019) C. Planchette, F. Marangon, W.-K. Hsiao, and G. Brenn, Breakup of asymmetric liquid ligaments, Physical Review Fluids 4, 124004 (2019).
  * Pierson  _et al._ (2020) J.-L. Pierson, J. Magnaudet, E. J. Soares, and S. Popinet, Revisiting the Taylor-Culick approximation: Retraction of an axisymmetric filament, Physical Review Fluids 5, 073602 (2020).
  * Keller (1983) J. B. Keller, Breaking of liquid films and threads, Physics of Fluids 26, 3451 (1983).
  * De Corato  _et al._ (2022) M. De Corato, D. Tammaro, P. L. Maffettone, and N. Fueyo, Retraction of thin films coated by insoluble surfactants, Journal of Fluid Mechanics 942 (2022).
  * Savva and Bush (2009) N. Savva and J. W. Bush, Viscous sheet retraction, Journal of Fluid Mechanics 626, 211 (2009).
  * Koplik and Banavar (2000) J. Koplik and J. R. Banavar, Molecular simulations of dewetting, Physical Review Letters 84, 4401 (2000).
  * Koplik and Banavar (1993) J. Koplik and J. R. Banavar, Molecular dynamics of interface rupture, Physics of Fluids A: Fluid Dynamics 5, 521 (1993).
  * Zhao  _et al._ (2019) C. Zhao, J. Sprittles, and D. Lockerby, Revisiting the rayleigh–plateau instability for the nanoscale, Journal of Fluid Mechanics (2019).
  * Diemand  _et al._ (2014) J. Diemand, R. Angélil, K. K. Tanaka, and H. Tanaka, Direct simulations of homogeneous bubble nucleation: Agreement with classical nucleation theory and no local hot spots, Physical review E 90, 052407 (2014).
  * Kadau  _et al._ (2010) K. Kadau, J. L. Barber, T. C. Germann, B. L. Holian, and B. J. Alder, Atomistic methods in fluid simulation, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 368, 1547 (2010).
  * Todd and Daivis (2017) B. D. Todd and P. J. Daivis,  _Nonequilibrium molecular dynamics: theory, algorithms and applications_ (Cambridge University Press, 2017).
  * Kono  _et al._ (2014) S. Kono, T. Kaneko, and I. Ueno, Elevation of the temperature of liquid films caused by rapid rupturing, Physical Review E 90, 051004 (2014).
  * Tammaro  _et al._ (2021) D. Tammaro, V. Chandran Suja, A. Kannan, L. D. Gala, E. Di Maio, G. G. Fuller, and P. L. Maffettone, Flowering in bursting bubbles with viscoelastic interfaces, Proceedings of the National Academy of Sciences 118, e2105058118 (2021).
  * Manikantan and Squires (2020) H. Manikantan and T. M. Squires, Surfactant dynamics: hidden variables controlling fluid flows, Journal of Fluid Mechanics 892 (2020).
  * Smith (2013) E. Smith, On the Coupling of Molecular Dynamics to Continuum Computational Fluid Dynamics, PhD thesis, Imperial College London (2013).
  * (36) See supplemental material at […], for more details on the methodology, effect of initial conditions, temperature field, surface tension calculation and derivation of the mathematical framework newly developed.
  * Frankel and Mysels (1969) S. Frankel and K. J. Mysels, Bursting of soap films. II. Theoretical considerations, The Journal of Physical Chemistry 73, 3028 (1969).
  * Vrij (1966) A. Vrij, Possible mechanism for the spontaneous rupture of thin, free liquid films, Discussions of the Faraday Society 42, 23 (1966).
  * Vrij and Overbeek (1968) A. Vrij and J. T. H. G. Overbeek, Rupture of thin liquid films due to spontaneous fluctuations in thickness, Journal of the American Chemical Society 90, 3074 (1968).
  * Ivanov  _et al._ (1970) I. B. Ivanov, B. Radoev, E. Manev, and A. Scheludko, Theory of the critical thickness of rupture of thin liquid films, Transactions of the Faraday Society 66, 1262 (1970).
  * Taylor and Michael (1973) G. I. Taylor and D. H. Michael, On making holes in a sheet of fluid, Journal of Fluid Mechanics 58, 625 (1973).
  * Pandit and Davidson (1990) A. B. Pandit and J. F. Davidson, Hydrodynamics of the rupture of thin liquid films, Journal of Fluid Mechanics 212, 11 (1990).
  * Sanjay  _et al._ (2022) V. Sanjay, U. Sen, P. Kant, and D. Lohse, Taylor–Culick retractions and the influence of the surroundings, Journal of Fluid Mechanics 948, A14 (2022).
  * Vaynblat  _et al._ (2001) D. Vaynblat, J. R. Lister, and T. P. Witelski, Rupture of thin viscous films by van der waals forces: Evolution and self-similarity, Physics of Fluids 13, 1130 (2001).
  * Neél and Villermaux (2018) B. Neél and E. Villermaux, The spontaneous puncture of thick liquid films, Journal of Fluid Mechanics 838, 192 (2018).
  * Villermaux (2020) E. Villermaux, Fragmentation versus cohesion, Journal of Fluid Mechanics 898, P1 (2020).
  * Debrégeas  _et al._ (1995) G. Debrégeas, P. Martin, and F. Brochard-Wyart, Viscous bursting of suspended films, Physical Review Letters 75, 3886 (1995).
  * Roth  _et al._ (2005) C. Roth, B. Deh, B. Nickel, and J. Dutcher, Evidence of convective constraint release during hole growth in freely standing polystyrene films at low temperatures, Physical Review E 72, 021802 (2005).
  * Müller  _et al._ (2007) F. Müller, U. Kornek, and R. Stannarius, Experimental study of the bursting of inviscid bubbles, Physical Review E 75, 065302 (2007).
  * Trittel  _et al._ (2013) T. Trittel, T. John, K. Tsuji, and R. Stannarius, Rim instability of bursting thin smectic films, Physics of Fluids 25, 052106 (2013).
  * Mysels and Stikeleeather (1971) K. J. Mysels and J. A. Stikeleeather, The Bursting of Soap Films III. Spontaneous and Induced Bursting, Journal of Colloid and Interface Science 35, 159 (1971).
  * Moody and Attard (2001) M. Moody and P. Attard, Curvature dependent surface tension from a simulation of a cavity in a Lennard-Jones liquid close to coexistence, The Journal of Chemical Physics 115, 8967 (2001).
  * Wen  _et al._ (2021) J. Wen, D. Dini, H. Hu, and E. R. Smith, Molecular droplets vs bubbles: Effect of curvature on surface tension and tolman length, Physics of Fluids 33, 072012 (2021).
  * Irving and Kirkwood (1950) J. H. Irving and J. G. Kirkwood, The statistical mechanical theory of transport processes. IV. The equations of hydrodynamics, The Journal of Chemical Physics 18, 817 (1950).
  * Yeh and Hummer (2004) I.-C. Yeh and G. Hummer, System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions, The Journal of Physical Chemistry B 108, 15873 (2004).
  * Kirkwood and Buff (1949) J. G. Kirkwood and F. P. Buff, The statistical mechanical theory of surface tension, The Journal of Chemical Physics 17, 338 (1949).
  * Rayleigh (1880) L. Rayleigh, On the stability, or instability, of certain fluid motions, Proceedings of the London Mathematical Society 9, 57 (1880).
  * Lhuissier and Villermaux (2009a) H. Lhuissier and E. Villermaux, Soap films burst like flapping flags, Physical Review Letters 103, 054501 (2009a).
  * Lhuissier and Villermaux (2009b) H. Lhuissier and E. Villermaux, Destabilization of flapping sheets: The surprising analogue of soap films, Comptes Rendus Mecanique 337, 469 (2009b).
  * Brenner and Gueyffier (1999) M. P. Brenner and D. Gueyffier, On the bursting of viscous films, Physics of Fluids 11, 737 (1999).
  * Mysels and Vijayendran (1977) K. J. Mysels and B. Vijayendran, The bursting of soap films. 8. rim velocity in radial bursting, The Journal of Physical Chemistry 81, 731 (1977).
  * Nosé (1984) S. Nosé, A molecular dynamics method for simulations in the canonical ensemble, Molecular Physics 52, 255 (1984).
  * Hoover (1985) W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Physical Review A 31, 1695 (1985).
  * Oratis  _et al._ (2020) A. T. Oratis, J. W. Bush, H. A. Stone, and J. C. Bird, A new wrinkle on liquid sheets: Turning the mechanism of viscous bubble collapse upside down, Science 369, 685 (2020).
  * Hulshof (1901) H. Hulshof, The direct deduction of the capillarity constant $\gamma$ as a surface tension, Annalen der Physik 4, 165 (1901).
  * Mecke  _et al._ (1997) M. Mecke, J. Winkelmann, and J. Fischer, Molecular dynamics simulation of the liquid–vapor interface: The lennard-jones fluid, The Journal of Chemical Physics 107, 9264 (1997).