Generation of droplet arrays with rational number spacing patterns driven by a periodic energy landscape

TL;DR

This paper presents an experimental system that uses a periodic magnetic energy landscape to generate droplet arrays with tunable spacing and size, including rational number patterns, and develops a dynamical model to understand these phenomena.

Contribution

The study introduces a novel method for controlling droplet patterns using a time-varying magnetic energy landscape and provides a dynamical systems model for analysis.

Findings

Droplet arrays can be tuned to have rational number spacing patterns.

Droplet sizes can be mono- or multidisperse depending on parameters.

A discrete-time dynamical model reproduces observed patterns.

Abstract

The generation of droplets at low Reynolds numbers is driven by non-linear dynamics that give rise to complex patterns concerning both the droplet-to-droplet spacing and the individual droplet sizes. Here we demonstrate an experimental system in which a time-varying energy landscape provides a periodic magnetic force that generates an array of droplets from an immiscible mixture of ferrofluid and silicone oil. The resulting droplet patterns are periodic, owing to the nature of the magnetic force, yet the droplet spacing and size can vary greatly by tuning a single bias pressure applied on the ferrofluid phase; for a given cycle period of the magnetic force, droplets can be generated either at integer multiples (1, 2, etc.), or at rational fractions (3/2, 5/3, 5/2, etc.) of this period with mono- or multidisperse droplet sizes. We develop a discrete-time dynamical systems model not only to reproduce the phenotypes of the observed patterns but also provide a framework for understanding systems driven by such periodic energy landscapes.

Tables

Table 1. Table 1: Table with nominal dimensions of ‘T’ and ‘I’ bars

Dimensions of ‘T’ and ‘I’ bars
	($\mu m$)
$l_1$	1012.5
$l_2$	1125
$l_3$	1125
$w$	213.5
$t$	25
$g$	70

Equations

