\section{Introduction}

Classical novae are the results of thermonuclear explosions of hydrogen on
the surfaces of white dwarfs (WDs) that accrete H-rich material from their low-mass MS or red-giant binary companions.
Depending  on the initial parameters of the binary system, the primary component can end up its evolution
as a CO or a more massive ONe WD.
The accreted material is heated by gravitational compression and, when the pressure at its base reaches
a critical value, H is ignited under partially degenerate conditions, which leads to a thermonuclear runaway (TNR).
The rapidly increasing temperature and energy generation rate trigger convection in the envelope, where
the hot CNO cycle transforms the stable CNO isotopes into their $\beta^+$-unstable counterparts, mainly
$^{13}$N, $^{14}$O, $^{15}$O, and $^{17}$F. From this moment on, the energy generation rate is limited by
the half-life times of the unstable isotopes, of order $10^2$\,--\,$10^3$ seconds.
The degenerate conditions are lifted soon after the temperature starts to rise during the early stages of
the TNR. Mass loss is driven by the TNR and results from the conversion of energy released by nuclear reactions into
kinetic energy, leading to $v_{\rm env} > v_{\rm esc}$. This will be supplement by additional mass loss from a wind (mainly).
The observed enrichment of the ejecta of novae in heavy elements is believed to be a signature of
mixing between the accreted envelope and WD.
Recent two- and three-dimensional nuclear-hydrodynamic simulations of a nova outburst have shown that
this mixing is most likely to be driven by the hydrodynamic instabilities
and shear-flow turbulence induced by steep horizontal velocity gradients at the bottom of the convection zone that
develops during the explosion \cite{cea11a,cea11b}. Therefore, we call it the convective boundary mixing (CBM).
The large mass fractions, on average 30\% \cite{gea98}, of CNO and Ne nuclei in the ejecta
confirm that nova outbursts and CBM occur on both CO and ONe WDs.