$$U(t)=\{\begin{array}{ll}{B}_0{M}_{d}Qtsin(wt+\phi )& ,\text{leavevmode}sin(wt+\phi )\ge 0\\ 0& ,\text{leavevmode}sin(wt+\phi )<0\end{array}$$

$$q{\tau }_{i+1}f({\tau }_{i+1},\text{leavevmode}{\phi }_i)=1$$

Full text

**Generation of droplet arrays with rational number spacing patterns driven by a periodic energy landscape

**

Anatoly Rinberg1, Georgios Katsikis2& Manu Prakash3∗

1School of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA, 2Department of Mechanical Engineering, Stanford University, Stanford, California 94305, USA

3Department of Bioengineering, Stanford University, Stanford, California 94305, USA

∗To whom correspondence should be addressed; E-mail: trinberg@gmail.com

**Keywords

**droplet generation, magnetic energy, phase diagrams, discrete-time dynamical systems, droplet patterns

**Abstract

**The generation of droplets at low Reynolds numbers is driven by non-linear dynamics that give rise to complex patterns concerning both the droplet-to-droplet spacing and the individual droplet sizes. Here we demonstrate an experimental system in which a time-varying energy landscape provides a periodic magnetic force that generates an array of droplets from an immiscible mixture of ferrofluid and silicone oil. The resulting droplet patterns are periodic, owing to the nature of the magnetic force, yet the droplet spacing and size can vary greatly by tuning a single bias pressure applied on the ferrofluid phase; for a given cycle period of the magnetic force, droplets can be generated either at integer multiples (1, 2, etc.), or at rational fractions (3/2, 5/3, 5/2, etc.) of this period with mono- or multidisperse droplet sizes. We develop a discrete-time dynamical systems model not only to reproduce the phenotypes of the observed patterns but also provide a framework for understanding systems driven by such periodic energy landscapes.

Introduction

Discrete-time dynamical systems have been used to study physical phenomena such as population dynamics of predator-prey behavior [M. G. Neubert, M. Kot, and M. A. Lewis(1995), X. Liu and D. Xiao,(2007)], spatial ecological patterns [M. Kot(1992)], control theory [K. Ogata,(1995)] and chaotic electronic circuits [I. Campos-Canton, E. Campos-Canton,(2009)]. These systems often deal with recursive mathematical relations and use iterative maps to describe behaviors such as convergence to stable points, limit cycles and chaos [T. Yoshida, H. Mori, and H. Shigematsu(1983), S. H. Strogatz,(2014), R. M. May(1976)].

In fluidic systems, droplet generation can be thought of as a discrete event, corresponding to the moment when a droplet breaks free from the bulk phase, making droplet generation well positioned to be studied as a discrete-time dynamical system. Yet, little work has explored this connection so far [D. Sessoms, A. Amon, L. Courbin, and P. Panizza,(2010)], despite extensive studies of droplet generation either for technological purposes [J. A. Schwartz, J. V. Vykoukal, and P. R. Gascoyne,(2004), H. Song, D. L. Chen, and R. F. Ismagilov,(2006), S.-Y. Teh, R. Lin, L.-H. Hung, and A. P. Lee,(2008), T. Schneider, J. Kreutz, and D. T. Chiu,(2013)], or fundamental physical understanding. With regards to the latter, there is a conceptual connection between discrete-time dynamical systems and the study of droplet pattern generation, which relates the size of the droplets to their spacings, often revealing asymmetries even at low Reynolds numbers under laminar flow [T. Thorsen, R. W. Roberts, F. H. Arnold, and S. R. Quake,(2001), P. Garstecki, M. J. Fuerstman, H. A. Stone, and G. M. Whitesides,(2006), D. Link, S. L. Anna, D. Weitz, and H. Stone,(2004), P. Garstecki, M. J. Fuerstman, and G. M. Whitesides,(2005)]. Such patterns further enable self-organization phenomena where generated droplets are driven into ordered structures [K. Kita, M. Ichikawa, and Y. Kimura,(2008), J. V. Timonen, M. Latikka, L. Leibler, R. H. Ras, and O. Ikkala,(2013)].

Unlike microchannel configurations for droplet generation [G. F. Christopher and S. L. Anna,(2007)] that induce shearing between the two phases through T-junctions [T. Thorsen, R. W. Roberts, F. H. Arnold, and S. R. Quake,(2001), P. Garstecki, M. J. Fuerstman, H. A. Stone, and G. M. Whitesides,(2006), J. D. Tice, H. Song, A. D. Lyon, and R. F. Ismagilov,(2003)] or flow focusing [S. L. Anna, N. Bontoux, and H. A. Stone,(2003), W. Li, E. W. Young, M. Seo, Z. Nie, P. Garstecki, C. A. Simmons, and E. Kumacheva,(2008)], where the droplet formation timescales arises from balance of viscous forces and capillary pressure, in this work, we report a novel microfluidic system with an intrinsic driving frequency determined by the time-varying magnetic energy landscape with a two-phase immiscible mixture of water-based ferrofluid (FF) and silicone oil. The magnetic energy landscape generates an oscillatory force that produces the droplet arrays whose patterns depend on the energy of breakup, the oscillation frequency and a bias flow-rate. The same concept of magnetic energy landscapes has previously been utilized to synchronously manipulate water-based FF droplets and, through droplet-to-droplet interactions, perform physical logic operations [G. Katsikis, J. S. Cybulski, and M. Prakash,(2015)]. In this letter, we use this platform to demonstrate control over periodic droplet patterns, characterized by different droplet-to-droplet spacing and droplet sizes, and develop a discrete-time dynamical systems model to explain the dynamics driving the formation of these patterns.

Experimental methods

We supply the FF through an inlet tubing (diameter ${d}_{tube}=300\leavevmode\nobreak\ \mu m$) that is placed at a distance $d=50-200\leavevmode\nobreak\ \mu m$ from a substrate covered with a $3-5\leavevmode\nobreak\ mm$ thick film of silicone oil (Fig. 1a, side view; Supp. Information). The FF reservoir is held at a height $h_{ff}$ from the substrate, that creates a differential pressure $\Delta P={\rho}_{ff}\leavevmode\nobreak\ g\leavevmode\nobreak\ h_{ff}$, where ${\rho}_{ff}=1.28\leavevmode\nobreak\ g/{cm}^{3}$ is the density of the FF and $g=9.81\leavevmode\nobreak\ m/{s}^{2}$ is the acceleration of gravity. Due to this pressure, $\Delta P$, there is flow of bulk FF with a rate $Q$.

The droplets are generated through the interaction of the bulk FF with soft-magnetic (permalloy) tracks (characteristic length $\sim 1\leavevmode\nobreak\ mm$) on the substrate via exposure to two magnetic fields. The first magnetic field, $\mid{B}_{z}\mid=250\leavevmode\nobreak\ G$, is perpendicular to the substrate, has a fixed magnitude and polarizes the bulk FF in a uniform manner. The second magnetic field, $\mid{B}_{xy}\mid=40\leavevmode\nobreak\ G$, is in the plane of the substrate, is rotating with a radial frequency $\omega$ and polarizes the tracks. As a result, these magnetic fields create a dynamic, spatiotemporal magnetic energy landscape, where the FF will be driven towards the minimization of its potential energy. To accomplish this, the lower end of the bulk is subject to a magnetic force $\overrightarrow{F}_{mag}$ that extracts sub-millimeter diameter droplets (Fig. 1a, side view). For this study, we restrict ourselves to tracks that have shapes of ‘T’ and ‘I’ bars that ensure that they can be polarized effectively by the $\overrightarrow{B}_{xy}$ and suffice not only to generate droplets but also to propagate them along the tracks (Fig. 1a, top view). For a fixed position of the inlet tube, we show both droplet generation and propagation along the tracks of the substrate (Fig. 1b, top view; Supp. Video 1). To avoid overcrowding the substrate with droplets, we use outlet lines connected to a negative pressure line that remove the droplets from the substrate (Fig. 1a and S1).

Experimental observations

For given magnetic fields and fixed positions of the inlet and outlet tubes, we apply pressures in the range $\Delta{P}=0.5-8\leavevmode\nobreak\ kPa$ and observe that the generated droplet arrays converge to a steady-state pattern within 2-3 cycles of $\overrightarrow{B}_{xy}$. The system is in a constant flow and pressure regime (Supplementary Information). In steady-state, there is a minimum of one full cycle of $\overrightarrow{B}_{xy}$ required to generate a single droplet ($C_{d}=1$). Here, we define $C_{d}$ as the number of cycles of rotating magnetic field per released drop. For decreasing $\Delta{P}$, more cycles are required for the generation of single droplet ($C_{d}\geq 1$), resulting in droplet arrays that are less tightly spaced (Fig. 2 and Supp. Video 2). Interestingly, we observe that the spacings between consecutive droplets can be non-constant, nonetheless still having a periodic pattern; for example there can be periodic alternation between one and two cycles per droplet (i.e. ${C}_{d}$ exhibits sequence ‘..1-2-1-2..’), resulting in an average of $\bar{C_{d}}=3/2$ (Fig. 2b). In these cases, the volumes of the droplets can also be different.

Model

To explain the different droplet-to-droplet spacing and individual droplet volumes in our generated arrays (Fig. 2b), we develop a theoretical model. We write a tractable expression for the magnetostatic energy of the droplet, which theoretically is defined as $U=\int{-\overrightarrow{M}\cdot\overrightarrow{B}_{bar}dV}$, where $\overrightarrow{M}$ is the magnetization of the droplet, $\overrightarrow{B}_{bar}$ is the magnetic field generated by the bars and $V$ is the volume of the droplet. To simplify the complicated expression for $U$ (Supplementary Information), we base our model on the following five assumptions: First, we consider the droplet as a point mass and write $U=-\overrightarrow{M}\cdot\overrightarrow{B}_{bar}V$. Second, we assume that $\overrightarrow{M}=M\hat{z}$ (Fig. 1a) with $V$ increasing linearly over time $t$ for a given flow rate $Q$, allowing us to write the magnitude of the magnetic moment $\overrightarrow{\mu}=\overrightarrow{M}V$ as $\mu(t)=M_{d}Qt$. Third, we assume that $\overrightarrow{B}_{bar}={B}_{bar}\hat{z}$ with ${B}_{bar}$ varying as a sine wave over time, consistent with the oscillatory nature of $\overrightarrow{B}_{xy}$, and thus write ${B}_{bar}(t)=B_{0}sin(\omega t+\varphi_{i})$ where $B_{0}$ is positive and is the maximum amplitude of $\overrightarrow{B}_{bar}$, $\omega$ is the angular frequency, and $\varphi_{i}$ is the phase of $\overrightarrow{B}_{xy}$. Fourth, we assume that a droplet breaks up from the bulk when its energy $U$ is minimized to a threshold $U_{breakup}$ which is constant and does not depend on droplet volume. We base this assumption on the fact that, for droplets that are roughly the diameter of the inlet tube or larger, $U_{breakup}$ is set by the cross-sectional area of the inlet tube and the surface energies of the fluids. In our experiments, the radius of the smallest droplet was measured to be $r_{min}=240\mu m$, suggesting a constant $U_{breakup}$. Additionally, for the rest of this work, we will refer to the absolute value of the energy $U$. Fifth, we assume that droplet breakup can occur only in the attractive phase of the oscillation when $sin(wt+\varphi)>0$ and ${B}_{bar}(t)>0$. In the repulsive phase, the droplet is pushed away from the magnetized bar, which then reduces the applied magnetic force on the droplet, preventing breakup from occurring.

Combining all five of these assumptions, we write down the equation for the magnetostatic energy of the model as:

[TABLE]

Once a droplet is released, only the phase of $\overrightarrow{B}_{xy}$ at the previous breakup is needed to determine the time to next breakup. This allows us to write equation (1) as a recursive formula; assuming that a droplet ‘i’ is generated at time $t=\tau_{i}$ when $\overrightarrow{B}_{xy}$ is at angle ${\varphi}_{i}$, then the next droplet ‘i+1’ will be generated at time $t=\tau_{i+1}$ and ${\phi}_{i+1}$, which occurs when the droplet magnetic energy is equal to $U({\tau}_{i+1})={U}_{i+1}={U}_{breakup}$ (Fig. 3b). Without loss of generality, we reduce equation (1) by setting $B_{0}M_{d}Q=q$ (${s}^{-1}$), $\omega=2\pi$ ($rad/sec$) and ${U}_{breakup}=1$, and write the recursive expression as:

[TABLE]

Where $f$ is the waveform of the magnetic field relevant for breakup, and is given by $f(t,\leavevmode\nobreak\ \varphi_{i})=\textrm{max}(sin(2\pi t+\varphi_{i}),\leavevmode\nobreak\ 0)$ (Fig. 3a, gray field magnitude curve). Next, we solve equation (2) to reproduce the phenotype of the droplet arrays generated experimentally (Fig. 2b). For given $q$ and angles $\varphi_{i}$ in the range $[0,\pi]$, we find the corresponding values of $\tau_{i+1}$. We restrict our parameter range for $\varphi_{i}$ to an upper bound of $\pi$ since no breakup can occur from $\pi$ to $2\pi$. Then, we calculate both the angle $\varphi_{i+1}$ based on equation $\varphi_{i+1}=\textrm{mod}_{2\pi}(2\pi\tau_{i+1}+\varphi_{i})$ and the number of cycles $C_{d}$ required to generate a droplet ‘i+1’ based on $C_{d}=\textrm{quotient}_{2\pi}(2\pi\tau_{i+1}+\varphi_{i})$, therefore generating phase maps for specific $q$ values that relate $\varphi_{i}$ to $\varphi_{i+1}$ (Fig. 3b,c). For $q=1.1$, $\varphi_{i}$ converges to a single steady-state angle $\varphi_{s}$ (Fig. 3b) resulting in monodisperse droplets with $\bar{C}_{d}=1$ (Fig. 3c) independent of the initial $\varphi_{0}$. In other regimes, for example at $q=0.9$ and $q=0.864$, $\varphi_{i}$ periodically alternates respectively between two and three steady-state angles (Fig. 3c) resulting in multidisperse droplets with $\bar{C}_{d}=3/2$ and $\bar{C}_{d}=5/3$ in qualitative agreement with experiments (Fig. 2b).

To study the stability and pattern-space of the model, we conduct a parameter sweep of $q$ in the range: $[0.15,1.5]$ (Fig. 4). The phase-stability map reveals domains of single steady-state points, where $\varphi_{i}=\varphi_{i+1}$ with integer $\bar{C_{d}}$ values (Fig. 4a; black lines), that are interrupted by domains of multiple steady-state points and non-integer $\bar{C_{d}}$ values (Fig. 4a; red lines). These multidisperse transitional domains occur at discontinuous boundaries in the phase map (qualitatively as in Fig. 3c, $\bar{C_{d}}=3/2,\leavevmode\nobreak\ 5/3$). In addition, for the explored parameter range, we find that given any initial $\varphi_{i}$ value, the subsequent $\varphi_{i+1}$ is always narrowed to a band of $[0.509,\leavevmode\nobreak\ 1.771]\leavevmode\nobreak\ rad$ (Fig. 4a, colorbar; Supp. Info).

Furthermore, to illustrate the richness in potential droplet spacing and volume patterns, we calculate droplet volume over discrete cycle intervals at different $q$ values (Fig. 4b). The pattern-space includes monodisperse and multidisperse droplet sequences at $\bar{C_{d}}$ (Fig. 4c) values observed experimentally (Fig. 2c).

Comparison of experiment and model

To understand the relationship between the droplet volume and pressure, we study one configuration at an in-plane frequency of $2\leavevmode\nobreak\ Hz$, describe the measured physical quantities in detail and test our analytic model by comparing to the experimental results (Fig. 5).

Decreasing pressure down from $8\leavevmode\nobreak\ kPa$, we find monotonically increasing $C_{d}$ values (Fig. 5a). For a given $\bar{C_{d}}$ value, average droplet volume decreases with decreasing pressure. As $\bar{C_{d}}$ transitions from 1 to 2, 2 to 3/2 and 3/2 to 3, droplet volumes jump abruptly to higher values before decreasing again (Fig. 5b). We find that the average minimal droplet volume for all integer $\bar{C_{d}}$ is $V_{min}=0.059\mu l$ (Fig. 5b, dashed line). Plotting the flow rate, $Q=V_{drop}f/\bar{C}_{d}$, as function of $\Delta P$ gives a linear relationship with a slope of $26.2*10^{-3}\,\mu l\,s^{-1}\,kPa^{-1}$ (Fig. 5c, $R^{2}=0.997$). The linearity of this relationship confirms that the magnetic, capillary and hydrodynamic forces at the exit of the tube are much smaller than the force driving the FF flow.

Given the experimentally determined $V_{min}$, we can reevaluate equation (2) and compare theory to experiment, by parameterizing $q=Q/V_{min}$ and setting $w=2*2\pi$. We use the recursive equation (2) to numerically solve exact values of $\bar{C_{d}}$ (Fig. 5a, gray line) and the droplet volumes (Fig. 5b) for $V_{min}=0.059\leavevmode\nobreak\ \mu l$ over a range of $q$. With $V_{min}$ as the single-parameter fit, we find good qualitative agreement between experiment and theory, particularly in the transitions between different $\bar{C_{d}}$. For $\bar{C_{d}}=3/2$, we find a difference in expected droplet volumes suggesting that there may need to be important corrections made to the $B_{xy}$ waveform.

Conclusions

In summary, we have demonstrated an experimental platform in which a periodic force generates droplet arrays with complex patterns of droplet spacings and sizes. We have developed a discrete-time dynamical systems model to explain the observed patterns, and found good agreement with experimental measurements. More broadly, this work may suggest a new formalism to study droplet generation under time-dependent force using iterative phase maps and other discrete-time dynamical systems approaches.

Acknowledgements

The authors thank all members of the Prakash Lab for useful discussions. AR acknowledges support from the NSF GRFP. GK acknowledges support from the Onassis Foundation and A. G. Leventis Foundation. MP acknowledges support from Pew Foundation and Keck Foundation.

1   Experimental methods

Fabrication of fluidic chips $-$ The ‘T’ and ‘I’ bars are fabricated by etching permalloy foils that are epoxy-bonded on glass substrates, using an protocol identical to reference [Katsikis2015]. The ‘T’ and ‘I’ bars (Fig.S1) have millimeter-size dimensions (Table 1). The permalloy bars are coated with teflon and the fluidic chips do not have a top cover.

Magnetic fields $-$ The magnetic fields are generated using the system of electromagnetic coils described in reference [Katsikis2015]. The ratio between the magnitudes of the magnetic fields is $\mid{B}_{z}\mid/\mid{B}_{xy}\mid\geq 5$, thus ensuring that the induced magnetization of the generated droplets is along the z-axis (Fig. 1A). However, the induced magnetization of the metallic bars is always in the x-y plane as they are too thin (for example ${t}/{l}_{1}=1/40$) to support magnetization in the z-axis.

Two-phase mixture of fluids $-$ The mixture consists of two phases. The first phase is silicone oil (Sigma Aldrich, CAS: 63148-62-9, kinematic viscosity $5\leavevmode\nobreak\ cSt$, density $0.913\leavevmode\nobreak\ kg/{m}^{3}$) which is pipeted on the surface of the fluidic chip forming a film of thickness $h_{oil}=3-5mm$ beneath the open air-oil interface. The second phase is water-based ferrofluid (Ferrotec EMG 700, kinematic viscosity $5\leavevmode\nobreak\ cSt$, density $1.28\leavevmode\nobreak\ kg/{m}^{3}$) which is dispensed on the film using an inlet tubing (Fig. 1A).

Inlet tubing $-$ The inlet tubing is made of teflon (PTFE) with internal diameter $300\leavevmode\nobreak\ \mu m$ and length $1\leavevmode\nobreak\ m$. The first tip of the tubing is suspended at a height $d=50-200\leavevmode\nobreak\ \mu m$ above the permalloy bars. This height $d$ is always smaller than the thickness of the silicone oil film on the substrate, that is $d<{h}_{oil}$, thus making this tip completely immersed in the film. The second tip of the tubing is connected to a ferrofluid reservoir whose top surface is at a height $h_{ff}=10-80\leavevmode\nobreak\ mm$ above the permalloy bars. This height $h_{ff}$ creates a pressure difference $\Delta P$ generating flow that fills the tubing with ferrofluid and - via the first tip - dispenses it into the substrate. The pressure difference $\Delta P$ is adjusted by adjusting the height $h_{ff}$ of the ferrofluid reservoir. Furthermore, the inlet tubing is threaded through a glass capillary with internal diameter $500\mu m$, which is mounted on a three-axis translational stage for adjusting the position of the end of the inlet tubing relative to the permalloy bars. The height of the oil ($h_{oil}$) contributes an insignificant reduction in pressure and is not considered here.

Outlet tubing $-$ The outlet tubing is made of teflon, similar to the inlet tubing (Fig. 6a). At its lower end that is in proximity to the substrate, it also contains a blunt-tip pin made of stainless steel (23 gauge). The magnetic field $\overrightarrow{B}_{z}$ along the z-axis magnetizes the pin. The magnetized pin attracts the ferrofluid droplets and by also using an additional negative pressure difference across the outlet tubing, the droplets that reach the outlet tubing are removed from the substrate (Fig. 1 b and Fig. 6b).

Imaging $-$ Droplet volume measurements are performed by imaging the chip with a dSLR (Canon T3i, Canon EF 100mm f/2.8L IS USM Macro Lens).

PTFE-Oil-Ferrofluid Surface Energy. In order to estimate the volume of sessile droplets, by only imaging from the top, we measured the contact angle between ferrofluid, PTFE in silicone oil. We measured 11 droplets from the side, sessile on a PTFE surface, for an average surface angle of $\theta=24.86\pm 2.72$ (Fig. 7).

2   Data Analysis

Droplet Volume Measurement. For each measurement, droplets are first generated and then all B-fields are turned off, so that the droplets are in a sessile state on the chip surface. Droplet radii are measured using contrast-based object detection in Matlab ( YORGOS - add any comments here). For maximal droplet volumes of $V_{droplet}\approx 0.12ul$, the Bond number is $\approx 0.25$ ($\Delta\rho=\rho_{ff}-\rho_{oil}=0.2g/ml$; $\gamma=3mN/m$ [Flament1996]), therefore justifying the spherical cap assumption in calculating the volumes of the droplets, where the $V_{cap}(r,\theta)=(\pi r^{3}/6)(1-cos{\theta})(3sin{\theta}^{2}+(1-cos{\theta})^{2})$. In our system, the FF-teflon-oil surface contact angle is measured to be $\theta\approx 25\degree$.

3   Droplet Generation Videos

Video 1. Monodisperse droplet generation at $B_{xy}$ frequency of $2Hz$.

Video 2. Pressure sweep of droplet generation at $B_{xy}$ frequency of $2Hz$, displaying various $\bar{C_{d}}$ regimes.

3.1   Model and Fits

Computational Solution. MATLAB R2014a was used to numerically solve the recursive equation (2). The recursive process is as follows: after the i-th droplet is generated, time is reset to $t=0$ and $\varphi_{i}$ is propagated to the subsequent iteration. We next solve for the time, $\tau_{i+1}$, that it takes for the energy to reach $U_{breakup}$. To plot phase maps, we solve the recursive equation for a range of $\varphi$ from [math] to $\pi$ in increments of at least $0.001$.

Fitting. Linear fitting was done using the first-order Polyfit function in Matlab. $R^{2}$ value was then calculated as an estimator of linearity.