Recently, we have demonstrated that the state-of-the-art stellar evolution code
MESA\footnote{MESA is an open source code; it is described in detail in \cite{pea11} and can be downloaded
from {\bf http://mesa.sourceforge.net}.}
can successfully be used for 1D numerical simulations of nova outbursts occurring on CO WDs \cite{dea12}.
In a forthcoming paper, we will show that MESA can model the evolution of ONe novae
as well. In addition, the multi-zone parallel code MPPNP, which is a part of the NuGrid
research collaboration ({\bf www.nugridstars.org}), is used in the second paper
for post-processing computations of nucleosynthesis in CO and ONe novae. Our
MESA nova models serve as a background that provides MPPNP with necessary stellar structure and mixing parameter profiles.
In both MESA and NuGrid codes, convection is not approximated by instantaneous mixing but instead treated as
a diffusion process with diffusion coefficients proportional to local values of convective velocity and mixing length.
As a result of our cooperation with JINA and TRIUMF,
we have combined MESA and NuGrid tools into a framework that allows to produce up-to-date models of nova outbursts and
nucleosyntheis. The interaction between different parts of the network is controlled by user-friendly shell scripts.
It can be used for research and educational purposes.
In this preliminary report, we present the most interesting new results related
to CO and ONe nova outbursts that have been obtained with the Nova Framework.

\section{Results Related to CO Novae}

In the first paper \cite{dea12}, we have used the MESA stellar evolution code to construct
multicycle nova evolution sequences with CO WD cores. We have explored a range of WD
masses ($0.65\,M_\odot$, $0.85\,M_\odot$, $1.0\,M_\odot$, $1.15\,M_\odot$, and
$1.2\,M_\odot$) and accretion rates (from $10^{-11}M_\odot/\mbox{yr}$ to
$10^{-9}M_\odot/\mbox{yr}$) as well as the effect of different cooling times (different initial WD temperatures
and luminosities) before the onset of accretion. In addition, we have studied the dependence on the
elemental abundance distribution of accreted material and convective boundary mixing (CBM)
at the core-envelope interface. In the first case, which is commonly adopted in 1D nova
simulations, the CBM was mimicked by assuming that the accreted envelope had been
pre-enriched (pre-mixed) with 30\% of WD's material. In the second case, the CBM was implemented using a
diffusion coefficient that was exponentially decreasing with a distance below the bottom of
the convective envelope. This 1D approximation for CBM is supported by 3D hydrodynamic
simulations of He-shell flash convection in AGB stars \cite{hea11}. Our nova models with
such CBM display an enrichment of the accreted envelope with C and O from the underlying white
dwarf that is commensurate with observations.  We have compared our results with the previous
work and investigated a new scenario for novae with the $^3$He-triggered convection. We have
found that $^3$He can be produced \emph{in situ} in H-rich envelopes accreted with slow rates
($\dot{M} < 10^{-10}M_\odot/\mbox{yr}$) by cold ($T_{\rm WD} < 10^7$ K) CO WDs, and
that convection is triggered by $^3$He burning before the nova outburst in this case.

\section{Results Related to ONe Novae}

In the forthcoming paper,
we show that for ONe novae the CBM prescription with the exponentially decaying
diffusion coefficient produces maximum temperature evolution
profiles and nucleosynthesis yields in good agreement with those obtained using the pre-mixed
nova model (Fig.~\ref{fig:f1}) and that this is not true for CO novae.

We demonstrate that in ONe novae $^3$He can also be produced \emph{in situ}
in H-rich envelopes accreted with slow rates ($\dot{M} < 10^{-10}\,M_\odot/\mbox{yr}$) by cold ($T_{\rm WD} < 10^7$ K)
WDs, and that convection is triggered by $^3$He burning before the nova outburst in this case (Fig.~\ref{fig:f2}).

Additionally, we find that the interplay between the $^3$He
production and destruction in the H-rich envelope accreted with the intermediate rate,
$\dot{M} = 10^{-10}\,M_\odot/\mbox{yr}$, by the $1.15\,M_\odot$ ONe WD with the central temperature
$T_{\rm WD} = 15\times 10^6$ K
leads to the formation of a thick radiative buffer zone that separates the bottom of the convective
envelope from the WD surface (Fig.~\ref{fig:f3}).

\begin{figure}
\centering
\includegraphics[bb=60 225 480 695]{ne_nova_1.15_12_X_123_mixed_ovsh.eps}
\caption{A comparison of $T_{\rm max}$ evolution profiles (the upper panel) and final abundances (the lower panel) from two of our
         $1.15\,M_\odot$ ONe nova simulations with
         $T_{\rm WD} = 12$ MK and $\dot{M} = 2\times 10^{-10}\,M_\odot/\mbox{yr}$.
         The solid red curve and red circles are the results obtained with the 50\% pre-mixed
         initial abundances in the accreted envelope (a mixture of equal amounts of WD's and solar-composition
         materials), while the dashed blue curve and blue squares represent the case
         with the solar-composition accreted material and convective boundary mixing (CBM) modeled using the same prescription
         (the exponentially decaying diffusion coefficient) that we used in the $1.2\,M_\odot$ CO nova model in \protect\cite{dea12}.}
\label{fig:f1}
\end{figure}

\begin{figure}
\centering
\includegraphics[bb=75 200 480 695]{plot_3_ne_nova_1.3_7_B_3Heconv.eps}
\caption{A snapshot of profiles of various stellar structure parameters in the envelope of our $1.3\,M_\odot$ ONe nova model with
         $T_{\rm WD} = 7$ MK and $\dot{M} = 10^{-11}\,M_\odot/\mbox{yr}$ at the moment following the $^3$He ignition that
         has triggered convection. Like in the CO nova model with the cold WD and low accretion rate
         discussed in \protect\cite{dea12}, the low $T_{\rm WD}$ and $\dot{M}$ values favour the accumulation of $^3$He in a slope
         adjacent to the core-envelope interface followed by its ignition at a relatively low $T$,
         before the major nova outburst ensues.}
\label{fig:f2}
\end{figure}

\begin{figure}
\centering
\includegraphics[bb=75 200 480 695]{plot_3_ne_nova_1.15_15_A_buffer.eps}
\caption{A snapshot similar to that shown in Fig.~\protect\ref{fig:f2}, but for our $1.15\,M_\odot$ ONe nova model with
         $T_{\rm WD} = 15$ MK and $\dot{M} = 10^{-10}\,M_\odot/\mbox{yr}$. The interplay between the $^3$He production
         and destruction shifts the peak temperature away from the core-envelope interface, located at the leftmost steps on
         the abundance profiles (the middle panel), to a place where the $^3$He
         burning generates the maximum energy. This leads to the formation of a convective zone separated from
         the interface by a buffer zone.}
\label{fig:f3}
\end{figure}

\acknowledgments
The speaker appreciates helpful discussions with Lars Bildsten, Barry Davids,
Ami Glasner, Bill Paxton, Chris Ruiz, and Michael Wiescher.


\section{Introduction}
A central objective of non-equilibrium statistical mechanics is the description of the time evolution of thermodynamic properties of macroscopic systems that are prepared far from equilibrium. In principle one should think that it is always possible to derive an effective equation of motion, e.g. a Fokker-Planck equation, that governs the non-equilibrium behavior of thermodynamic quantities, from the underlying microscopics. In general this is, however, a formidable task and even more so in isolated systems where there is no apparent system-bath decomposition or an immediately manifest hierarchy in the coupling strengths between two or more subsystems. Nevertheless, these isolated systems often exhibit a generic relaxation dynamics, in the sense that for a subset of observables the non-equilibrium evolution leads to a steady state which is well-described by a thermodynamic ensemble.

For classical systems the underlying microscopic Hamilton equations of motion are non-linear in the dynamical variables. Generically, this leads to chaos such that the system ultimately explores the entire phase space. The occurrence of thermalization is then made plausible by invoking the Ergodic Hypothesis \cite{re:85}. In contrast, the microscopic equation describing the dynamics of a closed quantum system is a linear wave equation for a complex valued wavefunction and observables are represented by Hermitian operators. Due to these mathematical differences it is not obvious which equivalent arguments would lead to thermalization in this quantum context.

Given this, a first approach to study thermalization of quantum systems would involve a direct calculation of their microscopic dynamical evolution. However, determining the exact quantum dynamics of an interacting many-body system is a formidable task. The microscopic approach of solving Schr\"odinger's equation is usually restricted to very small systems, since the exponential growth of the Hilbert space with the number of particles also entails a growth in the mathematical complexity of the problem that quickly becomes unmanageable. Circumventing this problem requires sophisticated numerical tools and clever mathematical approaches. It is fair to say that our understanding of the equilibration and thermalization of closed quantum systems is still far less developed than for their classical counterparts.

Recently, there is an ever growing interest in addressing the question of when and how these systems relax towards an equilibrium state when prepared in a non-stationary initial state \cite{baal+:06,posh+:09,leol+:10,chki:10,biko+:10,gomu+:11,pode+:11,buda+:11,baci+:11,pose+:11,yu:11,sces:12,geho+:12,olgo+:12,atga+:12}. This development is, on the one hand, sparked by the availability of state of the art experiments - especially in the field of ultracold gases - where the investigated systems are extremely well isolated from the thermal environment and thus constitute almost ideal realizations of closed systems \cite{kiwe+:06,hole+:07,wibe+:10,trch+:12,grku+:12}. These experiments do not only constitute ideal test-beds for new theoretical ideas and concepts but they also point to the fact that there is growing interest for a quantum theory of the dynamics of closed systems to adequately interpret experimental observations made in these thermally isolated setups.

On the other hand, recent successes in the theory of equilibration and thermalization of closed quantum systems have fueled the interest in this topic. A very insightful idea was formulated with the Eigenstate Thermalization Hypothesis (ETH) that links thermalization to the spectral properties of the quantum Hamiltonian \cite{de:91,sr:94}. In essence, ETH states that in a thermalizing system any many-body eigenstate contains a thermal state, i.e., that the expectation value of an observable $\mathcal{O}$ calculated in an eigenstate with energy $\varepsilon$ coincides with the thermal average taken in the microcanonical ensemble at same energy. This conjecture has numerically been verified for a number of model systems \cite{ridu+:08,ecko:08,pahu:10}.

Despite its success the ETH does not provide an answer to how relaxation takes place and does not provide any information on the non-equilibrium evolution which eventually drives the system into a steady state. Here we take a closer look at the equilibration of a closed, interacting many-body quantum system in the time domain, expanding on a recent work \cite{atga+:12}. To this end, we consider a simple yet generic one-dimensional spin model, which belongs to the class of Ising models in a transverse field and can be mapped on a system of quantum hard rods that are coherently removed from/deposited on a lattice. For this system we show that it is possible to derive analytically a Master equation which captures the non-equilibrium dynamics of the number of hard rods. We undertake a comparison of this effective description with numerically exact simulations of the fully quantum problem and find excellent agreement. In this way, the results of this work also contribute to recent efforts to derive effective diffusion equations that describe the non-equilibrium dynamics in closed quantum systems \cite{atga+:12,tiva+:12p,ge:12pc}.

This paper is organized as follows: In  \sref{sec:hamiltonian} we introduce a spin Hamiltonian that represents a quantum version of a lattice gas of hard rods. Subsequently, in \sref{sec:config_network}, we introduce a specific graphical representation of the Hilbert space of the quantum model - the configuration network - to conveniently illustrate the dynamics of this system. Using the insights from a statistical analysis of the network, we present the derivation of a Master equation that describes the dynamics of the number distribution function of hard rods in \sref{sec:master_equation}. In particular, we show that the rate coefficients of this Master equations can be calculated analytically. Moreover, we numerically investigate the time evolution and steady state distribution of the density of hard rods. Following this analysis, we compare in \sref{sec:numerics} the predictions of the Master equation with  the results obtained by numerically solving the many-body Schr\"odinger equation for our model system and find excellent agreement. Conclusions and an outlook are provided in  \sref{sec:outlook}. Unless stated otherwise, we set $\hbar =1 $ in this work.

\section{The system} \label{sec:hamiltonian}
The quantum system in which we are interested here is a chain of $L$ equally spaced, interacting spin-$1/2$ particles at nearest neighbor distance $a$ in a transverse magnetic field of strength $\Omega$. Its Hamiltonian is given by
\begin{eqnarray}
H_{\mbox{\scriptsize spin}} = H_\Omega+H_V = \Omega \sum_{j=1}^L \sigma_j^x + V \sum_{j=1}^L \sum_{i=1}^\lambda n_j n_{j+i},
\label{eqn:spin_H}
\end{eqnarray}
where, $\sigma_j^x = \left(\left|\uparrow \right> \left< \downarrow \right| + \left|\downarrow  \right> \left< \uparrow \right| \right)_j$ is a Pauli spin matrix effectuating spin flips of the $j$-th spin and $n_j = \left( \left|\uparrow  \right> \left< \uparrow \right| \right)_j$ is a projector on the up-state (excited state) of the $j$-th spin. The spin-spin interaction of strength $V$ depends on the spin-state as well as on the separation of the spins:  two spins interact (i) if both of them are simultaneously in the spin-up state and (ii) if they are separated by a distance that does not exceed a critical distance $\lambda$, which we will refer to as \emph{blockade radius}. We use this terminology because in this work we focus specifically on the regime $V/\Omega \to \infty$. This means that it is energetically forbidden to have two or more up-spins located within the critical radius $r_{\mathrm{c}}=\lambda\,a$, i.e. an up-spin blocks the excitation of spins to the up-state in its vicinity. The model defined by the Hamiltonian (\ref{eqn:spin_H}) is in fact of practical relevance, as it has been shown to capture the quantum dynamics of strongly interacting, laser-driven Rydberg atoms trapped in an optical lattice \cite{viba+:11,scch+:12p}. For more details we refer the reader to references \cite{olmu+:10,le:11,atle:12,homu+:12p}.

The Hilbert space of the system is spanned by all classical spin configurations compatible with the blockade condition. In \fref{fig:System}(a) we provide examples of permitted spin configurations for the cases $\lambda=1$ and $\lambda=2$. Instead of using spins we employ a description of the basis states in terms of hard rods. Each configuration of spins can be uniquely mapped into and arrangement of hard rods of length $\lambda +1$, as illustrated in \fref{fig:System}(a). A spin-flip from the down- to the up-state corresponds to the deposition a hard rod, while the opposite process removes a hard rod from the lattice. Our aim is to analyze the non-equilibrium dynamics of the hard rods effectuated by $H_{\mathrm{spin}}$ after a quench. In particular we seek an effective evolution equation for the probability density $p_n (t)$, which describes the probability of having the lattice occupied by $n$ hard rods at time $t$. The first step towards such equation is to gain an understanding of the Hilbert space structure. Here it is convenient to employ a graphical representation of the Hilbert space as a configuration network. We will see that the structure of the effective evolution equation for the $p_n$ will crucially depend on the statistical properties of the network. In the following, we will construct the network and analyze its properties in detail.

\begin{figure}
\includegraphics[width=1\textwidth]{figure1a_1b.png}
\caption{(a) Representation of configurations in terms of spins (arrows) and the hard rods (grey rectangles). Each up spin inhibits further excitations on $\lambda$ neighboring sites on both its left and right. This blockade effect can equivalently be illustrated by hard rods on the lattice, each occupying $\lambda+1$ sites, and spin flips correspond to placing/removing hard rods on the lattice. (b) An example of the configuration network for $L=8$ and $\lambda =1$. Each node represents a classical arrangement of $n$ hard rods $| n \mathcal{C}_n \rangle$ (see text). The links between nodes connect configurations that can be converted into one another by the deposition/removal of a hard rod. Since we are interested in the dynamics in number space we do not resolve individual states but consider the time evolution of the operator $P_n$ that projects on all basis states contained in the $n$-th column. The equation of motion for $\langle P_n \rangle$ depends directly on the structure of the graph, specifically on the abundance of loops (green), reflections (orange) and transmissions (red).}
\label{fig:System}
\end{figure}

\section{Configuration network}\label{sec:config_network}
\subsection{Structure}
To construct the configuration network we introduce the basis states $| n \mathcal{C}_n \rangle$. Each of these states represents a specific classical arrangement $\mathcal{C}_n$ of $n$ hard rods. These states satisfy the completeness $\sum_{n\mathcal{C}_n}\left|n\mathcal{C}_n\right>\left<n\mathcal{C}_n\right|= \dblone$ and orthonormality relation $\left<n\mathcal{C}_n |m\mathcal{K}_m\right> = \delta_{nm}\delta_{\mathcal{C}_n\mathcal{K}_m}$. We interpret each of the microstates $| n \mathcal{C}_n \rangle$ as a node of a network. By grouping configurations containing the same number of hard rods into columns we obtain the network structure depicted in  \fref{fig:System}(b). The time-evolution of the system can then be imagined as a temporal change in the occupation of these nodes. Dynamics is introduced through the Hamiltonian $H_\Omega$ which leads to transitions between microstates that are represented as edges of the network. Since $H_\Omega$ causes only single spin flips, nodes in neighboring columns are only linked directly, if their corresponding microstates can be converted into one another by the removal/deposition of one hard rod. For example, setting $\lambda=1$ the state $\left|\downarrow\uparrow\downarrow\downarrow\uparrow\downarrow\downarrow\uparrow\downarrow\right>$ in the $n=3$ column is directly linked with $\left|\downarrow\uparrow\downarrow\downarrow\downarrow\downarrow\downarrow\uparrow\downarrow\right>$, but not with $\left|\downarrow\downarrow\uparrow\downarrow\downarrow\downarrow\uparrow\downarrow\right>$ in the $n=2$ column. In the following we will analyze in detail the properties of the configuration network.

\subsection{Properties of the configuration network}
The most basic properties that define the structure of our configuration network are the number of columns and the number of nodes within each column. Fixing the length of the system to $L$ sites, and applying periodic boundary condition, the maximum number of hard rods that each occupy $\lambda+1$ sites (i.e., the blockade radius is $\lambda$), which can be placed on the lattice is $\lfloor L/(\lambda+1) \rfloor$, where $\lfloor x \rfloor$ denotes the closest integer smaller or equal to $x$. The index counting the number of hard rods can thus take the values $n = 0,1,..., \lfloor L/(\lambda+1) \rfloor$. The number of microstates $\nu_{n}$ contained in the $n$-th column is given by the number of ways in which one can distribute $n$ indistinguishable hard rods of length $\lambda+1$ over $L$ lattice sites. This is a standard combinatorial problem with solution
\begin{equation}
\nu_{n} = \frac{ L (L-1-\lambda n)!}{n!(L-(\lambda+1)n)!} \, .
\end{equation}
Having determined the properties of the ``backbone'' of the network we now turn to assessing the linkage of the nodes. In particular, we calculate the mean number of different possibilities $T_{n\to n\pm 1}$ to go from a state with $n$ hard rods to one in the adjacent columns. This quantity can be expressed as $T_{n\to n\pm 1} = c_{n, n \pm 1}/\nu_n$, where $c_{n, n \pm 1}$ denotes the \emph{total} number of links between columns $n$ and $n \pm 1$ and hence  $c_{n \pm 1 , n} = c_{n , n \pm 1}$. (Note that this symmetry does not hold for the $T_{n \to m}$.) Moreover, we know that $T_{n \to n-1} = n$, as in a configuration of $n$ hard rods there are $n$ possibilities for removing one hard rod reaching a state with $n-1$ hard rods. Using these relations the total number of links between two columns evaluates to
\begin{eqnarray}
c_{n , n-1} &= &T_{n\to n-1}\nu_n =  \frac{ L (L-1-\lambda n)!}{(n-1)!(L-(\lambda+1)n)!}, \nonumber \\
c_{n+1, n } &= &T_{n+1\to n}\nu_{n+1} = \frac{ L (L-1-\lambda (n+1))!}{n!(L-(\lambda+1)(n+1))!},\nonumber \\
\end{eqnarray}
which can then be used to calculate $T_{n \to n+1}$.

\begin{figure}
\centering
\includegraphics[width=0.5\columnwidth]{figure2a_2c.png}
\caption{(a,b) Graphical illustration of the reflection processes. Two configurations can at most be connected by two reflection pathways (a). However, since rods must not overlap, the path containing an additional hard rod during the intermediate step is often forbidden (b). (c) For the numerical study we fix the physical system length $l$ and the critical radius $r_\mathrm{c}$. The parameter $\lambda$ is varied by increasing the number of lattice sites.}\label{fig:reflections}
\end{figure}

Let us continue by analyzing the second order processes shown in \fref{fig:System}(b), i.e. loops, reflections and transmissions, as they will enter in the derivation of the Master equation. Selecting a node from the network, the number of loop transitions from that node equals the number of links that this node has with other nodes. Therefore, the mean number of loop transitions from a state with $n$ hard rods is given by
\begin{equation}
\overline{N_{\mathrm{loop}}^{(n)}} = T_{n\to n+1} + T_{n \to n-1}.
\label{eq:diagonal}
\end{equation}
Reflections connect two configurations that contain the same number of hard rods but differ in the position of exactly one hard rod [cf. \fref{fig:reflections} (a,b)]. Two configurations that are randomly selected from one column will typically differ by the positioning of more than one hard rod and are thus not connected by a reflection. If two microstates happen to be connected there are at most two paths as illustrated in \fref{fig:reflections}: a deexcitation followed by an excitation or vice versa. For $n \gg 1$ (at high density) there is often even only one path available [\fref{fig:reflections}(b)]. Similar considerations can also be made for transmission diagrams, i.e., the average number of paths connecting two microstates containing $n$ and $n\pm2$ hard rods is $\sim 1$. In fact the mean number of  reflections (transmissions) $\overline{N_{\mathrm{refl}}^{(n)}}$ ($\overline{N_{\mathrm{trans}}^{(n)}}$) between two \emph{randomly} selected states can be calculated analytically:
\begin{eqnarray}
\overline{N_{\mathrm{refl}}^{(n)}} &=& \frac{T_{n+1 \to n}(T_{n+1\to n}-1)\nu_{n+1}}{\nu_n(\nu_n-1)} +\frac{T_{n-1 \to n}(T_{n-1\to n}-1)\nu_{n-1}}{\nu_n(\nu_n-1)}, \nonumber\\
\overline{N_{\mathrm{trans}}^{(n)}} &= & \frac{T_{n+2 \to n+1}T_{n+1\to n}+T_{n-2\to n-1}T_{n-1\to n}}{\nu_n}.\nonumber \\
\label{eqn:average_transitions}
\end{eqnarray}
To illustrate that loop transitions are far more abundant than reflections and transmissions we present some numerical examples in \tref{tab:loop}. Here we compare $\overline{N_{\mathrm{loop}}^{(n)}}$, $\overline{N_{\mathrm{refl}}^{(n)}}$ and $\overline{N_{\mathrm{trans}}^{(n)}}$ for a number of lattice and hard rod sizes. This leads to two observations. First, the relative weight of loop transitions largely increases with increasing systems size $L$, as due to the larger dimension of the Hilbert space each state can have more connections to other configurations. Second, the probability that two states within a column of the network are connected by a reflection is vanishingly small in the ''bulk'' of the network ($1 \ll n \ll \lfloor L/(\lambda+1) \rfloor$). The same is also true for transmissions between columns $n$ and $n \pm 2$. Note, that near the boundaries of the network, i.e. columns close to the maximum/minimum $n$-value, the condition $\overline{N_{\mathrm{loop}}^{(n)}}\gg \overline{N_{\mathrm{trans}}^{(n)}},\overline{N_{\mathrm{refl}}^{(n)}}$ are less well satisfied. However, this concerns only an exponentially small subset of states forming the configuration network. These two observations on the statistics of the configuration network are of central importance in the derivation of the effective Master equation for $p_n(t)$ that we are going to present in the following section.

\begin{table}
\centering
  \begin{tabular}{|c |c|c|| c| c|c|}
  \hline
     L& $\lambda$ & $n$ & $ \overline{N_{\mathrm{loop}}^{(n)}}$ & $\overline{N_{\mathrm{refl}}^{(n)}}$ & $\overline{N_{\mathrm{trans}}^{(n)}}$  \\ \hline
     & & 2 & 8.22 & 0.54 & 2.22\\
     \cline{3-6}
     \raisebox{1ex}[0pt]{12} &  \raisebox{1ex}[0pt]{1} & 4 & 5.71 & 0.17& 0.51 \\
     \hline
      &   & 2 & 296 & 0.026& $2$ \\ \cline{3-6}
      &  \raisebox{1ex}[0pt]{1}  & 75 & 174.78 & $1.27\times10^{-57}$& $1.39\times10^{-57}$ \\ \cline{2-6}
      &   & 2 & 280.1 & 0.026 & $2$\\ \cline{3-6}
      &  \raisebox{1ex}[0pt]{5}  & 25 & 93.39 & $1.64\times10^{-27}$& $3.18\times10^{-27}$ \\ \cline{2-6}
       &   & 2 & 244.78 & 0.025&$2$ \\ \cline{3-6}
      \raisebox{7ex}[0pt]{300}  &  \raisebox{1ex}[0pt]
      {14}  & 10 & 68& $3.11\times10^{-13}$&$ 1.73\times10^{-12}$ \\      \hline
  \end{tabular}
  \caption{Average number of connection between nodes established by the three types of second order processes depicted in \fref{fig:System}(b). We have chosen a relatively small lattice size ($L=12$) and a large one ($L=300$) for illustration.}   \label{tab:loop}
\end{table}

\section{Master Equation}\label{sec:master_equation}
In this section we derive an equation of motion which describes the evolution of number of hard rods on the lattice as function of time. This means that we are not interested in the actual population of individual nodes of the graph shown in \fref{fig:System}(b) but rather in the probability $p_n(t)$ of the system to reside in a specific column $n$ at time $t$. We will see that this eventually leads to a Master equation that has a steady state $p_n(t \to \infty)$ which is proportional to the number of classical configurations contained in a specific column. This strongly suggests that this steady state corresponds to a microcanonical equilibrium state in which all arrangements of hard rods occur with equal probability.

\subsection{From the von-Neumann equation to the Master equation}\label{sec:master_equation_deriv}
The probability $p_n$ is defined as $p_n=\mbox{Tr} \rho(t) P_n$, where $\rho(t)$ is the density matrix of the system and $P_n = \sum_{\mathcal{C}_n} \left|n\mathcal{C}_n\right>\left<n\mathcal{C}_n\right|$ is a projector which projects onto the subspace spanned by all microstates contained in the $n$-th column of the configuration network [see \fref{fig:System}(b)]. Throughout this work we are interested in a situation where the initial state of the system $\rho(0)$ contains a fixed number of hard rods, i.e. $\left[P_n,\rho(0)\right]=0$.

To begin with the formal derivation of the Master equation, let us momentarily return to the case of finite $V$ and transform the Hamiltonian (\ref{eqn:spin_H}) into the interaction picture with respect to $H_V$. As shown in \ref{sec:integrated_v_N_equation} one can then write the evolution equation for $p_n = \left<P_n\right> = \mbox{Tr} \rho(t)P_n$ as
\begin{eqnarray}
\partial_t {{\left<P_n\right>}}_t &=&   -\int_0^t ds\, \mbox{Tr} \{ P_n H_I(t) H_I(s) \rho(s)    +H_I(s)H_I(t)P_n\rho(s)
\nonumber \\ &&- H_I(s)P_n H_I(t)  \rho(s)  -H_I(t)P_n H_I(s) \rho(s)\},
\label{eqn:integral_von-Neumann_expand}
\end{eqnarray}
with
\begin{eqnarray}
H_I(t) &=& \sum_{n\mathcal{C}_n, m\mathcal{K}_m} \left<n\mathcal{C}_n|H_\Omega|m\mathcal{K}_m\right> e^{-i(\omega_{m\mathcal{K}_m}-\omega_{n\mathcal{C}_n})t}\left|n\mathcal{C}_n \right> \left<m\mathcal{K}_m\right|,
\label{eqn:inter_pic}
\end{eqnarray}
where $\omega_{n\mathcal{C}_n}=\left<n\mathcal{C}_n|H_V|n\mathcal{C}_n\right>$ denotes the configuration energy of the spin configuration $|n\mathcal{C}_n \rangle$. In order to evaluate the integral in the still exact evolution \eref{eqn:integral_von-Neumann_expand} we make the replacement  $\rho(s)\to\rho(t)$, which effectively amounts to a second order approximation. At this point the quality of this approximation is not clear but we will provide a numerical justification of this step \emph{a posteriori}.
Returning to the ideal blockade regime ($V/\Omega \to \infty$) the fact that each allowed spin configuration has the same configuration energy, $\omega_{n\mathcal{C}_n} = 0$, removes the exponential factor of the interaction picture Hamiltonian making it time-independent ($H_I(t) \to H_{\Omega}$). Thus, the time-integration in \eref{eqn:integral_von-Neumann_expand} can be carried out and simply amounts to a multiplication by $t$. The steps result in the second order equation
\begin{equation}
\partial_t {{\left<P_n\right>}}_t = -t \mathrm{Tr} \left[
P_n H^2_{\Omega} + H^2_{\Omega} P_n - 2 H_{\Omega} P_n H_{\Omega}
\right] \rho(t),
\label{eq:after_integration}
\end{equation}
that depends on products of the form
\begin{equation}
\sum_ {m\mathcal{K}_m } \left< n\mathcal{C}_n|H_\Omega |m\mathcal{K}_m\right>\left< m\mathcal{K}_m|H_\Omega |p\mathcal{L}_p\right> = \left<n\mathcal{C}_n\right| H_\Omega^2 \left|p\mathcal{L}_p\right>  \nonumber
\label{eqn:matrix_elements}.
\end{equation}
These matrix elements are the second order processes discussed in the previous section. Diagonal elements correspond to loops and off-diagonal ones to either reflections or transmission. From the analysis of the graph we know that loop transitions are far dominant and that $H_\Omega^2$ is thus approximately diagonal, hence
\[
\left<n\mathcal{C}_n\right| H_\Omega^2 \left|p\mathcal{L}_p\right> \approx
\left<n\mathcal{C}_n\right| H_\Omega^2 \left|n\mathcal{C}_n\right> \delta_{np} \delta_{\mathcal{C}_n \mathcal{L}_p}.
\]
The final approximation consists of neglecting the variation of the matrix elements of $H^2_{\Omega}$ with $\mathcal{C}_n$ within a given column of the configuration network, by replacing the matrix element with its average taken over all microstates within a column
\begin{equation}
\left<n\mathcal{C}_n\right| H_\Omega^2 \left|n\mathcal{C}_n\right> \to
\Omega^2 \overline{N_{\mathrm{loop}}^{(n)}} =
\Omega^2 \left( T_{n \to n + 1} + T_{n  \to n-1 } \right).
\end{equation}
Inserting these approximations into \eref{eq:after_integration} and writing $p_n (t) = \mathrm{Tr} P_n \rho(t) = \sum_{\mathcal{C}_n} \mathrm{Tr} \left|n\mathcal{C}_n\right> \left<n\mathcal{C}_n\right| \rho(t)$ we arrive at the desired Master equation for the dynamics of the probability distribution $p_n$ of finding $n$ hard rods in the system at time $t$
\begin{eqnarray}
\partial_t {p}_{n} \left(t\right) &=& 2{\Omega}^2 t \left[ T_{n+1\rightarrow n} p_{n+1}(t) + T_{n-1\rightarrow n} p_{n-1}(t)\right] \nonumber \\
 &&- 2{\Omega}^2 t \left[ T_{n\rightarrow n-1}  + T_{n\rightarrow n+1} \right] p_{n}(t) .
\label{eqn:master_eqn}
\end{eqnarray}
This Master equation has a stationary solution despite the explicit appearance of the time variable on the right hand side. This can be seen by transforming to a new time variable $\tau = \Omega t^2$ which removes the explicit time-dependence. The rate coefficients of the Master equation (\ref{eqn:master_eqn}) obey detailed balance, i.e., $p_n^{\mathrm{eq}}T_{n \to n+1} = p_{n+1}^{\mathrm{eq}}T_{n+1 \to n}$, where  $p_n^{\mathrm{eq}}$ denotes the steady state distribution
\begin{equation}
p_n^{\mathrm{eq}}= \frac{\nu_n}{\sum_{n=0}^{L/(\lambda+1)}\nu_n}.
\label{eqn:counting_dist}
\end{equation}
This distribution is proportional to the number $\nu_n$ of arrangements of $n$ hard rods, suggesting that each node of the network is populated with equal probability, as one would expect from a "maximum entropy state".

As shown in reference \cite{atga+:12}, \eref{eqn:master_eqn} is a discretized Fokker-Planck equation in excitation number space. Its drift coefficient is proportional to $T_{n\rightarrow n-1}  - T_{n\rightarrow n+1}$ and its diffusion coefficient proportional to $T_{n\rightarrow n-1}  + T_{n\rightarrow n+1}$. The Master equation thus describes the excitation dynamics of the system as diffusion between the columns of the configuration network shown in \fref{fig:System}(b).

\subsection{Time Evolution and Steady State of the Master Equation} \label{sec:time_ev_and_steady_state}
\begin{figure}
\centering
\includegraphics[width=0.5\columnwidth]{figure3.png}
\caption{Time evolution of the distribution function $p_n(t)$ as a function of the hard rod density $ n\lambda /L$ for $\lambda=1$ (nearest neighbor blockade). The system considered has a lattice length to blockade radius ratio $l/r_{\mathrm{c}}=120$. From left to right $p_n$ is shown at: $\Omega t=0.02 \mathrm{ (red)},0.1,0.2,0.3,0.4,0.5,0.6,0.8$ and $\Omega t=2$ (blue, corresponding to the steady state distribution). The initial distribution function used is $p_0 = \delta_{n,0}$.
}
\label{fig:Master_n.n.}
\end{figure}

\begin{figure}
\includegraphics[width=1\columnwidth]{figure4a_4d.png}
\caption{
Dependence of the time evolution and steady state properties on $\lambda$.
(a) Time-evolution of the average hard rod density $ \langle n \rangle \lambda /L$ starting from the configuration without hard rods for different  numbers of blockaded sites: $\lambda =1$ (blue), $\lambda =9$ (green) and $\lambda =59$ (red). (b) Steady state values of $\langle n \rangle \lambda /L$ as function of $\lambda$. (c) Steady state distribution $p_n$ as function of the hard rod density $ n \lambda /L$ for the three values of the number of blockaded sites used in (a). (d) Steady state value of the parameter Q as function of the number of blockaded sites. For more details, see text.
}
\label{fig:Master_diff_lambda}
\end{figure}

We will now study the time evolution of the distribution function $p_n$ by numerically solving \eref{eqn:master_eqn}. We perform the following analysis by fixing both the physical lattice length $l$ and the critical radius $r_{\mathrm{c}}$. The ratio $l/r_{\mathrm{c}}$ is then also fixed. We choose a value of $l/r_{\mathrm{c}}=120$ in the numerical examples of this subsection, which means that at most $120$ hard rods can be placed on the lattice. What we change is the number of lattice sites occupied by a hard rod, i.e. $\lambda+1$, which is done by varying the total number $L$ of lattice sites. This is illustrated in \fref{fig:reflections}(c). We note that for $\lambda \to \infty$ this procedure results in taking the continuum limit $a/r_{\mathrm{c}} \to 0$.

\Fref{fig:Master_n.n.} shows the time evolution of $p_n$ as a function of the hard rod density $n\lambda /L$ for $\lambda=1$ starting from the initial state with no hard rods, i.e., $p_n (0) = \delta_{n,0}$. The time evolution of the distribution function towards its steady state (indicated in blue) shows the characteristic features expected from a system following a Fokker-Planck evolution: $p_n$ gets broader due to diffusion between the columns of the configuration network and drifts towards its equilibrium position, where the rate coefficients for further excitation ($T_{n\to n+1}$) and de-excitation ($T_{n\to n-1}$) become equal.

\Fref{fig:Master_diff_lambda} summarizes the behavior of the solution of the Master equation for varying $\lambda$. Panel (a) shows the time evolution of the mean hard rod density $ \langle n \rangle \lambda /L$ starting from an empty lattice. For short times the hard rod density exhibits a quadratic time dependence, $ \langle n \rangle \lambda /L \propto t^2$. This behavior is induced by the explicit linear time dependence of the right hand side of the Master equation. For long times $ \langle n \rangle \lambda /L$ saturates to its steady state value, which is reached faster the larger $\lambda$. Furthermore, the hard rod density in the steady state gets larger as $\lambda$ increases [\fref{fig:Master_diff_lambda}(b)]. Interestingly, however, the fluctuations around the mean density $ \langle n \rangle \lambda /L$ decrease with increasing $\lambda$, which can be seen in \fref{fig:Master_diff_lambda}(c), where the full steady state distribution function $p_n$ is depicted for three different values of $\lambda$. This behavior is further quantified in \fref{fig:Master_diff_lambda}(d), where we show the steady state value of the ratio $Q \equiv (\langle n^2 \rangle - \langle n \rangle^2)/\langle n \rangle$, which relates the width of the distribution function to its mean, as a function of $\lambda$. The ratio asymptotically approaches zero with increasing $\lambda$, i.e. increasing number of lattice sites occupied by a hard rod [cf.\ \fref{fig:Master_diff_lambda}(d)]. This surprising feature is the result of an entropic effect which has recently been discussed in reference \cite{atle:12}.

\section{Comparison with the Exact Quantum Evolution}\label{sec:numerics}

\begin{figure}
\centering
\includegraphics[width=1\columnwidth]{figure5a_5b.png}
\caption{
Comparison of the exact quantum dynamics with the predictions of the Master equation (\ref{eqn:master_eqn}) for different $\lambda$. (a) Time evolution of the hard rod density and (b) the distribution functions at $\Omega t =20$.  The vertical dashed lines indicate the maximum possible hard rod density for each $\lambda$, which is given by $\lambda/(\lambda+1)$. From top to bottom, the blockade radius is increased as $\lambda = 1,2,3$ while the ratio $l/r_{\mathrm{c}}=10$ remains fixed.  In all plots, the black curves are the solutions of the Master equation. The colored curves show the solution of the Schr\"odinger equation using three randomly chosen initial states $| n_0 \mathcal{C}_{n_0} \rangle$ with the same number $n_0$ of initial hard rods (from top to bottom $n_0 = 2,2,3$).
}
\label{fig:blockade_dependence}
\end{figure}

After having discussed the main features of the time evolution and the steady state of the Master equation we will now compare its predictions to the exact quantum dynamics of the system. To this end we have numerically solved  Sch\"odinger's equation with the Hamiltonian (\ref{eqn:spin_H}) in the limit of $V/\Omega \to \infty$ for up to $L=30$ sites. The left column of \fref{fig:blockade_dependence} shows the numerically exact quantum evolution of the hard rod density together with the prediction of the Master equation (in black). The ratio of system length to critical radius is fixed to  $l/r_{\mathrm{c}}=10$ and we choose $\lambda = 1, 2, 3$ from the top to the bottom panel. Note that the dimension of the Hilbert space increases from top to bottom. The differently colored curves in each panel show solutions to the Schr\"odinger equation starting from randomly chosen initial states $| n_0 \mathcal{C}_{n_0} \rangle$ with the same number of hard rods $n_0$ but different spin configurations $\mathcal{C}_{n_0}$. The initial number of hard rods is chosen such that the initial state lies in a region of the configuration network with large connectivity, i.e., where the statistical assumptions underlying the derivation of the Master equation are well met. The right column shows the corresponding probability distributions $p_n$ at $\Omega t =20$. For both $ \langle n \rangle \lambda /L$ and $p_n$ the agreement between the results of the exact quantum calculation and the prediction of the Master equation is remarkably good. In particular, for long times the results of the full quantum calculation and solution of the Master equation only differ by roughly one per cent. For short times the quadratic time dependence of $ \langle n \rangle \lambda /L$ as well as its dependence on $\lambda$ are well reproduced (see insets). For longer times the full quantum solutions exhibit oscillations around the equilibrium value of $ \langle n \rangle \lambda /L$. Being a simple rate equation our Master equation does not reproduce these. However, the quantum oscillations decrease with increasing dimension of the Hilbert space. For $\lambda =1$  this behavior was also reported in reference \cite{atga+:12}.

At long times we observe a small systematic offset of $ \langle n \rangle \lambda /L$ obtained from the exact numerics from the steady state values predicted by the Master equation. In the cases shown in \fref{fig:blockade_dependence}(a), where the initial state contains fewer hard rods than the equilibrium state, the quantum calculations suggest a sightly lower value of the average hard rod density at long times. In contrast, for initial states with a higher hard rod density than the equilibrium value the quantum results lie sightly above the prediction of the Master equation. Due to this systematic dependence on the initial state, we attribute this small offset to the presence of memory effects in the quantum dynamics that were completely neglected in the derivation of the Master equation. Furthermore, the data of the full quantum calculation for $\lambda=3$  exhibit a very slow drift of $ \langle n \rangle \lambda /L$ towards the steady state of our rate equation. This effect is seen more clearly in \fref{fig:long_time_evolve}, where we follow the time evolution of the hard rod density for $\lambda =4$ to much longer times. The observed shift might be indicative of a pre-equilibration process in the quantum system, in which $\langle n \rangle \lambda /L$ quickly reaches a quasi-equilibrium state, which then very slowly equilibrates to the ``true'' steady state. However, the drift might also stem from a long wavelength oscillation present in the full quantum calculation due to the finite size of the system. Since the full quantum calculations are limited to small system sizes it is difficult to further explore this effect, which seems to be more pronounced with increasing $\lambda$.
\begin{figure}
\centering
\includegraphics[width=1\columnwidth]{figure6.png}
\caption{Long time quantum evolution of the hard rod density starting from a randomly chosen initial state with $n_0 =3$ hard rods for $\lambda = 4$ and $l/r_{\mathrm{c}} = 9$. The evolution predicted by the Master equation is plotted in black. To highlight the shift, the time evolution from $\Omega t=0.5$ to $\Omega t=20$ is enlarged and shown in the left inset, and the evolution from $\Omega t=80$ to $\Omega t=100$ is enlarged and shown in the right inset.}
\label{fig:long_time_evolve}
\end{figure}

Let us finally return to the observation made in \fref{fig:blockade_dependence}(b) that distribution functions $p_n$ calculated fully quantum mechanically and using the Master equation agree very well at long times. In order to quantify the degree of agreement  we use the following overlap measure \cite{bezy:06},
\begin{equation}
\mathcal{D}=1-\frac{1}{2}\sum_{n=0}^{L/(\lambda+1)} |\bar{p}_n-p_n^{\mathrm{eq}}| \, .\label{eq:overlap}
\end{equation}
Here $p_n^{\mathrm{eq}}$  denotes the steady state solution of the Master equation and $\bar{p}_n$ is the equilibrium distribution obtained from the solution of the Schr\"odinger equation, time averaged over a time interval $\Delta t$, i.e.,  $\bar{p}_n = \int_{t}^{t+\Delta t} \mathrm{d} \tau \, p_n (\tau) / \Delta t$. The distribution functions are identical when $\mathcal{D}=1$, and  $\mathcal{D}=0$ for distribution functions that are completely non-overlapping. For the three situations discussed in \fref{fig:blockade_dependence} we have selected $100$ randomly chosen initial spin configurations for $L=10, \lambda=1$ and $500$ randomly chosen initial spin configurations for $L=20, \lambda=2$ and $L=30, \lambda=3$, and have used them as initial states of the quantum evolution. The time average in order to compute $\bar{p}_n$ was taken over the interval $\Delta t = [20/\Omega, 40/\Omega]$.
\begin{figure}
\centering
\includegraphics[width=0.9\columnwidth]{figure7.png}
\caption{
Histogram of the parameter $\mathcal{D}$, defined in \eref{eq:overlap}, for the three parameter sets used in \fref{fig:blockade_dependence}. The numbers in the parentheses show the dimension of the Hilbert space and the number of initial states used to create the histogram, respectively. For the bottom panel we also show the distribution of the number of hard rods contained in the initial states as an inset. The black dashed line with blue cross marks in the inset shows the microcanonical steady state distribution of the number of hard rods for comparison.
}
\label{fig:histogram}
\end{figure}
We have collected the results of these simulations in the histograms shown in \fref{fig:histogram}. Here we see that for the vast majority of initial conditions $\mathcal{D}$ is close to one with only a few outliers. This indicates that indeed equilibration is largely independent of the initial state.
In order to demonstrate that the used set of initial conditions is representative we have looked at the distribution of the number of hard rods contained in the initial states. As an example the inset in \fref{fig:histogram} shows this distribution for the parameters of the bottom panel. Comparing this to the microcanonical steady state distribution we see that indeed each number is represented with the correct weight. Finally, looking at the histograms shown in \fref{fig:histogram} as a function of the dimension of the Hilbert space we see that the system equilibrates better with increasing number of available microstates, confirming similar observations made in \cite{atga+:12} for the case $\lambda=1$. In the thermodynamic limit  the Master equation (\ref{eqn:master_eqn}) therefore indeed provides an excellent description of the non-equilibrium dynamics in particle number space.

\section{Conclusions and Outlook}\label{sec:outlook}
\begin{figure}
\centering
\includegraphics[scale = 0.7]{figure8.png}
\caption{Configuration network of a system with
$L=6$ sites, $\lambda=1$ and finite ratio $V/\Omega$. Different layers
represent different interaction energies and adjacent layers
are separated by $V$. The lowest layer corresponds to a network
similar to the one shown in \fref{fig:System}(b). Due to the structure of
the Hamiltonian transitions take place only between adjacent layers and
layers that are separated by $2V$. The timescale for
inter-layer-transitions is $\sim \Omega^2/V$ while
intra-layer-relaxation takes place on a time $\sim
\Omega^{-1}$. The index $j$ labels the number of adjacent nearest neighbor excitation pairs.} \label{figure3}
\end{figure}

In this work we have investigated the non-equilibrium dynamics and equilibration of a strongly interacting closed quantum system of hard rods of varying size. We showed that the steady state and the non-equilibrium relaxation towards it are well captured by a Master equation whose numerical coefficients can be analytically calculated using the combinatorics of hard rods placed on a lattice. The dependence of the relaxation dynamics on the initial state has been systematically investigated through numerical simulations of the full quantum dynamics. It has been found that relaxation into the steady state is achieved for the majority of initial conditions. Expectation values taken in the steady state are compatible with the assumption of an equal population of each microstate, suggesting that the system here is indeed well described by a microcanonical ensemble.

To conclude, let us finally comment qualitatively on large but finite interaction $V$, i.e. where the system no longer can be described in terms of hard rods. In this more general case we can still employ a representation of the state space in terms of a graph, but as shown in \fref{figure3} (for $\lambda=1$) a third dimension has to be added which accounts for the interaction energy $V$ of different combinations. The lowest layer corresponds to the set of nodes in the $V/\Omega\rightarrow\infty$ limit centered around zero energy with no contiguous excitations, cf. \fref{fig:System}(b). All higher layers contain $j$ pairs of adjacent excitations and hence are centered around the energy $\epsilon_j=j\times V$. Each of these layers can be thought of as an energy shell for which an intra-layer diffusive Fokker-Planck dynamics can be derived.  In reference \cite{leol+:10} it was numerically shown that the steady state population of these energy shells follows a Boltzmann law $p_j=p_{\epsilon_j} \propto e^{-\beta j V}$, with an inverse temperature $\beta$, provided that $V\gg\Omega$. This can be understood as follows: Choosing an initial state in the lowest energy shell the diffusive dynamics establishes a `maximum entropy state' within this shell on a timescale $\Omega^{-1}$. Transitions to higher energy shells take place on a slow timescale $\sim\Omega^2/V$. Single spin flips as effectuated by the Hamiltonian can, however, only couple energy shells with a distance of at most $2V$, i.e. a single spin flip can at most create two pairs of adjacent excitations. That means that the population of $j$-th energy shell can be estimated by $(\Omega^2/V^2)^j$ which approximates a Boltzmann distribution with temperature such that $\beta\Omega\approx -2(\Omega/V)\log (\Omega/V)$. It remains an intriguing question whether also in this more complex system the derivation of a Master equation is possible from first principles.

\ack
We wish to acknowledge discussions with J. Gaskell. C.A. acknowledges support through a Feodor-Lynen Fellowship of the Alexander von Humboldt Foundation. This work was funded in part by EPSRC Grant no. EP/I017828/1 and Leverhulme Trust grant no. F/00114/BG.


\section*{Acknowledgements}

This material is based upon work supported by the National Science Foundation under Grant No. PHY 1212472.


\section{Introduction\label{sec:Introduction}}

Non-linear continuous media can produce finite-time singularities,
where a physical quantity nearly blows up in a fixed amount of time.
In practice, the system never becomes fully singular; new forces become
important at small scales and cut off the singularity. Examples of
physical singularities include black hole formation, supernova explosion,
crack propagation, and fluid pinch-off. Here, we study a singularity
that occurs when an underwater air bubble pinches off from its air
source. Despite the small energies involved, the pinch-off produces
a small region of water moving at speeds of 10s of meters per second
\citep{plesset77}. Many of the interesting features of bubble pinch-off
apply to implosion-type singularities in general.

Bubble pinch-off also has practical applications itself. For example,
\citet{mccormick2009drag} discovered that injecting air bubbles under
the hull of a ship reduces skin drag. Only recently has the idea been
implemented on commercial ships, resulting in an overall 5-10\% fuel
savings \citep{ferrante2004physical,mitsubishi2011}. The drag reduction
depends on the size distribution of the bubbles, which in turn depends
on the specifics of pinch-off \citep{longuet-higgins1991,PhysRevLett.98.144503}.
Also, the bubbles can pinch apart in a turbulent flow, generating
acoustic emissions \citep{frizell1987noise} and dispersing into a
distribution of smaller bubbles, which depends on the details of pinch-off
\citep{Lasheras2002247}. The mechanisms controlling bubble pinch-off
are size-independent and also govern the breach of a submarine hull
\citep{ikeda2009}, and the violent whoosh behind an object plunged
into water \citep{PhysRevLett.96.154505}. This paper focuses specifically
on the effects of up--down asymmetry in bubble dynamics. Up--down
asymmetry is present in nearly all bubbles due to the hydrostatic
pressure gradient. In general, the bottom portion of the neck is more
slender than the top due to the hydrostatic pressure. This effect
can be exacerbated if the bubbles are generated in a dynamic stream
of water \citep{oguz93}. In fact, this research has potential implications
for the jet of water that shoots upward after such an impact, called
the Worthington jet after \citet{worthington1897}, \citep{CambridgeJournals:7921450}.
\S\ref{sec:Conclusion} discusses how the initial conditions of the
cavity determine the bubble shape just before the impact that spawns
the Worthington jet.

\begin{figure}
\begin{centering}
\includegraphics[clip]{long1}
\par\end{centering}

\caption{\label{fig:experimental_images}Bubble (dark) pinching off from a
nozzle submerged in water (light). Bright spots are optical artefacts.
(a) Initially, the bubble neck shape near the minimum has a generic
quadratic profile. Zooming in (b-c), just before pinch-off the shape
becomes two cones connected at the vertex by a short segment. A characteristic
vertical length scale $L_{z}$ for the neck is the distance from the
minimum to a height where the neck radius is $2R_{\mathrm{min}}$
(b). Even at microscopic scales (c) the neck is relatively squat $(L_{z}/R_{\mathrm{min}}\approx3.2,$
$\theta^{+}\approx28^{\circ}),$ still far from the theoretical end-state
$(L_{z}/R_{\mathrm{min}}\rightarrow\infty$ and $\theta\rightarrow0^{\circ}).$
Neck vibrations eventually cut off the dynamics before the singularity
is reached. (Images courtesy of N. C. Keim and S. R. Nagel).}

\end{figure}

Few people realise that they are creating a singularity when they
blow a bubble into a glass of water using a drinking straw. In experimental
set-ups, the effect can be controlled by precisely pumping air into
a pressurised nozzle (see figure \ref{fig:experimental_images}a).
As the air rises, surface tension competes with buoyancy to create
a quadratic neck shape. Then, the neck becomes unstable and water
rushes in from all sides to close it off. This sets up the regime
we study in this paper, where the inertial flows diverge in speed
where they are closing the shrinking neck, dominating all other forces
and creating a singularity. If this regime were allowed to continue
unaltered, the neck would approach a long-and-slender cylinder, with
purely radial flows \citep{longuet-higgins1991}. Before this occurs,
however, fast air-flows develop in the neck, potentially altering
the collapse scaling \citep{PhysRevLett.95.194501,PhysRevLett.101.214502}.
Eventually, the air even goes supersonic \citep{PhysRevLett.104.024501}.
Just before the airflow effect, some authors including \citet{PhysRevLett.95.194501}
and \citet{PhysRevE.83.056325} measure a small residual asymmetry
that influences the air flow. This is consistent with our results
that the initial asymmetry decays as a power law, but is still measurably
appreciable after at this scale. Finally, \citet{PhysRevLett.97.144503}
found that if the bubble is not quite axisymmetric to begin with,
the deviation actually becomes more pronounced as the neck approaches
pinch-off. The outcome is a side-to-side contact before the average
radius reaches $0,$ cutting off the singularity. These vibrations
were confirmed and studied in detail in later theory \citep{schmidt2009},
experiments \citep{PhysRevE.83.056325,enriquez2011}, and simulations
\citep{Turitsyn2009,lai:102106}.

This singularity is unusual in that it reaches the cut-off point before
the neck shape has a chance to reach its predicted end-state shape,
which is a long cylinder (connecting slender cones on the top and
bottom, see \S2). Most singularities with a unique end-state shape
(\textit{universal} singularities) rapidly approach that shape with
power-law dependence \citep{shi1994}. In those cases, the transient
effects leading up to the universal shape are short-lived and ignorable.
In contrast, the air bubble neck shape evolves very gradually as a
function of the neck radius $R_{\mathrm{min}},$ only going as $\log\left(R_{\mathrm{min}}\right)$\citep{PhysRevLett.95.194501,egg07,PhysRevE.80.036305}.
Indeed, the transition is so slow that on observable length scales
the neck never appears cylindrical, but rather as two relatively squat
cones, smoothly connected (see figure \ref{fig:experimental_images}b-c).
Although technically universal, this singularity spends its lifetime
in the transient phase. Therefore, this paper explores the various
transients that can occur in evolution.

We begin by exploring the effect of highly distorted initial neck
shapes, including necks with large curvature and up--down asymmetry.
We found that the asymmetry is truly transient in every sense: it
vanishes after only several factors of $10$ decrease in the minimum
radius. The neck undergoes this transition not by flexing the cones,
but by shifting the vertical position of the pinch-off. By shifting
instead of readjusting in place, the (new) collapse region rapidly
becomes nearly up--down symmetric. Then, the neck goes through a phase
of being symmetric but also squat, with significant (symmetric) vertical
flows \citep{PhysRevE.84.026313}. Depending on where the neck shifts
vertically, the symmetric squat neck inherits its dynamics from either
its top or bottom side, or a blend of both. Interestingly, we found
that by tuning the initial vertical flow near the neck, we could adjust
the degree of this blending. This effect may explain some of the late-stage
effects seen in experiments that tune the density of the inner gas
\citep{PhysRevLett.101.214502,PhysRevE.83.056325}. When the bubble
is in the formation stage, gas rises through the neck. Denser gasses
impart more momentum to the neck area, giving the neck an initial
vertical velocity component. We predict that this should cause the
neck to inherit its aspect ratio from the top side, which is more
squat. \citet{PhysRevLett.101.214502} do in fact see more squat end-state
neck shapes for higher gas densities (up to their observed critical
density). \citet{PhysRevE.83.056325} observes that the very early
stage air flow during formation sets up a water velocity field, which
then evolves independent of the air during the inertial stage studied
in this paper. He also confirms our observations that the water flow
affects the up--down symmetry and neck position later in the evolution.
Neither experiment, however, attempted to manipulate the initial neck
velocity field purposely, except by altering the gas density. Given
our simulation results, it should be worthwhile to explore more deliberate
ways of controlling the initial airflow through the neck, for example
by blowing the bubble with a top and bottom nozzle simultaneously.

\S\ref{sec:background} summarises the previous research on the inertial
stage of the bubble singularity, and the need for looking at up--down
asymmetry separately. \S\ref{sec:methods} gives an overview of the
computational methods used to perform the studies. The results are
divided into two sections for ease of reading: \S\ref{sec:ResultsI_shape_asymmetry}
gives a full description of a typical collapse. \S\ref{sec:Results-part-II}
describes how the initial flow field, specifically the vertical velocity
gradient, can be used to tune the details of the final collapse. \S\ref{sec:Conclusion}
concludes with a discussion of how the results could be seen in experiment,
and problems that remain unsolved.

\section{\label{sec:background}Background}

The first stage, described in \S\ref{sec:Introduction}, sets up
the initial conditions for the implosion, and involves the interplay
between surface tension and buoyancy. After this point, surface tension
becomes dwarfed by inertial forces. This is the point where we begin
our simulations, and we assume inviscid flow with only inertial and
pressure forces.

We use the common assumption that the water's velocity field $\mathbf{u}$
is incompressible $(\boldsymbol{\nabla\cdot u}=0),$ irrotational
$(\boldsymbol{\nabla}\times\mathbf{u}=\mathbf{0}),$ and decays to
zero far away.

We add another assumption that the air in the bubble is dynamically
passive, with a uniform pressure $P_{\mathrm{air}}(t)$ whose value
ensures constant bubble volume in time. Experiments show that after
the neck begins to collapse, airflow is unimportant until the neck
radius reaches roughly 100 \textmu{}m \citep{PhysRevLett.101.214502,PhysRevE.83.056325}.
Around that point, airflow through the neck becomes relevant and our
assumption breaks down. Soon after, neck vibrations dominate. Here,
we will only study what happens up to this point, assuming an axisymmetric
bubble with uniform air pressure.

Being curl-free, the velocity can be described by a scalar potential,
defined as $\boldsymbol{\nabla}\phi\equiv\mathbf{u}.$ Since the exterior
flow is incompressible, $\boldsymbol{\nabla\cdot}\left(\boldsymbol{\nabla}\phi\right)=0$
implies \begin{equation}
\mathbf{\nabla}^{2}\phi=0.\label{eq:Laplace_equation}\end{equation}
That is, Laplace's equation governs the potential in the bulk. This
allows $\mathbf{u}$ to be solved everywhere in the exterior given
$\phi$ on the interface. Equating the relevant stresses on the interface
gives a non-linear differential equation for $\phi$ on the surface
of the bubble, first-order in time and space: \begin{equation}
P_{\mathrm{air}}+\rho\left(\frac{\partial\phi}{\partial t}+\frac{1}{2}\left|\boldsymbol{\nabla}\phi\right|^{2}\right)=0.\label{eq:stressboundary}\end{equation}

where $\rho$ is the water density. The previous two equations determine
the evolution of the velocity field. The velocity field in turn determines
how the interface evolves. Surface elements are advected by the local
velocity on the surface, $S:$

\begin{equation}
\left.\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}\right|_{S}=\left.\boldsymbol{\nabla}\phi\right|_{S}.\label{eq:kinematic_condition}\end{equation}

Previously, the solutions to equation \ref{eq:stressboundary} were
only studied carefully in the up--down symmetric, long-and-slender
limit \citep{PhysRevLett.95.194501,egg07,PhysRevE.80.036305}. These
solutions produce elegant, analytic results showing the neck approaching
a cylinder logarithmically slowly as $R_{\mathrm{min}}\rightarrow0$
. There are differences between the theories, but each assumes a (time-dependent)
purely radial flow. This is equivalent to saying that the potential
distribution is proportional to the zeroth order cylindrical harmonic:

\begin{eqnarray}
\phi_{cyl}\left(t\right) & = & -q\left(t\right)\log\left[r/r_{\infty}\left(t\right)\right],\nonumber \\
\mathbf{u}_{cyl}\left(t\right) & = & -q\left(t\right)\boldsymbol{\hat{r}}/r,\label{eq:u_cyl}\end{eqnarray}
 where $r_{\infty}\left(t\right)$ is a length scale related to the
system size, and $q\left(t\right)$ is related to the rate of inflow,
and is determined by equation \ref{eq:stressboundary}. For a perfectly
cylindrical neck with constant $P_{\mathrm{air}}$ and $r_{\infty},$
the reader may verify that $q\left(t\right)=\left[\left(P_{\mathrm{air}}R_{\mathrm{min}}^{2}+\mathrm{const}\right)/\log\left(r_{\infty}/R_{\mathrm{min}}\right)\right]^{\frac{1}{2}}$
satisfies the governing equations. In most cases, the pressure term
is negligible, leaving $q\left(t\right)\sim\left[\log\left(r_{\infty}/R_{\mathrm{min}}\right)\right]^{-\frac{1}{2}},$
which is nearly constant. This defines a flow, $\boldsymbol{u}=-\left[\log\left(r_{\infty}/R_{\mathrm{min}}\right)\right]^{-\frac{1}{2}}r^{-1}\boldsymbol{\hat{r}},$
that evolves the neck gradually towards a cylinder. If $q(t)$ were
instead constant, \citet{longuet-higgins1991} found that any parabola-of-rotation
would generically evolve into a hyperbola-of-rotation. Indeed, at
late stages the bubble neck shape does resemble a hyperbola-of-rotation
by the fact that the neck looks like two smoothly connected cones
(figure \ref{fig:experimental_images}c), which we classify by top
and bottom opening angles, $\theta^{+}$ and $\theta^{-}$ respectively.
For consistency, we actually measure an effective cone angle $\theta_{\mathrm{eff}}$
defined by $2R_{\mathrm{min}}=R_{\mathrm{min}}\sqrt{1+\left(\tan^{2}\theta_{\mathrm{eff}}\right)\left(L_{z}/R_{\mathrm{min}}\right)^{2}}.$
Taking into account the logarithmic dependence of $q$ on $R_{\mathrm{min}},$
the cone opening angles decrease very slowly to $0,$ approaching
a cylinder. Even at the small scale of figure \ref{fig:experimental_images}c,
the cone angle is approximately $28^{\circ}.$ Even under ideal conditions
the neck does not become cylindrical before neck vibrations tear off
the bubble.

\begin{figure}
\centering{}\includegraphics{long5}\caption{\label{fig:cone_angle_evolution_vs_cyl}The near-cylindrical theory
of \citet{PhysRevE.80.036305} (dotted lines) fails to properly reproduce
the simulation (solid lines, see \S\ref{sec:methods}) for a neck
shape that is initially up--down asymmetric and squat $(\theta^{+}=40^{\circ}$
and $\theta^{-}=5^{\circ}),$ with initial radial flow. The simulation
reveals that the asymmetry is a short-lived transient, whereas the
theory predicts that the top and bottom sides evolve slowly and separately
to $0^{\circ}$ cone angle.}

\end{figure}

There are two problems with analytic solutions of the form in equation
\ref{eq:u_cyl}. First, they only apply when the neck shape is long-and-slender
and up--down symmetric. Neither assumption holds for real bubble necks
(see figure \ref{fig:experimental_images}). The theory predicts qualitatively
incorrect behaviour for initially up--down asymmetric necks (figure
\ref{fig:cone_angle_evolution_vs_cyl}). In the theory, the neck shape
evolves as a stack of non-interacting cross sections, each one governed
by equation \ref{eq:u_cyl} (or similar). In reality, vertical interaction
plays a large role in the shape evolution. Without vertical flow,
the asymmetry decays logarithmically slowly. With vertical flow, we
find that the asymmetry vanishes as a power law (measured to go as
$R_{\mathrm{min}}^{-1}).$

\begin{figure}
\begin{centering}
\includegraphics{long6}
\par\end{centering}

\caption{\label{fig:Effective-cone-angle}Only when the top and bottom sides
fully interact does our simulation show the up--down asymmetry to
be short-lived (solid lines). For comparison, we also evolved two
symmetric shapes with $\theta^{+}=\theta^{-}=40^{\circ}$ and $\theta^{+}=\theta^{-}=5^{\circ},$
respectively (dashed lines). Taken together, these two simulations
show how the asymmetric shape would evolve if the vertical position
of the neck were fixed (and no flow were allowed through the $z=0$
plane). In that case, the asymmetry decreases very gradually. This
comparison shows that the asymmetry vanishes quickly because the top
and bottom interact. Specifically, the vertical position of the neck
shifts. }

\end{figure}

In fact, for up--down asymmetric shapes, the most important effect
comes from the global interaction between the top and bottom sides.
If we artificially forbid flow through the mid-plane, $z=0,$ the
up--down asymmetry no longer vanishes quickly (figure \ref{fig:Effective-cone-angle}
dashed lines). \S\ref{sec:ResultsI_shape_asymmetry} describes the
mechanism in detail. In short, the singularity shifts vertically in
the direction of the smaller cone angle, and in doing so equalises
the cone angles. \S\ref{sec:Results-part-II} shows that this effect
can be manipulated. By introducing a vertical flow to the initial
conditions, the vertical position of the neck can be tuned. In doing
so, the resulting symmetric neck inherits a different blend of characteristics
from the top and bottom sides. This effect should be realisable in
experiments by adjusting the air flow rate through the neck during
the initial set-up phase.

\section{\label{sec:methods}Methods}

To Solve equations \ref{eq:Laplace_equation}-\ref{eq:kinematic_condition},
we use the standard procedure of reformulating the equations as a
boundary value problem, where the bulk differential equation becomes
a surface integral equation. This requires two steps. First, the ordinary
derivative in equation \ref{eq:stressboundary} is replaced by the
co-moving derivative $\frac{D\phi}{Dt}=\frac{\partial\phi}{\partial t}+\boldsymbol{u\cdot\nabla}\phi,$
and will apply to the fluid parcels that constitute the surface. Second,
equation \ref{eq:Laplace_equation} is reformulated using Green's
integral:

\begin{equation}
2\pi\phi\left(\mathbf{x}_{0}\right)=\oint\limits _{\mathbf{x}_{1}\in S}\left(\frac{\left(\mathbf{x}_{1}-\mathbf{x}_{0}\right)\phi\left(\mathbf{x}_{1}\right)}{\left|\mathbf{x}_{1}-\mathbf{x}_{0}\right|^{3}}+\frac{\boldsymbol{\nabla}\phi\left(\mathbf{x}_{1}\right)}{\left|\mathbf{x}_{1}-\mathbf{x}_{0}\right|}\right)\boldsymbol{\cdot}\mathbf{\hat{n}}\mathrm{d}S,\label{eq:green}\end{equation}
 where $\mathbf{x}_{0}$ is a fixed point on $S,$ $\mathbf{x}_{1}$
parameterizes the surface, and $\mathbf{\hat{n}}$ is the unit surface
normal at $\mathbf{x}_{1}$ pointing into the bubble.

With both governing equations defined on the interface, the formulation
reduces to 2 dimensions. We further reduce to 1 dimension by assuming
axisymmetry. At each time step, we use cubic splines to interpolate
the $N(t)$ discrete surface nodes $\mathbf{x}_{i}$. Then, for each
$i\in\left\{ 1\ldots N\right\} $ we allow $\mathbf{x}_{i}$ to take
the place of $\mathbf{x}_{0}$ in equation (\ref{eq:green}). We use
Gaussian quadrature to perform the integral over segments connecting
adjacent spline midpoints $j\in\left\{ 1\ldots N\right\} ,$ with
the first and last segments being only half splines. Thus, the $N$
instances of equation \ref{eq:green} can be schematically written
as $\phi_{i}=\boldsymbol{\mathsf{\mathit{A}}}_{ij}\phi_{j}+\boldsymbol{\mathsf{\mathit{B}}}_{ij}\left(u_{\perp}\right)_{j},$
where $\boldsymbol{\mathsf{\mathit{A}}}$ and $\boldsymbol{\mathsf{\mathit{B}}}$
are $N$-by-$N$ matrices and $u_{\perp}\equiv\boldsymbol{\nabla}\phi\boldsymbol{\cdot}\mathbf{\hat{n}}.$
Therefore, given $\phi$ on the surface we have $u_{\perp}=\boldsymbol{\mathsf{\mathit{B}}}_{ij}^{-1}\left(\phi_{i}-\boldsymbol{\mathsf{\mathit{A}}}_{ij}\phi_{j}\right)$
and $u_{\parallel}=\frac{\mathrm{d}\phi}{\mathrm{d}s},$ where $\mathrm{d}s$
is an infinitesimal arc length. Using $\Delta t\sim R_{\mathrm{min}}/\dot{R}_{\mathrm{min}},$
we advect the N nodes using $\mathbf{u}$ and evolve $\phi$ using
equation \ref{eq:stressboundary}, completing the cycle \citep{pozrikidis,oguz93}.
The pressure $P_{\mathrm{air}}$ can be solved for explicitly after
evolving the shape and before updating $\phi_{i}$ (formula and derivation
available upon request).

\begin{figure}
\begin{centering}
\includegraphics{long2}
\par\end{centering}

\caption{\label{fig:shape_evolution} (a) The initial surface is comprised
of an upper half with opening angle $\theta^{+}=40^{\circ},$ and
a lower half with opening angle $\theta^{-}=5^{\circ}.$ The initial
flow is radial. (b) Zooming in $50\times,$ four snapshots show the
neck becoming rapidly more symmetric via a descent of the neck position
(indicated by black curves). In this egregiously asymmetric example,
a second neck actually forms causing a jump in the vertical neck position
(seen as a tiny kink in the black curve in b). (c) Zooming in another
$50\times,$ the simulation evolves the bubble to length scales smaller
than physically possible. Here, the shape is nearly cylindrical and
the neck shift is insignificant. Time steps between curves are: (b)
$0.19\, R_{0}/u_{0r},$ and (c) $4.4\cdot10^{-4}\, R_{0}/u_{0r}.$
Each innermost profile becomes the outermost in the next image. Every
tenth node is shown. }

\end{figure}

In order to accurately resolve the pinch-off dynamics, we found it
important to use a node distribution scheme that maintains a gradual
variation in the spacing between node points and that continually
adds nodes in the neighbourhood of the minimum (see figure \ref{fig:shape_evolution}).
At each moment, if the vertical distance $\Delta z$ between nodes
exceeds a maximum spacing $\Delta z_{\mathrm{max}}(t,\, z)$ in the
region where the neck radius is less than $4R_{\mathrm{min}},$ we
add a new node at the midpoint of that spline. After experimenting
with several functions for $\Delta z_{\mathrm{max}},$ we found $\Delta z_{\mathrm{max}}=\max\left[0.15\min\left(\left|z\right|,L_{z}\right),0.08L_{z}\right]$
allows the simulation to accurately track the dynamics %
\footnote{For the up--down asymmetric case, we define the current height of
the neck minimum as $z_{0},$ $R\left(z_{0}\pm L_{z}^{\pm}\right)=2R_{\mathrm{min}},$
and $\Delta z_{\mathrm{max}}=\max\left[0.15\min\left(\left|z-z_{0}\right|,L_{z}^{+}\right),0.08\min\left(L_{z}^{+},L_{z}^{-}\right)\right]$
for $z\geq z_{0}$ (reverse superscripts for $z<z_{0}).$ %
}.

We begin the simulation by specifying an interface shape and velocity
field. In general there are an infinite number of each. Here, we examine
generic smooth distributions. To expedite the computation, we do not
prescribe a purely quadratic shape profile at $t=0,$ but instead
use the result from previous studies that a slender quadratic neck
evolves into a hyperbola, and prescribe:

\[
R(z)=\sqrt{\left(R_{0}^{2}+z^{2}\tan^{2}\theta\right)-\left(z^{4}/h^{4}\right)\left(R_{0}^{2}+h^{2}\tan^{2}\theta\right)},\]
where $2\theta$ is the opening angle of the cone, $R_{0}\equiv R_{\mathrm{min}}\left(t=0\right),$
and $h=140R_{0}$ is height from the neck minimum to the end of the
bubble.

To initialise the velocity field, most previous research uses a pure
radial collapse, because it is the counterpart of the cylindrical
end-state solution. We use this as a starting point. Equivalently,
we prescribe a line of point sinks along the z-axis:

\begin{eqnarray*}
\Phi_{l.s.}\left(r,z\right) & \equiv & \intop_{-L}^{L}\left(r^{2}+\left(z-\zeta\right)^{2}\right)^{-\frac{1}{2}}\mathrm{d}\zeta\\
 & = & \ln\left[\frac{z+L+\sqrt{r^{2}+\left(z+L\right)^{2}}}{z-L+\sqrt{r^{2}+\left(z-L\right)^{2}}}\right].\end{eqnarray*}

We use $L=100,$ much larger than the neck region, to avoid end effects.
For \S\ref{sec:Results-part-II}, we also add a vertical velocity
component to the initial conditions of the neck region. One trivial
possibility is a uniform vertical component. That would be the same,
however, as just boosting the neck region to a new reference frame
$z=u_{0z}t,$ and the pinch-off would be identical, only boosted.
Instead, we use the following form, which consists of a line of sinks
(sources) along the z-axis that vary in strength proportional to height
$\zeta:$

\begin{eqnarray}
\Phi_{z}\left(r,z\right) & \equiv & \intop_{-L}^{L}\zeta\left(r^{2}+\left(z-\zeta\right)^{2}\right)^{-\frac{1}{2}}\mathrm{d}\zeta\nonumber \\
 & = & d_{-}-d_{+}+z\log\left[L-z+d_{-}\right]-z\log\left[-L-z+d_{+}\right],\label{eq:phi_z}\end{eqnarray}

where $d_{\pm}\equiv\sqrt{r^{2}+\left(L\pm z\right)^{2}},$ and the
limit $r\rightarrow0$ is taken for points on the z-axis. This field
acts as a stretched out dipole. The total initial potential field
consists of a weighted sum of the line sink and the dipole field:

\begin{equation}
\Phi_{tot}=\left(\frac{R_{0}d_{0}}{2L}\right)u_{0r}\Phi_{l.s.}+\left(-\frac{2L}{d_{0}}+\log\left[\frac{d_{0}+L}{d_{0}-L}\right]\right)u_{0z}\Phi_{z},\label{eq:phi_line_sink}\end{equation}
where $d_{0}\equiv\sqrt{R_{0}^{2}+L^{2}}.$ The terms in parentheses
are normalisation factors, so that the initial velocity at the neck
minimum $(z=0)$ has radial and vertical components $u_{0r}$ and
$u_{0z},$ respectively. This initial velocity field has the added
benefit of resembling the effect of vertical airflow through the neck
on the surrounding water.

\section{Results\label{sec:Results}}

Our simulations show that the evolution of the neck towards singularity
depends critically on both the initial neck shape and the initial
velocity field. We show the two effects separately.

\S\ref{sec:ResultsI_shape_asymmetry} details how an up--down asymmetric
shape evolves with a simple initial velocity field: pure radial flow.
The velocity field immediately develops a pressure ring just outside
the neck region, but below the neck minimum itself (in general, towards
the side with smaller cone angle). As the pressure peak moves inward,
it bulges the shape inward, shifting the vertical location of the
neck minimum and making the neck symmetric. The symmetric neck then
progresses as a symmetric collapse, with the upper cone angle adjusted
to match the bottom.

\S\ref{sec:Results-part-II} adds the effect of initial vertical
flow around the neck. A small upward flow is enough to oppose the
neck from shifting downward. After becoming symmetric, the neck is
a blend of the top and bottom cone angles. This gives a possible experimental
method for adjusting the cone angle of the neck at pinch-off.

\begin{figure}
\centering{}\includegraphics{long3}\caption{\label{fig:pressure}The characteristic pressure peak that drives
the singularity develops below neck minimum of the initial shape (original
neck shape: dark dashed curves; original neck minimum: thin horizontal
dotted line), (left). The pressure peak lies just outside the new
neck minimum, with steep gradients between the peak and the surface,
along which the pressure is fixed to be $P_{\mathrm{air}}\left(t\right),$
(right). (Note: this point of evolution corresponds to the centre-most
profile in figure \ref{fig:shape_evolution}b).}

\end{figure}

\subsection{Shape Asymmetry\label{sec:ResultsI_shape_asymmetry}}

We studied many different combinations of initial top and bottom cone
angles $(\theta_{+}$ and $\theta_{-}$ respectively), with each angle
ranging between $5^{\circ}$ and $60^{\circ},$ along with different
initial velocity distributions. The results are qualitatively the
same for any combination of top and bottom cone angles, but the asymmetric
effects become more pronounced for larger asymmetry. In this section,
we use a neutral radial inflow (equation \ref{eq:phi_line_sink} with
$u_{0z}=0$), and an initial $40^{\circ}$ top cone and $5^{\circ}$
bottom cone.

The initial shape consists of the neck, a hyperbola-of-rotation with
the given cone angles, rounded off by large end-caps (figure \ref{fig:shape_evolution}a).
Immediately after the simulation starts, a ring of high pressure develops
in the water surrounding the neck, but slightly below the initial
neck minimum. This causes the shape to bulge in at a height below
the original neck minimum (figure \ref{fig:shape_evolution}b), and
the position of the neck minimum to descend. In our extreme example,
the descent occurs so fast that a new neck actually forms, which appears
as a jump in the neck's vertical position (figure \ref{fig:shape_evolution}b).
Meanwhile, the pressure peak continues to grow in strength and become
more symmetric (figure \ref{fig:pressure}). The neck minimum shifts
downward at roughly a constant rate (lab frame) until pinch-off (figure
\ref{fig:shape_evolution}b). This constant vertical rate is soon
overcome by the large radial implosion flow, and becomes insignificant
at later times (figure \ref{fig:shape_evolution}c). By this time,
the neighbourhood of the neck is completely symmetric, and the system
undergoes the previously-studied symmetric collapse. In this process,
the top side of the shape performs most of the accommodation to match
the bottom side, and the large initial upper angle has little effect
after the transient (figure \ref{fig:Effective-cone-angle}).

\begin{figure}
\centering{}\includegraphics{long7}\caption{\label{fig:vertical_velocity_map}(Left) The pressure peak from figure
\ref{fig:pressure} creates a symmetric stretching flow centred below
the original neck minimum, with positive vertical velocity component
$\boldsymbol{u\cdot\hat{z}}$ above the peak and negative below. This
flow leads the shape into a symmetric, vertically-offset collapse.
(Right) Initially, the vertical flow is prescribed with no vertical
flow near the neck minimum (innermost 1). Immediately, the pressure
peak (figure \ref{fig:pressure}) develops, mainly generating positive
vertical velocities on the top side (2). The pressure peak shifts
downward, forming a new, symmetric neck, below the original minimum
(3). Out of this neck, an up--down symmetric straining flow develops
(4-7). Curve 5 (thick) corresponds to the left figure. When the straining
flow begins to develop, it first is only a local maximum on the top
surface ({*}), but continues to grow (6-7), eclipsing the original
peak on the top side, which stagnates and becomes irrelevant. The
straining flow is slightly offset in $u_{z},$ signifying that the
neck continues shifting at a constant rate in the lab frame.}

\end{figure}

The process can be seen in more detail by studying the evolution of
the velocity field, specifically the vertical velocity. Initially,
no vertical flow is present (figure \ref{fig:vertical_velocity_map}
curve 1), but immediately an asymmetric straining flow develops (figure
\ref{fig:vertical_velocity_map} curve 2). Soon afterwords, the origin
of the straining flow shifts downward (figure \ref{fig:vertical_velocity_map}
curve 3). Driven by the pressure ring, the straining flow becomes
symmetric (figure \ref{fig:vertical_velocity_map} curve 4), which
in turn makes the neck shape symmetric. Eventually, the vertical velocity
forms a local maximum in $u_{z}$ on the upper side (figure \ref{fig:vertical_velocity_map}
curve 5, indicated by '{*}'). From then on, the straining flow grows
asymptotically strong, leaving behind only vestiges of the original
top side, which stagnate and become insignificant (figure \ref{fig:vertical_velocity_map}
curves 6-7). The transition from curve 2 to curve 4 shows that the
top part of the straining flow adjusts to match the bottom half through
a downward shift of the straining flow centre.

\begin{figure}
\centering{}\includegraphics{long8}\caption{\label{fig:phase_diagram} The two sides of the asymmetric evolution
(solid curves) evolve drastically differently from their symmetric
counterparts when viewed in phase space (velocity ratio $\left(\max\left|u_{z}\right|/\mathrm{max}\left|u_{r}\right|\right),$
versus the neck aspect ratio $L_{z}/R_{\mathrm{min}}).$ Due to the
strong asymmetry, the two sides see each other and their phase trajectories
attract. The $40^{\circ}$ side is stretched as the neck rolls down,
while the $5^{\circ}$ side is squashed. Meanwhile, the $40^{\circ}$
side is prevented from reaching its full potential vertical velocity.
This is because the singular straining flow is generated on the $5^{\circ}$
side (see figure \ref{fig:vertical_velocity_map}). The newly formed
local maximum vertical velocity (figure \ref{fig:vertical_velocity_map},
'{*}') is represented here also by a '{*}'. This peak eventually prevails
as the global maximum. The two phase curves finally meet, as the neck
region becomes symmetric, then slowly becomes more slender going towards
pinch-off.}

\end{figure}

The shape and velocity field are interdependent, and it is informative
to show how they co-evolve. Figure \ref{fig:phase_diagram} condenses
this information into a phase space. The vertical axis gives the ratio
of the maximum surface vertical velocity component to the maximum
surface radial velocity component (two distinct points). The horizontal
axis gives the shape aspect ratio $L_{z}/R_{\mathrm{min}},$ i.e.,
when the trajectory moves to the right, it corresponds to the neck
shape becoming more long-and-slender. The top and bottom sides each
have separate trajectories, and for comparison, we show each side's
symmetric counterpart (which shows what would happen if no flow were
allowed through the plane $z=0).$ Each trajectory starts on the horizontal
axis, corresponding to the prescribed zero initial vertical velocity,
at a position depending on the initial shape. For example, the $40^{\circ}$
initialised side begins with a corresponding aspect ratio of $L_{z}/R_{0}=\sqrt{3}\cot40^{\circ}=2.06.$
After an initial transient, all the curves coalesce onto a universal
curve \citep{PhysRevE.84.026313} and slowly progress to the universal
form: long-and-slender (large shape aspect ratio) and radial flow
(zero vertical velocity). The progression along the universal curve
has been studied. We focus here on the transient.

Each side of the asymmetric shape differs significantly from its symmetric
counterpart, indicating that the vertical flow through the mid-plane
(interaction between the top and bottom) $z=0$ is important. The
flow through the mid-plane allows the neck to shift until the straining
flow is symmetric. The result is that the trajectories of the top
and bottom sides {}``attract'' one another in phase space. In this
case, the direction of the neck shift is downward. The downward flow
of the neck shift opposes the straining flow on the $40^{\circ}$-initial
top side, and augments the straining flow on the $5^{\circ}$-initial
bottom side. The result is that the peak surface vertical velocity
ratio undershoots its symmetric counterpart on the top side (figure
\ref{fig:phase_diagram} top solid vs. dashed), whereas it overshoots
its symmetric counterpart on the bottom side (figure \ref{fig:phase_diagram}
bottom solid vs. dashed). Similarly, the squat top side shape becomes
more slender and the slender bottom side becomes more squat. In this
way, the two trajectories are brought together in a short time. Some
important features seen in previous plots can also be seen in figure
\ref{fig:phase_diagram}. The local maximum vertical velocity ('{*}'
in figure \ref{fig:vertical_velocity_map}) is shown as a separate
trajectory in figure \ref{fig:phase_diagram} (also indicated by '{*}').
Also, figure \ref{fig:phase_diagram} shows that the bottom side does
not have to detour much in order to meet up with the top side trajectory.
This effect was also seen in figure \ref{fig:Effective-cone-angle}.
In short, since the neck shifts down in this case, the singularity
mostly inherits the characteristics of the bottom side. It is the
top side that adjusts to mirror the bottom side, and the singularity
ignores the fact that the top was initially squat. This effect can
be reversed, as seen in the next section.

\subsection{Tunable cone angle\label{sec:Results-part-II}}

\begin{figure}
\centering{}\includegraphics{long11}\caption{\label{fig:cone_angle_evolution_with_uz}Adding some vertical flow
$\left(u_{0z}=0.4u_{0r}\right)$ to the initial conditions causes
a significant change later in the evolution. Instead of following
the trajectory of a $5^{\circ}$-initial shape (see figure \ref{fig:Effective-cone-angle}),
the evolution more closely follows the $40^{\circ}$-initial shape. }

\end{figure}

To understand the context for this section, refer back to figure \ref{fig:Effective-cone-angle}.
There, after a short transient, the top and bottom sides (solid lines)
converge to a single trajectory that closely follows the trajectory
of a $5^{\circ}$--symmetric initial shape (lower dashed line). Here,
by adding a vertical velocity to the initial conditions (equation
\ref{eq:phi_line_sink} with non-zero $u_{0z}$), the picture changes.
A prescribed initial vertical velocity 40\% in strength to the initial
radial velocity, $u_{0z}=0.4u_{0r},$ causes the shape to evolve into
a symmetric shape that more closely follows the top-side $40^{\circ}$--symmetric
initial trajectory (figure \ref{fig:cone_angle_evolution_with_uz},
compare to figure \ref{fig:Effective-cone-angle}). The added vertical
velocity opposes the tendency for the neck to shift down in the early
evolution. As a result, when the neck becomes symmetric, it inherits
the characteristics of the top side, rather than the bottom side.
This modest change in initial condition causes a qualitatively different
outcome.

\begin{figure}
\centering{}\includegraphics{long10}\caption{\label{fig:uz_effects}For neutral and downward initial neck velocities,
the resulting shape after the transient (fixed here as $R_{\mathrm{min}}=0.001R_{0})$
is the same as if the initial shape were $5^{\circ}$--$5^{\circ}$
symmetric (lower line). That is, the upper, $40^{\circ}$ side has
no effect on the late-stage dynamics. For upward initial neck velocities,
i.e. $u_{0z}>0,$ a new effect occurs. The upward flow disrupts the
downward neck shift, and the post-transient neck shape is a blend
of the initial top and bottom sides. Evidently, a modest initial vertical
flow around the neck can tune the properties of the end-state singularity.}

\end{figure}

The effect is not bimodal, but continuous (see figure \ref{fig:uz_effects}).
For $u_{0z}/u_{0r}$ between the range of about $0.0$ to $0.5,$
the resulting trajectory smoothly varies between the bottom-side dynamics
and the top-side dynamics. For initial downward flow $(u_{0z}<0.0),$
the system is insensitive to the exact value; any downward velocity
just results in the singularity inheriting the bottom-side dynamics
(there is no overshoot effect). Interestingly, the dynamics never
fully approach the top-side trajectory, but seem to asymptote slightly
below it. We do not have an explanation for this disparity.

\begin{figure}
\centering{}\includegraphics{long9}\caption{\label{fig:phase_diagram_uz}Adding a small vertical flow to the initial
conditions has a big consequence on the transient. Here, the initial
shape is the same $\theta_{+}=40^{\circ},$ $\theta_{-}=5^{\circ},$
but now there is an added vertical flow given by equation \ref{eq:phi_line_sink}
with $u_{0z}/u_{0r}=0.4.$ The main effect is that the point where
the two solid lines meet (the point where the shape becomes symmetric)
is significantly further left (corresponding to a larger effective
cone angle). The upward initial velocity opposes the downward neck
shift, and the resulting singularity retains more of the squatness
of the upper side. This effect is tunable by adjusting $u_{0z}/u_{0r}.$ }

\end{figure}

With an added vertical velocity $u_{0z}/u_{0r}=0.4,$ it is the top
side trajectory that dictates the evolution, with the bottom side
going out of its way to accommodate (figure \ref{fig:phase_diagram_uz}).
The point on the phase diagram where the two trajectories merge (becoming
symmetric) is shifted to the left and up, corresponding to a more
squat aspect ratio with larger straining flow. Without the downward
neck shift, the top side evolves almost as if it does not see the
bottom (figure \ref{fig:phase_diagram_uz}, similarity between top
solid and top dashed). The initial vertical flow keeps the neck position
elevated, and the singularity develops within the context of the top
side. The resulting singularity inherits its characteristics from
the top side, nearly ignoring the bottom side completely. A small
initial vertical velocity near the neck can change whether the singularity
inherits from the top side, the bottom side, or a blend of each.

\section{\label{sec:Conclusion}Conclusion}

In general, bubble pinch-off begins with some up--down asymmetry.
We found that the neck adjusts the vertical position of the pinch-off
point to make the neck shape and velocity field symmetric. If the
initial flow field around the neck is purely radial, then the neck
shifts (down) towards the more slender cone. The new neck forms with
a top that emulates the bottom. The resulting singularity depends
only on the initial $\theta^{-}.$ The situation is reversed if a
strong upward velocity is introduced to the neck during the creation
stage. Then, the vertical flow prevents the neck from descending,
and the new neck depends only on $\theta^{+}.$ For intermediate upward
velocities, the shape inherits a blend of the two angles.

If the singularity were allowed to proceed all the way to $R_{\mathrm{min}}=0,$
this effect would not matter since all symmetric necks, even squat
ones, progress to a cylinder eventually \citep{PhysRevE.84.026313}.
In reality, however, the asymptotic form progresses so slowly that
other dynamics pre-empt the singularity before the neck can become
slender. Therefore, it is the post-transient shape that sets up the
post-pinch-off effects. For example, we expect that tuning the post-transient
cone angle will play a significant role in the later effects, including
neck air flow, azimuthal vibrations, and the Worthington jet, but
more research is needed.

This work was supported by NSF No.~CBET-0967282 and the Keck initiative
for ultra-fast imaging (University of Chicago). The author thanks
Justin Burton, Michelle Driscoll, Nathan Keim, Lipeng Lai, Sidney
Nagel, Thomas Witten, and Wendy Zhang for discussions and suggestions.