# Temperature Dependence of Inertial Pumping in Microchannels

Pavel E. Kornilovitch<EMAIL_ADDRESS>HP Inc., Corvallis, Oregon
97330 USA  Tyler Cochell  HP Inc., Corvallis, Oregon 97330 USA  Alexander N.
Govyadinov  HP Inc., Corvallis, Oregon 97330 USA

###### Abstract

Inertial pumping is a promising new method of moving fluids through
microchannels but many of its properties remain unexplored. In this work,
inertial flow rates are investigated for different channel lengths, operating
temperatures, and resistor pulse energies. The flow in closed channels is
visualized by adding fluorescent tracer beads to the test fluid (pure water).
A robust methodology of extracting flow rates from high-resolution video
recordings is developed. Flow rates are found to scale inversely with the
channel length. The observed dependence is explained based on a simple
phenomenological “kick” model of inertial pumping. Flow rates are also fitted
to the more fundamental one-dimensional model of inertial pumping from which
the intrinsic drive bubble strength is extracted. The measured flow rates vary
strongly with temperature. For well-developed drive bubbles, flow rates at
$T=70\,^{\circ}\\!{\rm C}$ are about 12$\times$ higher than at
$T=30\,^{\circ}\\!{\rm C}$. Three separate effects contribute to increasing
flow rates at high temperatures: (i) lower viscosity of the test fluid, (ii) a
stronger drive bubble, and (iii) increasing mechanical efficiency of the pump,
i.e., better conversion of the drive bubble strength to unidirectional post-
collapse kick. Relative contributions of the three effects are quantified. The
energy dependence of flow rates exhibits a clear saturation behavior. The
bubble strength is fitted to a phenomenological saturation model. In the end,
a complete predictive length-temperature-energy model of flow rates is
constructed. The observed strong temperature dependence of inertial pumping
should be considered when designing microfluidic workflows. It also highlights
the need for integrated flowmeters that could stabilize complex flow patterns
via sensory feedback.

##  I  Introduction

Inertial micropumps are an emerging microfluidic technology capable of
effectively pumping fluids in 10-500 $\mu$m-wide channels without need for any
external pressure source. Yuan1999 ; Yin2005a ; Yin2005b ; Torniainen2012 ;
Govyadinov2016  Based on the commercial thermal inkjet (TIJ) technology,
Stasiak2012  inertial micropumps can be fabricated in high densities up to
tens of thousands per square centimeter. They can drive extensive fluidic
networks performing complex workflows using only electrical power and can
eventually become a key enabler of the true “lab-on-a-chip,” as opposed to
“chip-in-a-lab,” devices. Physically, an inertial pump is a resistive
microheater (thermal inkjet, or TIJ, resistor) located inside a closed
microfluidic channel. A microsecond-long current pulse vaporizes a thin
interfacial layer of fluid and creates a vapor drive bubble, which performs
mechanical work through expansion. If the resistor is placed asymmetrically
between two reservoirs, the expansion-collapse dynamics of the drive bubble
results in a unidirectional flow from the shorter to the longer side of the
channel. Physics behind pump operation is complex. Yuan1999 ; Torniainen2012 ;
Govyadinov2016 ; Kornilovitch2013  A key role is played by vortices formed in
reservoirs, Torniainen2012  which dissipate fluid’s momentum and break the
expansion-collapse symmetry. Momentum dissipation proceeds faster on the
shorter side of the channel. As a result, the shorter side refills faster and
arrives at collapse with a higher mechanical momentum, which translates into
residual post-collapse flow. Transition from expansion to refill is delayed on
the longer side due to a larger mass of fluid, hence, the term “inertial
pumping.”

Inertial pumping was first theorized Yuan1999  in 1999 and later demonstrated
Yin2005a ; Yin2005b  by Prosperetti. In 2009, the Hewlett-Packard (HP) group
observed Torniainen2012  inertial pumping in microfluidic devices fabricated
on the same manufacturing platform as HP’s commercial printheads. A
phenomenological reduced-parameter one-dimensional model of inertial pumping
was developed in Ref. Kornilovitch2013, . It was based on the original work by
Yuan and Prosperetti. Yuan1999  The model was recently extended to two pumps
acting simultaneously by Hayes et al. Hayes2021  Careful analysis of single-
pulse dynamics using a million-frame-per-second imaging technique
Govyadinov2016  confirmed basic understanding of the pumping mechanism and
validated both phenomenological and first-principle CFD treatments of the
effect. Multiple pumps working within a branched fluidic network were
demonstrated in Ref. Hayes2018, .

Figure 1: (a) Fluidic architecture. Channel lengths are given in Table 1. The
red rectangle marks the inertial pump of channel 8. The “reservoir” is a
through-silicon slot that delivers test fluid from the back of the die. The
reservoir acts simultaneously as an inlet and outlet when the fluid circulates
in the anticlockwise direction. The scalebar size is 100 $\mu$m. (b) Optical
diagram. (c) Photograph of the optical system.

Being inertial in nature, the pumping is expected to work better at high
Reynolds numbers, at least up to centimeter-wide channels. At those
dimensions, instantaneous boiling is inefficient and thermal actuation needs
to be replaced with mechanical impacts or other means. Inertial pumping was
also demonstrated with laser Wang2004 ; Dijkink2008  actuation. However, it is
the low-Re limit that makes inertial pumping with thermal actuation stand out.
As channels become smaller, the resistor width decreases proportionally,
keeping the bubble power per die surface area constant. (Peak pressures inside
drive bubbles reach tens of bars.) However, since the channel height also
decreases, bubble power per cross-sectional area of the channel goes up. This
positive scaling partially offsets a fourth-power increase in channel’s
fluidic resistance and keeps inertial pumping effective down to channel widths
of 10 $\mu$m and less. Magnified by pulse frequencies of tens of kilohertz,
even a sub-picoliter-per-pulse pump can deliver flow rates sufficient for many
applications.

Many physical properties of inertial pumping remain uncharacterized. In this
work, we systematically investigate variation of flow rates with channel
length $L$, operating temperature $T$, and resistor pulse energy $E$. We
develop a robust method of extracting flow rates from raw videos based on
identifying the fastest tracer particle in Sec. II.3. To interpret results, we
derive a simple phenomenological kick model of inertial pumping, which is
appropriate for long channels in Sec. IV.2. Experimentally, we find flow rates
to scale $\propto 1/L$, consistent with a $\propto L$ increase in fluidic
resistance. Temperature variation of flow rates is found to be significant.
Between $T=30\,^{\circ}\\!{\rm C}$ and $T=70\,^{\circ}\\!{\rm C}$ flow rates
of pure water increase by as much as 12$\times$. The flow rates rise sharply
because of (i) falling viscosity, (ii) stronger drive bubbles, and (iii)
improving mechanical efficiency of the pump. By fitting experimental data to
phenomenological models and extracting intrinsic bubble strengths, we are able
to quantify the relative contributions of all three mechanisms. This is done
in Sec. IV. The energy dependence of flow rates shows a clear saturation
behavior indicating transition to fully developed drive bubbles. In the end, a
complete predictive Length-Temperature-Energy model of inertial pumping is
constructed. Our conclusions are presented in Sec. V.

##  II  System and Methods

###  II.1  Fluidic architecture

Microfluidic devices were fabricated using the same manufacturing processes as
HP’s commercial inkjet printheads. Stasiak2012  The starting substrates were
200 mm silicon wafers with CMOS driving electronics. Microheaters (TIJ
resistors) and metal contacts were fabricated by thin-film deposition and
patterning techniques. The fluidic layer was then formed on top of resistors
using a multilayer SU8 deposition and patterning process. Finally, slots were
etched through the silicon substrate to enable fluid delivery from the back of
the die.

The final architecture is shown in Fig. 1(a). It includes series of eight
U-shape channels connected to one common fluid reservoir. The channel width is
30 $\mu$m and channel height is 31 $\mu$m. Channel lengths are listed in Table
1. Each channel is equipped with an inertial pump – a $20\times 30$ $\mu$m2
resistive heater – in the east leg of the channel. When the resistor starts
firing, fluid begins to circulate through the channel in an anticlockwise
direction.

Channel | Length, $\mu$m | Channel | Length, $\mu$m
---|---|---|---
1 | 396 | 5 | 1196
2 | 466 | 6 | 1266
3 | 796 | 7 | 1596
4 | 866 | 8 | 1666
Table 1:  Total lengths of the U-shape channels studied.

###  II.2  Experiment

The test fluid was pure de-ionized water. To visualize flow, 2.0 $\mu$m
diameter fluorescent tracer beads (Polybead® Microspheres, Polysciences, Inc.)
were added. Bead density was adjusted to have between one and ten beads in a
channel at any given time, so that the tracers could be considered Lagrangian.
The top SU8 layer (channel “ceiling”) was optically transparent which allowed
for easy observation of the tracers. The pumps were driven by electric firing
pulses 1.0 $\mu$s in duration, with a total heat energy released in a single
pulse varied between 0.896 and 1.186 $\mu$J. The firing frequency was varied
between 10 and 50 Hz. Each firing created a vapor drive bubble that resulted
in a unidirectional flow pulse via the inertial pumping mechanism.
Torniainen2012 ; Kornilovitch2013  Individual flow pulses lasted less than 400
$\mu$s for all the conditions studied in this work. Thus, short pulses were
separated by long (20-100 ms in duration) waits in between. Since the pulses
were independent, the overall flow rate was simply the net-flow-per-pulse
times the firing frequency. Conversely, net-flow-per-pulse, which is the basic
intrinsic property of the pump, could be obtained by dividing average flow
rate by the firing frequency.

Figure 2: Particle tracking sequence. (a) Dry microchannels. (b) Channels
filled with water. The bright white spots are tracer particles added to water
to visualize flow. (c) Particle tracking software Fiji TrackMate Tinevez2017
identifies tracer spots within a region of interest (yellow rectangle) and
marks them with magenta circles. (d) By linking spots from different video
frames, the software constructs particle tracks. In this image, the tracks are
colored according to their mean velocity. Faster tracers are grouped near the
center of the channel, which is expected in developed laminar flow.

Our microfluidic chips included a built-in temperature control mechanism
comprising FET heaters, temperature-sensitive resistors (TSRs), and CMOS
circuitry that implemented a digital feedback loop. Die temperature could be
specified in the control software. Although the TSRs returned temperature near
the bottom of fluidic channels, computer modeling indicated that temperature
variation throughout the channel thickness (31 $\mu$m) was less than
$1\,^{\circ}\\!{\rm C}$ and could be neglected for our purposes.
Experimentally, die temperature was varied between $T=30\,^{\circ}\\!{\rm C}$
and $T=70\,^{\circ}\\!{\rm C}$, but some results were also collected at
$T=80\,^{\circ}\\!{\rm C}$. Characterizing temperature dependence of inertial
pumping was the main goal of this investigation.

The optical system utilized to observe moving tracers was similar to the one
described in Ref. Hayes2018, . The optical diagram is shown in Fig. 1(b) and
real setup in Fig. 1(c). Observations were performed using a custom-designed
inverted optical microscope with 10$\times$ magnification. The basic design
described earlier Torniainen2012  was modified to enable fluorescent imaging
capability by adding a set of filters and dichroic beam splitter using
illumination from a Lightspeed Technologies HPLS-36DD18B series single emitter
with a built-in driver for strobing a 478 nm blue LED. The frame rate of the
recording camera was synchronized with the pumps’ pulse rate, typically at a
1:1 ratio. Thus, tracers’ displacement observed over, say, 200 video frames
was the result of 200 consecutive pump pulses. Video records were processed
using Fiji software package Schindelin2012  and Fiji TrackMate plugin.
Tinevez2017  The result was a series of tracer tracks and their respective
mean velocities in pixels-per-frame as shown in Fig. 2. The region of interest
was limited to straight channel segments, which ensured that tracers moved
primarily along the $y$ axis and their $x$ coordinates were approximately
constant. Brownian motion of 2.0 $\mu$m tracers was unobservably small in our
experimental setup. The mean velocity of each track was, therefore,
interpreted as an estimate for the flow field $v(x)$ at the track’s $x$
position in the channel. Both $x$ and $v$ were converted from pixels to
micrometers using a known scale factor. The number of tracks extracted from
raw video records varied between 100 and 1000 for different experimental
conditions. Such statistics enabled reconstruction of the entire flow profile
and subsequent calculation of average flow rates with sufficient accuracy as
explained in the next section.

Figure 3: Flow rate measurement workflow illustrated on the loop 4,
$T=70\,^{\circ}\\!{\rm C}$, $E=1.053$ $\mu$J data set. (a) The starting point:
a scatterplot of tracer velocities vs position in the channel produced by Fiji
TrackMate, Tinevez2017  after removal of outliers. (b) Half of tracks are
randomly selected to form a subsample. (c) The sabsample tracks are grouped
into 15 bins of equal width. The number of bins $K$ should be chosen
carefully. If $K$ is too large, there will be too few tracks in each bin and
the subsample flow profile will fluctuate wildly. However, a small $K$ does
not resolve the flow profile enough, which reduces quality of the subsequent
fit. (d) Within $k$-th bin, a maximum velocity $v_{p}(x_{k})$ is computed,
which is then taken as an estimate of the flow profile at the bin’s midpoint
$x_{k}$. Collectively, $v_{p}(x_{k})$ forms a subsample flow profile. Then,
the set of pairs $\\{x_{k};v_{p}(x_{k})\\}$ is fitted to the theoretical
profile of Eqs. (28) and (29) using the $x$ coordinates of the two zeros and
the overall height as three adjustable parameters. The height of the optimum
pseudo-parabola is the $v_{\rm max}$ of the given subsample. Note that during
optimization, we used an asymmetric price function: points below the
theoretical curve were priced five times less than points above the curve.
This is because the points below the curve may result from insufficient data
(a bin may be missing a tracer with velocity close enough to the true profile,
which artificially depressed the sample profile at the corresponding $x_{k}$)
while the points above the curve may not. (e) Steps (b)-(d) are repeated 200
times producing a distribution of $v_{\rm max}$. The mean value of the
distribution is taken as the final estimate of the true maximal velocity, and
half-width as one standard deviation. (f) The original scatterplot of (a)
superimposed on the best-fit theoretical profile as well as profile’s upper
and lower bounds.

###  II.3  Flow rate measurements

Once a set of tracks is identified, one way to determine an overall flow rate
is to average tracks’ velocities. This method implicitly assumes that tracers
are distributed uniformly across the channel cross section. However, because
of several physical effects, this assumption is never valid. First, due to
their finite size, the tracers cannot get too close to the channel’s walls, so
the low-velocity parts of the flow field are undersampled. Second, the tracers
can settle or cluster. They also interact with the walls. Finally, the tracers
are disturbed by the drive bubble when flow over the pump, which moves them
across the channel in an unpredictable manner. Figure 3 shows that even a
high-statistics set of tracks clearly lacks sufficient uniformity to be sure
that simple averaging will produce an accurate result. Additionally, it is not
easy to assess the cross-sectional distribution experimentally, as the
tracers’ $z$ coordinates remain unresolved in our setup.

We found that a more robust method of calculating the overall flow rate relies
on identifying the fastest tracer. In a well-developed Poiseuille flow, there
exists a one-to-one correspondence between the maximal and average velocities.
For the simplest cross-sectional shapes, the ratio of the two is known
theoretically. For rectangular channels, the ratio is given in Appendix, Eqs.
(24)-(27). For 30/31 channels studied in this work, $v_{\rm max}/\langle
v\rangle=2.096$, see Fig. 11(a).

There is a legitimate question about applicability of fully developed flow
profiles to pulsatory flows generated by inertial pumps. All our prior studies
based on published Torniainen2012 ; Govyadinov2016  and unpublished CFD
simulations have shown that net fluid displacement after an inertial pulse
matches a fully developed flow profile with deviations less than 1% at any
point of the cross-section. Although the topic deserves further investigation,
the two flow types can be considered indistinguishable for practical purposes,
as statistical variations result in standard deviations of order 2%-3% as
shown below.

Figure 4: Inertial flow rates in pL/pulse vs channel length $L$. Each panel
corresponds to a particular resistor energy and multiple curves within a panel
correspond to different operating temperatures. The symbols are mean
experimental values, and error bars are three times the standard deviation
computed by the resampling procedure. Dividing flow rates by the channel
cross-sectional area of 930 $\mu$m2 yields an average linear displacement per
pulse. For example, a flow rate of 40 pL/pulse corresponds to an average
linear displacement of 43 $\mu$m/pulse. The dashed lines are best fits to the
1d model of inertial pumping, see Sec. IV.1. Phenomenological bubble strengths
extracted from the fits are given in the legends in units of kg m/s. The panel
with a dashed outline ($E=0.969$ $\mu$J) also includes data for
$T=75\,^{\circ}\\!{\rm C}$ and $80\,^{\circ}\\!{\rm C}$.

The calculation workflow is illustrated in Fig. 3. Once a set of tracks
$\\{v_{n}\\}$ is identified and cleaned from outliers, Fig. 3(a), taking the
max operation is all that is needed in the case of large and dense data sets.
Then, the flow rate is simply

$Q{\rm[pL/pulse]}=\frac{{\rm max}\\{v_{n}\\}[\mu{\rm m/pulse]}\cdot A[\mu{\rm
m}^{2}]}{2.096\cdot 1000[\mu{\rm m}^{3}/{\rm pL}]}\>,$ (1)

where $A=30\times 31=930$ $\mu$m2 is the cross-sectional area and the
coefficient 1000 in the denominator converts cubic micrometers to picoliters
that are more natural volume units for our system. This simple procedure is
justified by an expectation that within a large set of tracks, there is always
one that is sufficiently close to the channel center, and the sample maximum
is a good estimator of a true maximum velocity.

The described process does not work very well at low statistics when the
presence of a tracer near the channel center is not guaranteed. The max
operation might also be skewed by unobvious outliers near $v_{\rm max}$. The
method’s robustness can be improved further by considering the entire flow
profile $v_{\rm max}(x)$, that is the maximal flow velocity at a given point
$x$ across the channel. The function $v_{\rm max}(x)$ is also known
theoretically for any channel aspect ratio, see Appendix, Eqs. (28) and (29).
Anchoring the side arms of the profile enables robust determination of the
profile’s pinnacle even in the absence of well-defined tracks near the maximal
velocity.

Figure 5: Inertial flow rates in pL/pulse vs. operating temperature $T$. Each
panel corresponds to a particular resistor energy and multiple curves within a
panel correspond to different channel lengths. The symbols are mean
experimental values, and error bars are three times the standard deviation.
The panel with a dashed outline ($E=0.969$ $\mu$J) also includes data for
$T=75\,^{\circ}\\!{\rm C}$ and $80\,^{\circ}\\!{\rm C}$. The lines are guides
to the eye only.

To obtain a sample flow profile, the set of tracks was split into 12-15 bins
according to the tracks’ $x$ position in the channel, Fig. 3(c). Within each
bin, a maximal sample velocity $v_{p}(x_{k})$ was taken as an estimator of the
true velocity profile at an $x_{k}$ corresponding to the bin’s midpoint. Then,
the sample flow profile $\\{v_{p}(x_{k})\\}$ was fitted to the theoretical
shape of Eqs. (28) and (29) using the $x$ coordinates of the two zeros and the
height as three adjustable parameters, Fig. 3(d). A standard Matlab function
could accomplish a three-parameter optimization in about 1 s. Once the $v_{\rm
max}$ was determined from the fit, the overall flow rate was computed from Eq.
(1).

Finally, standard deviations of flow rates were computed by resampling the
tracks at a 50% sub-sample ratio, Fig. 3(b). For example, from a sample of 300
tracks, 150 tracks were selected at random, and a sub-sample flow rate was
computed using the procedure described above. The process was then repeated
200 times resulting in a distribution of flow rates. The mean value of the
distribution was taken as the final flow rate estimate for the given
experimental conditions. Half-width of the distribution was taken as flow
rate’s standard deviation, Fig. 3(e). The error bars shown in subsequent
figures are three-times the standard deviation obtained with the resampling
method.

Figure 6: Inertial flow rates in pL/pulse vs. resistor energy $E$. Each panel
corresponds to a particular channel length and multiple curves within a panel
correspond to different temperatures. The symbols are mean experimental values
and error bars are three times the standard deviation. The lines are guides to
the eye only.

##  III  Results

In this study, we measured inertial flow rates for different channel lengths,
operating temperatures, and resistor energies. There are multiple ways of
presenting flow rates vs three parameters. First, we show flow rate variation
with channel length $L$ in Fig. 4. Each panel corresponds to a fixed resistor
energy and each line within a panel corresponds to a fixed temperature. The
dashed lines are best fits to the 1D model of inertial pumping.
Kornilovitch2013  In our experimental conditions, net pumping effect is
dominated by a post-collapse flow, which is inversely proportional to total
fluidic resistance, as explained in Sec. IV. Since fluidic resistance grows
linearly with $L$, it explains the observed $1/L$ variation of flow rates.
Sharp flow rate increases with temperature are also apparent in the data.

Figure 7: Top lines: intrinsic bubble strength $q_{0}$ extracted from
experimental data by fitting the $L$ dependent flow rates of Fig. 4 to the 1D
model of inertial pumping. Kornilovitch2013  Bottom lines: kick strength $q$
extracted from experimental data by fitting the $L$ dependent flow rates of
Fig. 4 to the kick formula, Eq. (5).

Next, we show the $T$-dependence of flow rates in Fig. 5. Steep growth is
observed for all resistor energies and all channel lengths. For example, at
$E=1.186$ $\mu$J, the ratio of flow rates at $T=70\,^{\circ}\\!{\rm C}$ and
$T=30\,^{\circ}\\!{\rm C}$ is about a factor of 12. For one particular energy,
$E=0.969$ $\mu$J, we investigated even higher temperatures and found a further
increase up to at least $80\,^{\circ}\\!{\rm C}$. Beyond $80\,^{\circ}\\!{\rm
C}$ inertial pumping becomes less stable due to overheating and formation of
air bubbles mixed with water vapor which complicates measurement and analysis.
There are three physical effects that contribute to high flow rates at
elevated temperatures. First, it is the decreasing viscosity of the test
fluid. The dynamic viscosity of water at $70\,^{\circ}\\!{\rm C}$ is one half
of that at $30\,^{\circ}\\!{\rm C}$: $\eta=0.40$ vs $0.80$ mPa s. AntonPaar  A
smaller viscosity decreases fluidic resistance and increases the flow rate.
Thus, a factor of 2 out of 12 can be assigned to a changing viscosity. The
remaining factor of 6 is attributed to the pump getting stronger with
temperature. The latter effect can be further subdivided into two factors.
First, the drive vapor bubble gets stronger. As the fluid is preheated, a
progressively smaller fraction of energy released in a resistor pulse must be
spent to heat the fluid to the boiling temperature which leaves more energy to
produce vapor. Second, conversion of the initial bubble expansion to
unidirectional push, which is pump mechanical efficiency, also improves with
temperature. The analysis presented in Sec. IV enables to rigorously define
pump efficiency and separate the two contributions. It will be shown that the
bubble strength increases by about 2$\times$ and pump efficiency by about
3$\times$ between $30\,^{\circ}\\!{\rm C}$ and $70\,^{\circ}\\!{\rm C}$.

Finally, we show flow rate variation with resistor energy in Fig. 6. The data
show a clear saturation behavior. The flow rates initially increase but then
become independent of energy. It happens because of the vapor lock phenomenon.
At a critical energy flux, the entire fluid layer in direct contact with the
resistor vaporizes and isolates the hot resistor from the rest of the fluid.
Any further increase in the resistor temperature does not result in additional
heat transfer to the fluid. The pumping rate becomes independent of the
resistor temperature and released energy.

Figure 8: (a) Symbols: intrinsic bubble strength $q_{0}$ extracted from the
fits of Fig. 4, plotted vs. resistor energy and temperature. Lines: best fits
to the saturation model, Eq. (7). (b) Symbols: kick strength $q$ obtained by
fitting flow rates of Fig. 4 to the kick model, Eq. (5), plotted vs resistor
energy and temperature. Lines: best fits to the saturation model, Eq. (7). (c)
Mechanical efficiency of the inertial pump $\zeta=q/q_{0}$. The lines are
guides to the eye only.

##  IV  Analysis

###  IV.1  One-dimensional model of inertial pumping

To analyze the observed dependencies, we fitted the $L$-dependent flow rates
to the phenomenological 1d model of inertial pumping. Kornilovitch2013  Not
reproducing all model assumptions and equations here, it suffices to say that
three-dimensional nature of the channel flow is replaced by a one-dimensional
effective flow, both vapor-fluid interfaces are replaced by mathematical
points, and surface tension forces are neglected compared to inertial and
viscous forces. Boiling is assumed instantaneous, and resistor size is
neglected relative to the channel length. The model introduces one unknown
phenomenological parameter Kornilovitch2013  — “bubble strength” $q_{0}$ —
defined as mechanical momentum imparted on both sides of the channel at $t=0$
and measured in (kg m/s). By order of magnitude, $q_{0}$ is equal to the
product of peak vapor pressure in the bubble ($\sim 10^{6}$ Pa), time duration
of the high-pressure phase ($\sim 10^{-6}$ s) and channel cross-sectional area
($\sim 10^{-9}$ m2), providing an estimate $q_{0}=10^{-10}-10^{-9}$ kg m/s.
All other system’s parameters including the starting expansion point (middle
of the resistor), channel length, fluid viscosity (taken from tables), and
fluidic resistances of moving fluidic columns as functions of interface
positions, are all assumed to be independently known.

With only one unknown parameter, in theory a single flow rate measurement is
sufficient to determine $q_{0}$ from experimental data. In practice, this
approach leads to very noisy numbers as any experimental uncertainty directly
translates into a corresponding uncertainty in $q_{0}$. In this work, we
assume that drive vapor bubbles are local. That is, the intrinsic bubble
strength $q_{0}$ depends only on the resistor size and thermal properties of
the test fluid but is independent of the channel length and other parts of the
system away from the resistor. Flow rates in the eight channels should be
describable by one bubble strength. In other words, the eight flow rates
provide eight independent measurements of the same unknown parameter. By
fitting the entire $L$-dependence to the 1d model cancels experimental noise
of individual measurements and produces a more stable estimate of $q_{0}$.

Such $L$ fits are shown in Fig. 4 as dashed lines with best-fit $q_{0}$ values
listed in the legends. The same values are also plotted in Fig. 7 vs
temperature and in Fig. 8(a) vs resistor energy. In can be deduced from the
high-energy panels of Fig. 7 ($1.011-1.186$ $\mu$J) that the strength of well-
developed drive bubbles increases from about $2.5\times 10^{-10}$ kg m/s to
about $5.0\times 10^{-10}$ kg m/s between $30\,^{\circ}\\!{\rm C}$ and
$70\,^{\circ}\\!{\rm C}$. Thus, the bubble doubles in strength between the two
temperatures.

###  IV.2  Simplified kick model of inertial pumping and the pump mechanical
efficiency

The physical mechanism behind inertial pumping is complex. Torniainen2012 ;
Kornilovitch2013 ; Govyadinov2016  A key role is played by the reservoirs that
dissipate fluid’s momentum and provide quick refill on the shorter side of the
channel. Another important parameter is the asymmetry of resistor placement
within the channel. In general, total net flow consists of two contributions:
(i) primary pumping effect which is a shift of the bubble collapse point
relative to the starting expansion point, and (ii) secondary pumping effect
caused by unequal momenta of two colliding flows at bubble collapse and
resulting post-collapse flow. The full 3D CFD model and the 1d effective model
both capture all those effects well. Govyadinov2016

However, in effective inertial pumps the situation is simpler. In effective
pumps, the resistor is located much closer to one reservoir that to another.
The shorter side of the channel empties and refills very quickly so that
during that time the longer side barely moves due to its large mass and low
starting velocity. As a result, the primary pumping effect (shift of the
bubble collapse point relative to the resistor center) is relatively small. If
the channel also has a large cross section and low fluidic resistance, post-
collapse flow will persist for a long time, amplifying the secondary pumping
effect. Under these conditions, which are satisfied in our system, the initial
expansion-collapse phase can be neglected, in both time duration and pumping
effect, compared to the post-collapse flow. The situation is illustrated in
Fig. 9. The bubble collapses at 11.3 $\mu$s and the primary (linear) pumping
effect is 3.0 $\mu$m. After that, the post-collapse flow persists for another
400 $\mu$s resulting in a secondary displacement of 42 $\mu$m. One can see
that with a high degree of accuracy, the entire expansion-collapse phase can
be replaced by a nonzero mechanical momentum that the fluid receives when the
bubble collapses. Thus, a pump pulse can be thought to result from a kick the
fluid receives at time zero and then proceeds by inertia and decelerates by
friction. We will be calling this simplified picture of inertial pumping the
kick model. Its main advantage is extreme mathematical simplicity. Indeed, the
post-collapse flow is one with a constant fluid mass (there is no vapor bubble
anymore), so all the complications and nonlinearities associated with variable
mass motion Kornilovitch2013  disappear. We now derive the main equation of
the kick model.

Figure 9: Solutions of the one-dimensional inertial pumping equations
Kornilovitch2013  for loop 2 and $T=70\,^{\circ}\\!{\rm C}$. Notice how the
primary priming effect — shift of the bubble collapse point relative to the
starting expansion point ($x=63$ $\mu$m vs 60 $\mu$m) — is dominated by the
post-collapse flow. The overall pump dynamics can be accurately described by
Poiseuille flow after an initial kick.

Let $q$ be the mechanical momentum imparted on the entire fluid at bubble
collapse. The fluid’s mass is $M=\rho AL$, so the initial velocity is
$v_{0}=q/M$. The friction force acting on the fluid is $\kappa\eta vL$, where
$\kappa$ is the dimensionless geometric friction factor, see Appendix, Eqs.
(32) and (33). $\kappa_{\rm cyl}=8\pi$ for cylindrical channels, $\kappa_{\rm
sq}=28.45$ for 1:1 square channels, and $\kappa_{\rm sq}=28.47$ for 30:31
channels studied in this work. Neglecting the Bernoulli correction,
Kornilovitch2013  the equation of motion reads

$\rho AL\ddot{x}+\kappa\eta L\dot{x}=0\>,$ (2)

where $x(t)$ is the position of a reference point and the dot denotes a
derivative over time. Solving for the velocity and displacement one obtains

$\displaystyle v(t)$ $\displaystyle=$ $\displaystyle v_{0}\,e^{-t/\tau}\>,$
(3) $\displaystyle x(t)-x_{0}$ $\displaystyle=$ $\displaystyle
v_{0}\tau\left(1-e^{-t/\tau}\right)\>,$ (4)

where $\tau=(\rho A)/(\kappa\eta)$ is the characteristic time constant of the
post-collapse flow. In our system, $\tau\sim 50-100$ $\mu$s, depending on
temperature. Thus, the net flow at the end of a pulse is simply

$Q=\triangle x(\infty)A=v_{0}\tau A=\frac{qA}{\kappa\eta L}\>.$ (5)

This formula embodies a very transparent physics: The flow rate is kick
momentum per unit cross-sectional area $q/A$ divided by the channel’s fluidic
resistance $\kappa\eta L/A^{2}$. Note that the kick model predicts a precise
$1/L$ length dependence of flow rates, which is consistent with the
experimental data of Fig. 4.

Kick strength $q$ has a similar physical meaning as bubble strength $q_{0}$.
The difference is that $q_{0}$ is isotropic momentum with which the original
drive bubble pushes two columns of fluid in opposite directions, whereas $q$
is a unidirectional momentum acquired by the entire fluid as result of
inelastic collision of two flows at bubble collapse. Since both $q$ and
$q_{0}$ have the same physical units, their dimensionless ratio can be
interpreted as the pump’s mechanical efficiency $\zeta$,

$\zeta=\frac{q}{q_{0}}\>,\hskip 28.45274pt0\leq\zeta\leq 1\>.$ (6)

$\zeta$ describes how isotropic momentum of an original drive bubble is
converted to “useful” unidirectional kick momentum at bubble collapse.

To determine $q$ and $\zeta$, we fitted the experimental flow rates of Fig. 4
to the kick formula, Eq. (5). Best-fit $q$ values are plotted in Fig. 7
alongside $q_{0}$. Focusing again on well-developed bubbles ($1.011-1.186$
$\mu$J), the kick strength grows from about $0.3\times 10^{-10}$ to $2.0\times
10^{-10}$ kg m/s between $30\,^{\circ}\\!{\rm C}$ and $70\,^{\circ}\\!{\rm
C}$. Pump efficiency $\zeta$ increases by about 3$\times$ over the same
temperature interval. The kick strength and efficiency are also plotted vs.
resistor energy in Figs. 8(b) and 8(c).

###  IV.3  Threshold model of energy dependence

Next, we discuss the resistor energy dependence of bubble strength $q_{0}$ and
kick strength $q$. Experimental data are shown in Fig. 8. Both quantities
display a clear saturation behavior similar to that of flow rates, cf. Fig. 6.
For quantitative analysis, we fitted $q_{0}(E)$ and $q(E)$ to a three-
parameter phenomenological threshold model

$q_{0}(E),\;q(E)=S\cdot\tanh{\left(\frac{E-E_{T}}{E_{w}}\right)}.$ (7)

Best-fit values of scale $S$, threshold energy $E_{T}$ and transition width
$E_{w}$ are listed in Table 2. The resulting models are plotted in Figs. 8(a)
and 8(b) on top of experimental data.

Quantity | T, ∘C |  $S$, e-10 kg$\cdot$m/s |  $E_{T}$, $\mu$J |  $E_{w}$, $\mu$J
---|---|---|---|---
$q_{0}$ | 30 | 2.548 | 0.935 | 0.0609
40 | 3.288 | 0.911 | 0.0747
50 | 3.908 | 0.892 | 0.0767
60 | 4.375 | 0.871 | 0.0741
70 | 4.943 | 0.852 | 0.0567
$q$ | 30 | 0.355 | 0.956 | 0.0831
40 | 0.742 | 0.937 | 0.1006
50 | 1.160 | 0.922 | 0.0870
60 | 1.589 | 0.903 | 0.0818
70 | 2.170 | 0.878 | 0.0575
Table 2:  Parameters of the threshold model, Eq. (7), for different
temperatures.  Figure 10: Parameters of the threshold model, Eq. (7) and Table
2, plotted vs operating pump temperature. The $T=40\,^{\circ}\\!{\rm C}$ point
was excluded from fitting $E_{W}(T)$ for the kick strength $q$.

Finally, we discuss variation of $S$, $E_{T}$, and $E_{w}$ with temperature.
The values of Table 2 are plotted in Fig. 10. The scale parameter $S$ shows a
clear linear growth with $T$. This is a manifestation of the bubble getting
stronger with temperature, as discussed earlier. The threshold energy $E_{T}$
also shows a linear temperature dependence but with a negative slope. At
higher $T$, the drive bubbles reach a saturated state at an earlier energy
because less energy is spent on warming the fluid before boiling. This
explains the observed trend. Finally, the width parameter $E_{w}$ shows a non-
monotonic behavior with a maximum near $T=50\,^{\circ}\\!{\rm C}$. We do not
have an explanation for such a behavior.

The simplicity of the kick model allows for construction of a complete
predictive flow rate model for our microfluidic system. Combining the best
linear fits for $S(T)$ and $E_{T}(T)$, the best quadratic fit for $E_{w}(T)$,
cf. Fig. 10, the energy model for the kick strength $q(E)$, Eq. (7), and the
kick formula, Eq. (5), one obtains:

$Q(T,E,L)=\frac{S(T)A}{\kappa\eta(T)L}\tanh{\left(\frac{E-E_{T}(T)}{E_{w}(T)}\right)},$
(8)

where

$\displaystyle S(T)$ $\displaystyle=$ $\displaystyle(0.0448\,T-1.0351)\times
10^{-10}\;{\rm kg}\cdot{\rm m/s}\>,$ (9) $\displaystyle E_{T}(T)$
$\displaystyle=$ $\displaystyle(-0.0019\,T+1.0148)\;\mu{\rm J}\>,$ (10)
$\displaystyle E_{w}(T)$ $\displaystyle=$
$\displaystyle[-0.00005\,(T-50)^{2}-0.0006\,(T-50)$ (11)
$\displaystyle+0.0894]\;\mu{\rm J}\>,$

and temperature $T$ is specified in degrees ${}^{\circ}\\!{\rm C}$. Note,
however, that this model is specific to the particular resistors size, channel
cross-sectional area and test fluid investigated in this work. For a different
system, the model would have to be rebuilt using the methodology developed
herein.

##  V  Summary and conclusions

In this paper, we systematically investigated inertial pump flow rates versus
channel length $L$, operating temperature $T$, and TIJ resistor thermal
release energy $E$. We developed a robust methodology of measuring flow rates
based on identifying the fastest tracer particle and theoretical flow profiles
in rectangular channels, see Sec. II.3 and the Appendix.

Our experimental results are summarized in Figs. 4-6. The length dependence of
flow rates follows a close-to-perfect $1/L$ behavior, which is consistent with
the dominant contribution of post-collapse phase to overall net flow. The most
interesting findings are associated with temperature variation. We determined
that flow rates delivered by the same pump can grow by as much as a factor of
12 between $T=30\,^{\circ}\\!{\rm C}$ and $T=70\,^{\circ}\\!{\rm C}$. Only
partially such an increase can be accounted for by a lower viscosity while the
rest must be attributed to the pump getting stronger with temperature. Through
the phenomenological analysis of Sec. IV.1, we established that the drive
bubble grows stronger by about 2$\times$ in the same temperature interval. The
remaining contribution comes from an increasing mechanical efficiency of the
pump. The latter term was introduced in Sec. IV.2 as a consequence of a newly
developed kick model of inertial pumping. It is possible that the efficiency
enhancement is caused by a dropping viscosity in its own way, making the
overall viscosity impact on flow rates nonlinear. This topic deserves further
analysis.

We also observed a clear saturation behavior in the resistor energy dependence
of flow rates. Saturation can be understood based on vapor lock effect that
limits heat flux from a TIJ resistor beyond a certain threshold. Through
phenomenological modeling of Sec. IV.3, we determined threshold and transition
width parameters and established their variation with temperature. Finally, we
constructed a fully predictive model of our inertial pumps in Eqs. (8)-(11).

Integrated inertial pumps are expected to become key drivers of complex
microfluidic chips in the future. Hayes2018 ; Govyadinov2015  Temperature
variations are often an essential part of the workflows, for example, in chip-
scale Polymerase Chain Reaction thermocyclers. Our results demonstrate that
such system must be designed with care, as pumping rates might vary
significantly depending on what part of the workflow the device is executing.
Additionally, they point to the need of accurate flowmeters and temperature
sensors that could ensure stable flow rates via close-loop control systems.
Based on our results, we also expect significant temperature dependence in
other microfluidic operations enabled by TIJ resistors, for example in active
mixing. Hayes2018b

###### Acknowledgements.

The authors wish to thank Carson Denison, Mandana Hamidi Haines, Brandon
Hayes, Anand Jebakumar, Richard Seaver, Erik Torniainen for discussions on the
subject of this paper as well as Arun Agarwal, Paul Benning and Norman
Pawlowski for supporting this work.

##  Appendix A  Poiseuille flow in rectangular microchannels

Figure 11: (a) Maximal fluid velocity in the center of a rectangular channel
as a function of the aspect ratio. Specific values of $v_{\rm max}/\langle
v\rangle$ for (1:1) $=2.096$; (2:1) $=1.992$; (3:1) $=1.864$; and (10:1)
$=1.601$. (b) Velocity profile in channels with square cross-sections, $a=b$.
The red dashed line is a parabola. The thick blue line is the exact
theoretical formula, Eq. (28). (c) Velocity profile in a wide-and-thin
channel, $a=4b$. The red dashed line is a parabola. The thick blue line is the
exact theoretical formula, Eq. (28). Notice significant deviation from the
parabolic shape and flattening of the profile in the center of the channel.

For completeness, we derive here an analytical formula for the flow profile in
rectangular channels. Pressure-driven Poiseuille flow field $v_{y}(x,z)$ in a
channel of constant cross section satisfies a two-dimensional Poisson equation
Landau_v6 ; Bruus2008

$\frac{\partial^{2}v_{y}}{\partial x^{2}}+\frac{\partial^{2}v_{y}}{\partial
z^{2}}=-\frac{\triangle p}{\eta L}\>,$ (12)

where $L$ is the channel length and $\eta$ is fluid’s dynamic viscosity. The
channel is oriented along the $y$ axis and the cross section lies within the
$xz$ plane. $v_{y}=0$ at the channel walls due to the no-slip condition. Note
that the problem is mathematically equivalent to torsion of thin long elastic
beams. Timoshenko1970  Consider a rectangular channel of width $a$, $0\leq
x\leq a$, and height $b$, $0\leq z\leq b$. Separating variables, the solution
of (12) is given by the double series:

$v_{y}(x,z)=\frac{\triangle p}{\eta
L}\frac{16}{ab}\sum^{\infty}_{n,m=0}\frac{\sin(p_{2n+1}x)\sin(q_{2m+1}z)}{p_{2n+1}q_{2m+1}(p^{2}_{2n+1}+q^{2}_{2m+1})}\>,$
(13) $p_{2n+1}=\frac{\pi(2n+1)}{a},\hskip
14.22636ptq_{2m+1}=\frac{\pi(2m+1)}{b}\>.$ (14)

Integral of $v_{y}(x,z)$ over the cross section is equal to volumetric flow
rate $Q$, which in turn defines average velocity $\langle v\rangle$ as

$\langle v\rangle\equiv\frac{Q}{ab}=\frac{1}{ab}\int dx\,dz\>v_{y}(x,z)\>.$
(15)

Performing integration and partial summation Prudnikov1986  of Eq. (13), the
pressure gradient can be eliminated in favor of $\langle v\rangle$:

$v_{y}=\langle
v\rangle\frac{ab}{4S_{1}}\sum^{\infty}_{n,m=0}\frac{\sin(p_{2n+1}x)\sin(q_{2m+1}z)}{p_{2n+1}q_{2m+1}(p^{2}_{2n+1}+q^{2}_{2m+1})}\>.$
(16)

Here $S_{1}$ is a single series that can be expressed in two equivalent forms

$\displaystyle S_{1}(a,b)$ $\displaystyle=$
$\displaystyle\frac{a^{2}b^{4}}{3\cdot
2^{8}}\left\\{1-\frac{192}{\pi^{5}}\left(\frac{b}{a}\right)\cdot
S_{3m}(a,b)\right\\}$ (17) $\displaystyle=$
$\displaystyle\frac{a^{4}b^{2}}{3\cdot
2^{8}}\left\\{1-\frac{192}{\pi^{5}}\left(\frac{a}{b}\right)\cdot
S_{3n}(a,b)\right\\},$ (18)
$S_{3m}(a,b)=\sum^{\infty}_{m=0}\frac{1}{(2m+1)^{5}}\tanh{\frac{\pi
a(2m+1)}{2b}}\>,$ (19)
$S_{3n}(a,b)=\sum^{\infty}_{n=0}\frac{1}{(2n+1)^{5}}\tanh{\frac{\pi
b(2n+1)}{2a}}\>.$ (20)

The “flow profile” method of Sec. II.3 requires knowledge of the maximal fluid
velocity for a given $x$. It is clear from symmetry considerations that a
maximum velocity is reached at $z=b/2$ for any $x$. Setting $z=b/2$ in Eq.
(16) yields

$v_{\rm m}(x)=\langle
v\rangle\frac{ab}{4S_{1}}\sum^{\infty}_{n,m=0}\frac{\sin(p_{2n+1}x)(-1)^{m}}{p_{2n+1}q_{2m+1}(p^{2}_{2n+1}+q^{2}_{2m+1})}\>.$
(21)

The sum over $m$ can be calculated analytically Prudnikov1986

$v_{\rm m}(x)=\langle
v\rangle\frac{ab^{2}}{16S_{1}}\sum^{\infty}_{n=0}\frac{\sin(p_{2n+1}x)}{p^{3}_{2n+1}}\left(1-\frac{1}{\cosh{\frac{p_{2n+1}b}{2}}}\right).$
(22)

Out of the remaining two sums, the first one can also be calculated
analytically, Prudnikov1986  yielding the flow profile

$v_{\rm m}(x)=\langle
v\rangle\frac{a^{2}b^{2}}{2^{7}S_{1}}\left\\{x(a-x)-\frac{8}{a}\sum^{\infty}_{n=0}\frac{\sin(p_{2n+1}x)}{p^{3}_{2n+1}\cosh{\frac{p_{2n+1}b}{2}}}\right\\}\\!.$
(23)

Next, we calculate the absolute maximum velocity, which is reached in the
middle of a channel. Setting $x=a/2$ in Eq. (23), one obtains

$v_{\rm max}(a,b)=\langle
v\rangle\cdot\frac{3}{2}\frac{1-\frac{32}{\pi^{3}}\cdot
S_{4n}(a,b)}{1-\frac{192}{\pi^{5}}\left(\frac{a}{b}\right)\cdot
S_{3n}(a,b)}\>,$ (24)

where

$S_{4n}(a,b)\equiv\sum^{\infty}_{n=0}\frac{(-1)^{n}}{(2n+1)^{3}\cosh{\frac{\pi
b(2n+1)}{2a}}}\>,$ (25)

and $S_{3n}$ is defined in Eq. (20). Note that $v_{\rm max}$ can also be
expressed in an equivalent form

$v_{\rm max}(a,b)=\langle
v\rangle\cdot\frac{3}{2}\frac{1-\frac{32}{\pi^{3}}\cdot
S_{4m}(a,b)}{1-\frac{192}{\pi^{5}}\left(\frac{b}{a}\right)\cdot
S_{3m}(a,b)}\>,$ (26)
$S_{4m}(a,b)=\sum^{\infty}_{m=0}\frac{(-1)^{m}}{(2m+1)^{3}\cosh{\frac{\pi
a(2m+1)}{2b}}}\>,$ (27)

which reflects the fact that the maximum-to-average velocity ratio is
invariant under transformation $a\leftrightarrow b$ that leaves the channel’s
aspect ratio unchanged. It is possible to show that in both limits $a\gg b$
and $a\ll b$, $v_{\rm max}\rightarrow\frac{3}{2}\langle v\rangle$, the value
for an infinite slot. The full Eqs. (24) and (26) vs aspect ratio $(a/b)$ is
shown in Fig. 11(a).

Returning now to the task of calculating the velocity profile, divide Eq. (23)
by Eq. (24) to get, after transformations

$v_{\rm max}(x)=v_{\rm
max}\cdot\frac{4\left(\frac{x}{a}\right)\left(1-\frac{x}{a}\right)-\frac{32}{\pi^{3}}\cdot
S_{6}(x)}{1-\frac{32}{\pi^{3}}\cdot S_{4n}}\>,$ (28)
$S_{6}(x;a,b)\equiv\sum^{\infty}_{n=0}\frac{\sin{\frac{\pi
x}{a}(2n+1)}}{(2n+1)^{3}\cosh{\frac{\pi b(2n+1)}{2a}}}\>.$ (29)

The flow profile, Eq. (28), is the sum of a parabola plus a linear combination
of sine harmonics. The relative contribution of the sines depends on the
channel’s aspect ratio. In tall and skinny channels, $b\gg a$, the sines’
weights are exponentially small, and the flow profile is a parabola. As $a/b$
increases, the sine weights grow to be non-negligible and in square channels,
$a=b$, deviation from a pure parabola is already noticeable although it
remains small, see Fig. 11(b). However, in thin and wide channels, $a\gg b$,
many harmonics contribute rendering the flow profile substantially non-
parabolic, as illustrated in Fig. 11(c) for the case of $a/b=4$.

Calculating the viscous stress tensor

$\sigma_{ik}=\eta\left(\frac{\partial v_{i}}{\partial x_{k}}+\frac{\partial
v_{k}}{\partial x_{i}}\right),$ (30)

and integrating over the perimeter of the channel cross section, the total
friction force can be expressed as $F=\kappa\cdot\eta\langle v\rangle L$,
where $\eta$ is the dynamic viscosity and $L$ is the channel length. The
dimensionless coefficient $\kappa$ is

$\kappa_{\rm rectangle}=\frac{a^{3}b^{3}}{64\>S_{1}(a,b)}\>.$ (31)

Using Eqs. (19) and (20) yields the final result

$\displaystyle\kappa_{\rm rectangle}$ $\displaystyle=$
$\displaystyle\frac{12a}{b}\frac{1}{1-\frac{192}{\pi^{5}}\left(\frac{b}{a}\right)\cdot
S_{3m}(a,b)}\>,$ (32) $\displaystyle=$
$\displaystyle\frac{12b}{a}\frac{1}{1-\frac{192}{\pi^{5}}\left(\frac{a}{b}\right)\cdot
S_{3n}(a,b)}\>.$ (33)

For the channels studied in this work, $a/b=30/31$ and $\kappa=28.47$.

## References

  * (1) H. Yuan and A. Prosperetti, “The pumping effect of growing and collapsing bubbles in a tube,” J. Micromech. Microeng. 9, 402–413 (1999).
  * (2) Z. Yin and A. Prosperetti, “A microfluidic ‘blinking bubble’ pump,” J. Micromech. Microeng. 15, 643–651 (2005).
  * (3) Z. Yin and A. Prosperetti, “ ‘Blinking bubble’ micropump with microfabricated heaters,” J. Micromech. Microeng. 15, 1683–1691 (2005).
  * (4) E.D. Torniainen, A.N. Govyadinov, D.P. Markel, and P.E. Kornilovitch, “Bubble–driven inertial micropump,” Phys. Fluids 24, 122003 (2012).
  * (5) A.N. Govyadinov, P.E. Kornilovitch, D.P. Markel, and E.D. Torniainen, “Single–pulse dynamics and flow rates of inertial micropumps,” Microfluid. Nanofluid. 20, 73 (2016).
  * (6) J. Stasiak, S. Richards, and P. Benning, “Hewlett-Packard’s MEMS technology — Thermal inkjet printing and beyond,” in Microelectronics to Nanoelectronics: Materials, Devices & Manufacturability, edited by A. B. Kaul (CRC/Taylor & Francis, Boca Raton, FL, 2012).
  * (7) P.E. Kornilovitch, A.N. Govyadinov, D.P. Markel, and E.D. Torniainen, “One–dimensional model of inertial pumping,” Phys. Rev. E 87, 023012 (2013).
  * (8) B. Hayes, G.L. Whiting, and R. MacCurdy, “Modeling of contactless bubble–bubble interactions in microchannels with integrated inertial pumps,” Phys. Fluids 33, 042002 (2021).
  * (9) B.S. Hayes, A.N. Govyadinov, and P.E. Kornilovitch, “Microfluidic switchboards with integrated inertial pumps,” Microfluid. Nanofluid. 22, 15 (2018).
  * (10) G.R. Wang, J.G. Santiago, M.G. Mungal, B. Young, and S. Papandemetriou, “A laser induced cavitation pump,” J. Micromech. Microeng. 14, 1037–1046 (2004).
  * (11) R.J. Dijkink and C.D. Ohl, “Laser-induced cavitation based micropump,” Lab Chip 8, 1676–1681 (2008).
  * (12) J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S Saalfeld, B.Schmid, J.-Y. Tinevez, D.J. White, V. Hartenstein, K. Eliceiri, P. Tomancak, and A. Cardona, “Fiji: An open-source platform for biological-image analysis,” Nature Methods 9, 676–682 (2012).
  * (13) J.-Y. Tinevez, N. Perry, J. Schindelin, G.M. Hoopes, G.D. Reynolds, E. Laplantine, S.Y. Bednarek, S.L. Shorte, and K.W. Eliceiri, “TrackMate: An open and extensible platform for single-particle tracking,” Methods 115, 80–90 (2017).
  * (14) A.N. Govyadinov, E.D. Torniainen, P.E. Kornilovitch, and D.P. Markel, “Path to low cost microfluidics,” in Proceedings of NIP31, Portland, OR (IS&T 2015), pp. 490–497.
  * (15) See https://wiki.anton-paar.com/us-en/water/ for tables of water viscosity and density vs temperature.
  * (16) B. Hayes, A. Hayes, M. Rolleston, A. Ferreira, and J. Kirsher, “Pulsatory mixing of laminar flow using bubble-driven micro-pumps,” in Proceedings of the ASME 2018 International Mechanical Engineering Congress and Exposition (ASME 2018), Vol. 7.
  * (17) L.D. Landau and E.M. Lifshits, Course of Theoretical Physics, Volume 6: Fluid Mechanics, 2nd ed. (Elsevier Ltd., 1987).
  * (18) H. Bruus, Theoretical Microfluidics (Oxford University Press, 2008), Chap. 3.
  * (19) S.P. Timoshenko and J.N. Goodier, Theory of Elasticity, 3rd ed. (McGaw-Hill, Inc., 1970).
  * (20) A.P. Prudnikov, Yu.A. Brychkov, and O.I. Marichev, Integrals and Series, Volume 1: Elementary Functions (Taylor & Francis), Secs. 5.1 and 5.4.


††thanks: Current affiliation: National Research Council of Canada, 100 Sussex
Drive, K1N 5A2, Ottawa, Ontario, Canada

# A genetic algorithm for the response of
twisted nematic liquid crystals to an applied field

Alicia Sit  Nexus for Quantum Technologies, University of Ottawa, K1N 5N6,
Ottawa, ON, Canada  Francesco Di Colandrea<EMAIL_ADDRESS>Nexus for Quantum Technologies, University of Ottawa, K1N 5N6, Ottawa, ON,
Canada  Alessio D’Errico  Nexus for Quantum Technologies, University of
Ottawa, K1N 5N6, Ottawa, ON, Canada  Ebrahim Karimi  Nexus for Quantum
Technologies, University of Ottawa, K1N 5N6, Ottawa, ON, Canada

###### Abstract

When an external field is applied across a liquid-crystal cell, the twist and
tilt distributions cannot be calculated analytically and must be extracted
numerically. In the standard approach, the Euler-Lagrange equations are
derived from the minimization of the free energy of the system and then solved
via finite-difference methods, often implemented in commercial software. These
tools iterate from initial solutions that are compatible with the boundary
conditions, providing limited to no flexibility for customization. Here, we
present a genetic algorithm that outputs fast and accurate solutions to the
integral form of the equations. In our approach, the evolutionary routine is
sequentially applied at each position within the bulk of the cell, thus
overcoming the necessity of assuming trial solutions. The predictions of our
routine strongly support the experimental observations on different instances
of spatially varying twisted nematic liquid-crystal cells, patterned with
different topologies on the two alignment layers.

## Introduction

The response of liquid crystals to an applied field, whether magnetic Leslie
(1970) or electric Deuling (1974), can be modeled within the elastic continuum
theory De Gennes and Prost (1993). This predicts that the molecules tend to
align with the field direction. The equilibrium configuration can be found by
imposing the minimization of the total free energy. At the equilibrium, the
twist and tilt distributions of the liquid-crystal director thus satisfy a set
of Euler-Lagrange equations Deuling (1974); Preist  _et al._ (1989). These are
typically solved via finite-difference or finite-element methods Yang and Wu
(2014); Moser  _et al._ (2019), mostly embedded in modeling software such as
LC3D Wang  _et al._ (2004, 2005), DIMOS Peng  _et al._ (2016); Lee  _et al._
(2017) and COMSOL Guirado  _et al._ (2022); Gorkunov  _et al._ (2022); Sova
_et al._ (2023). Starting from trial distributions that satisfy the boundary
conditions, these routines converge to quasi-optimal solutions within a range
of iterations.

A recent work explored the application of genetic algorithms to determine the
director distribution in a few relevant cases, including hybrid and twisted
nematic liquid-crystal cells Yang and Collings (2020). Inspired by Darwin’s
evolutionary theory, the basic idea is to start from a population of random
guesses (individuals) and let the workflow of the genetic algorithm (GA)
select the individuals that better approximate the optimal solution of the
physical problem Holland (1992), formulated in terms of the minimization of a
cost function (fitness). The workflow is based on the application of an
iterative scheme composed of three sequential operators: selection, crossover,
and mutation. The joint usage of these operators allows the algorithm to
evolve a population of candidate solutions toward quasi-optimal solutions. A
selection operator picks the “best” individuals in the population, forming in
this way the mating pool. The individuals in the mating pool generate new
possible solutions, forming the offspring, through crossover and mutation
operators, which respectively emulate the recombination and the mutation of
individuals in a natural environment. In the end, the offspring replaces the
initial population of the algorithm, and a new iteration (generation) can
start. The algorithm ends when a certain termination criterion is satisfied
Yao (1993), for instance, when a maximum number of generations is reached or
when no significant variation of the fitness is detected between successive
iterations.

A natural choice for the cost function is the total free energy Yang and
Collings (2020). However, this still requires evolving populations of entire
twist and tilt distributions, parametrized in terms of the Cartesian
components of the director at various locations between the substrates. This
is essentially needed to approximate derivatives with finite differences,
wherein the number of sampled points can severely affect the performance of
the algorithm. Here, instead, we devise a GA to directly solve the bulk-
integral version of the Euler-Lagrange equations Deuling (1974); Preist  _et
al._ (1989) in correspondence with a discrete set of positions within the
cell, processed in a sequential routine. Remarkably, this approach allows us
to model the response of the liquid-crystal layer to arbitrary field
strengths.

In this paper, we focus on twisted nematic liquid-crystal (TNLC) cells, which
feature complex three-dimensional modulations of the molecular director;
however, the presented algorithm is not limited to TNLC configurations. First,
we review the theory of TNLC cells in the presence of an external electric
field, reporting the set of equations from which the twist and tilt
distributions can be retrieved. Then, we provide technical details of our GA
implementation and compute numerical results in a variety of settings. The
solutions are also used to infer semi-analytical expressions that can provide
good approximations for quantitative estimates. Finally, the performance of
our routine is experimentally validated on spatially varying TNLC cells,
specifically dual-$q$-plates Sit  _et al._ (2023), fabricated by designing
different topological charges on the two alignment layers.

## Theory

The free energy density of the system, given by the fundamental elastic
continuum equation, is De Gennes and Prost (1993),

$\begin{split}\mathcal{F}=\frac{1}{2}\\{&K_{1}[\divergence\mathbf{\hat{n}}(\mathbf{r})]^{2}+K_{2}[\mathbf{\hat{n}}(\mathbf{r})\cdot\curl\mathbf{\hat{n}}(\mathbf{r})]^{2}\\\
&+K_{3}[\mathbf{\hat{n}}(\mathbf{r})\times\curl\mathbf{\hat{n}}(\mathbf{r})]^{2}-\mathbf{D}\cdot\mathbf{E}-\mathbf{B}\cdot\mathbf{H}\\},\end{split}$
(1)

where the $K_{i}$ are the elastic constants for the splay, twist, and bend
distortions, respectively. The local director distribution
$\mathbf{\hat{n}}(\mathbf{r})$ is our object of interest, and can be expressed
in spherical coordinates,

$\displaystyle n_{x}$ $\displaystyle=$
$\displaystyle\cos{\phi(z)}\cos\theta(z),$ $\displaystyle n_{y}$
$\displaystyle=$ $\displaystyle\sin{\phi(z)}\cos\theta(z),$ (2) $\displaystyle
n_{z}$ $\displaystyle=$ $\displaystyle\sin\theta(z),$

where $\phi(z)$ and $\theta(z)$ are the twist and tilt distributions,
respectively, that minimize the free energy of the system. We limit the
following discussion to the TNLC geometry. For the field-free case, one
obtains a linear twist distribution ${\phi(z)=\alpha z/L}$, where $L$ is the
thickness of the cell, and ${\theta(z)=0}$ Yamauchi (2005). If an electric
field is applied in the direction normal to the cell,
$\mathbf{E}=E\mathbf{\hat{z}}$—which is equivalent to applying a voltage
$V=EL$ across the twisted cell—the resulting $\phi(z)$ and $\theta(z)$ must be
numerically found through a coupled set of integrals Deuling (1974). The twist
distribution $\phi(z)$ is given by,

$\phi(z)=\beta\int_{0}^{\theta(z)}\frac{\sqrt{1+\kappa\sin^{2}\theta}}{g(\theta)\cos^{2}\theta(1+\tau\sin^{2}\theta)}\text{d}\theta,$
(3)

where $\kappa=(K_{3}-K_{1})/K_{1}$, $\tau=(K_{3}-K_{2})/K_{2}$, and $\beta$ is
an integration parameter to be determined from the boundary conditions. The
tilt distribution $\theta(z)$ is determined implicitly through,

$\frac{z}{L}=\frac{1}{2}\int_{0}^{\theta(z)}\frac{\sqrt{1+\kappa\sin^{2}\theta}}{g(\theta)}\text{d}\theta\bigg{/}\int_{0}^{\theta_{m}}\frac{\sqrt{1+\kappa\sin^{2}\theta}}{g(\theta)}\text{d}\theta,$
(4)

where $\theta_{m}\equiv\theta_{m}(V)$ is the maximum tilt angle located at
$z=L/2$ and is itself a function of the applied field. The $g(\theta)$
function depends on $\beta$ and $\theta_{m}$, with the form,

$\begin{split}&g(\theta)=\bigg{\\{}\frac{\sin^{2}\theta_{m}-\sin^{2}\theta}{(1+\gamma\sin^{2}\theta)(1+\gamma\sin^{2}\theta_{m})}+\beta^{2}\frac{1+\kappa}{1+\tau}\times\\\
&\biggl{(}\frac{1}{(1+\tau\sin^{2}\theta_{m})\cos^{2}\theta_{m}}-\frac{1}{(1+\tau\sin^{2}\theta)\cos^{2}\theta}\biggl{)}\bigg{\\}}^{1/2},\end{split}$
(5)

where $\gamma=(\epsilon_{\parallel}-\epsilon_{\perp})/\epsilon_{\perp}$. Here,
$\epsilon_{\parallel}$ ($\epsilon_{\perp}$) denotes the dielectric constant
per unit volume that is parallel (perpendicular) to the local director. The
maximum tilt angle $\theta_{m}$ is found from,

$\frac{V}{V_{T0}}=\frac{2}{\pi}\int_{0}^{\theta_{m}}\frac{\sqrt{1+\kappa\sin^{2}\theta}}{(1+\gamma\sin^{2}\theta)g(\theta)}\text{d}\theta,$
(6)

where ${V_{T0}=\pi\sqrt{K_{1}/(\epsilon_{0}\Delta\epsilon)}}$, with
${\Delta\epsilon=\epsilon_{\parallel}-\epsilon_{\perp}}$, is the threshold
voltage for ${\phi(z)=0}$, when the Fréedericksz transition occurs in the
zero-twist configuration Deuling (1972). The integration parameter $\beta$ can
be found by evaluating Eq. (3) at $z=L/2$:

$\phi(L/2)=\frac{\phi_{m}}{2}=\beta\int_{0}^{\theta_{m}}\frac{\sqrt{1+\kappa\sin^{2}\theta}}{g(\theta)\cos^{2}\theta(1+\tau\sin^{2}\theta)}\text{d}\theta,$
(7)

where $\phi_{m}$ is the maximum twist angle located at the back plate
(${z=L}$) of the cell.

Equations (6)-(7) are a coupled set of integrals to be solved simultaneously.
The numerical integration of this system is not trivial due to the strongly
singular behaviour of the integrand functions. A better way to solve these
equations is to do it iteratively, by first setting ${\beta=0}$ in Eqs.
(5)-(6), which produces a first estimate for
${\theta_{m0}=\theta_{m}(\beta=0)}$. This can be used to determine
${\beta_{0}=\beta(\theta_{m0})}$ from Eq. (7), and so on until the desired
convergence is achieved. Once the parameters $\beta$ and $\theta_{m}$ have
been determined, we need to solve Eqs. (3)-(4) to eventually extract the twist
and tilt distributions. We have opted for an optimization approach based on
evolutionary methods, specifically a genetic algorithm.

The following two analytical expressions Deuling (1974) can be used as
indicators for whether our GA converges successfully. The threshold voltage
$V_{T}$ for a given total twist angle $\phi_{m}$ is,

$V_{T}(\phi_{m})=V_{T0}\left[1+\left(\frac{\phi_{m}}{\pi}\right)^{2}\left(\frac{K_{3}}{K_{1}}-2\frac{K_{2}}{K_{1}}\right)\right]^{1/2},$
(8)

and the value of the parameter $\beta_{T}$ at the threshold voltage for a
given $\phi_{m}$ is,

$\beta_{T}(\phi_{m})=\left[\left(\frac{\pi}{\phi_{m}}\right)^{2}+\frac{K_{3}}{K_{1}}-2\frac{K_{2}}{K_{1}}\right]^{-1/2}.$
(9)

These relations are useful as they describe limiting behaviours when
${V\rightarrow V_{T}}$, which our method must be able to reproduce.

## Numerical optimization

To approximately solve Eqs. (6)-(7), we implement two nested genetic
algorithms which evolve real-valued individuals $\theta_{mi}$ and $\beta_{i}$.
As prescribed by the iterative method mentioned above, the first GA is run to
determine $\theta_{m0}$, corresponding to the initial guess for ${\beta=0}$.
The second GA is then run to determine $\beta_{0}$, assuming
${\theta_{m}=\theta_{m0}}$, and so on. Each individual is a candidate to
provide an optimal approximation to the actual solutions $\theta_{m}$ and
$\beta$. By means of operators mimicking the natural selection mechanism, the
GAs select for reproduction those individuals which better minimize the
following cost functions:

$\mathcal{L}_{\theta_{m}}=\left|\frac{V}{V_{T0}}-\frac{2}{\pi}\int_{0}^{\theta_{m}}\frac{\sqrt{1+\kappa\sin^{2}\theta}}{(1+\gamma\sin^{2}\theta)g(\theta)}\text{d}\theta\right|^{2},$
(10)
$\mathcal{L}_{\beta}=\left|\phi_{m}-2\beta\int_{0}^{\theta_{m}}\frac{\sqrt{1+\kappa\sin^{2}\theta}}{g(\theta)\cos^{2}\theta(1+\tau\sin^{2}\theta)}\text{d}\theta\right|^{2},$
(11)

within the current generation. The numerical integrations required to evaluate
these cost functions are successfully performed using the IMT-rule Iri  _et
al._ (1970) to handle singularities in finite integration regions Davis and
Rabinowitz (2014).

The detailed sequence of operators used in our algorithm will now be
described. First, the well-known tournament selection mechanism Goldberg and
Deb (1991) is used as a selection operator. This consists of repeating the
following steps $N$ times, where $N$ is the population size:

  1. 1.

Randomly select a subset of $k$ individuals.

  2. 2.

Choose the best individual in the subset to be inserted in the mating pool.

For our purposes, the best are those individuals that better minimize the cost
functions in Eqs. (10)-(11). The blend crossover Eshelman and Schaffer (1993)
is applied to mate individuals in the mating pool. When two individuals
$\theta_{A}$ and $\theta_{B}$ reproduce, two newborn individuals $\theta_{1}$
and $\theta_{2}$ originate as random numbers belonging to the interval
${[\theta_{A}-c_{i}(\theta_{B}-\theta_{A}),\theta_{B}+c_{i}(\theta_{B}-\theta_{A})]}$,
where $c_{i}$ tunes the crossover, with ${i\in\\{1,2\\}}$, and we have assumed
${\theta_{B}\geq\theta_{A}}$. A similar reproduction occurs for two
individuals $\beta_{A}$ and $\beta_{B}$.

To explore a wider region of the parameter landscape, genetic mutations are
included in the workflow in the form of Gaussian noise with mean $\mu$ and
standard deviation $\sigma$, potentially affecting each newborn individual
Kramer (2017). Our GAs also include an elitism mechanism, i.e., the best
individual from the old population is carried over to the next one, replacing
the worst individual of the offspring. This pushes the algorithms to a faster
convergence toward the best solutions. To preserve the physical validity of
the final prediction for the maximum tilt angle $\theta_{m}$, a modulo-$\pi/2$
is performed after each operation on a $\theta$-individual. The maximum number
of generations $N_{\text{gen}}$ is used as the termination criterion.
Algorithm 1 presents the pseudo-code of the implemented GAs. Here, ${N=100}$,
$N_{\text{gen}}=50$, ${k=4}$, ${c_{1}=c_{2}=0.5}$, ${\mu=0}$, and
${\sigma=0.2}$. Blend crossover and mutation are non-deterministic operators
and occur with probability ${p_{c}=0.9}$ and ${p_{m}=0.01}$, respectively. The
desired convergence for $\theta_{m}$ and $\beta$ is typically achieved within
10 iterations. Adequate convergence is considered to be when the cost
functions are minimized with differences less than $10^{-15}$.

Algorithm 1 Pseudo-code of the implemented genetic algorithms

size of the population $pop\\_size$, tournament size $k$, crossover
probability $p_{c}$, $c$ for blend crossover, mutation probability $p_{m}$,
$\mu$ and $\sigma$ for Gaussian mutation, termination criterion $t$

the best solution $best$

$gen\leftarrow 0$

$pop\leftarrow$ generateRandomPopulation($pop\\_size$)

checkPhysicalConstraints($pop$)

evaluateFitness($pop$)

$best\leftarrow$ getBestIndividual($pop$)

while $gen<t$ do

$offspring\leftarrow$ executeTournament($pop,k$)

executeBlendCrossover(offspring, $p_{c}$, $c$)

checkPhysicalConstraints($offspring$)

executeGaussianMutation($offspring,p_{m},\mu,\sigma$)

checkPhysicalConstraints($offspring$)

evaluateFitness($offspring$)

$pop\leftarrow offspring$

$pop\leftarrow$ elitism($pop,best$)

$best\leftarrow$ getBestIndividual($pop$)

$gen\leftarrow gen+1$

end while

return $best$

Once the best estimates for $\theta_{m}$ and $\beta$ have been determined, a
similar evolutionary strategy is devised for retrieving the tilt distribution,
and then Eq. (3) can be used to find the twist distribution. The symmetry of
the director around the mid-plane of the cell restricts the optimization to
the first half of the cell: ${0<z<L/2}$. Half the cell thickness is divided
into small intervals—here, 50 intervals are used—and the GA is executed within
each slice. The set of solutions
$\\{\theta(0),\theta(z_{1}),\theta(z_{2}),...,\theta(L/2)\\}$ of Eq. (4) is
first determined by minimizing the following cost function within each
interval:

$\mathcal{L}_{\theta}=\left|\frac{z}{L}\int_{0}^{\theta_{m}}\frac{\sqrt{1+\kappa\sin^{2}\theta}}{g(\theta)}\text{d}\theta-\frac{1}{2}\int_{0}^{\theta(z)}\frac{\sqrt{1+\kappa\sin^{2}\theta}}{g(\theta)}\text{d}\theta\right|^{2}.$
(12)

At each propagation distance $z_{i}$, the corresponding solution
$\theta(z_{i})$ is then used to determine $\phi(z_{i})$ via Eq. (3). Since the
twist and tilt distributions are not expected to feature singular behaviours,
the solution found at a given position $z_{i}$ cannot be too different from
the solution associated with $z_{i-1}$ and $z_{i+1}$. Therefore, the initial
population of the GA performed at each position other than $z=0$—which is
known from boundary conditions—can be initialized from the solution found at
the previous position, perturbed with a uniform noise $\Delta$. Here, $\Delta$
is chosen between 0.05 and 0.1, depending on how close to $V_{T}$ the current
voltage is. This allows initializing the current GA very close to the actual
solution. Accordingly, fewer generations are needed to obtain an adequate
convergence, greatly reducing the computation time. The algorithms were
performed using MATLAB R2021B on a laptop with an 11th Gen Intel® CoreTM
i5-1145G7 CPU @ 2.60 GHz, 2611 MHz, 4 cores, and 8 logical processors. The
total calculation run time for a given set of initial conditions was about 2
min for the $\beta$–$\theta_{m}$ algorithm, and around 5 min for the
$\theta$–$\phi$ algorithm.

Figure 1: Integration parameter $\beta(\phi_{m},V)$. The dots correspond to
the numerically calculated values for maximum twist angles
${\phi_{m}=90^{\circ},67.5^{\circ},45^{\circ},22.5^{\circ},0^{\circ}}$. The
vertical gradient line is the analytical $\beta_{T}(\phi_{m})$ values from Eq.
(9), whereas the cascading solid-coloured lines are the lines of best fit
using Eq. (13).  Figure 2: Maximum tilt angle $\theta_{m}(\phi_{m},V)$. The
dots correspond to the numerically calculated values for maximum twist angles
${\phi_{m}=90^{\circ},67.5^{\circ},45^{\circ},22.5^{\circ},0^{\circ}}$. In the
bottom inset, the gradient line is the analytical $V_{T}(\phi_{m})$ values
from Eq. (8). The dashed lines in both insets connect the dots for visual
ease.

## Numerical Results

For our simulations, we use the material parameters for 6CHBT nematic liquid
crystals, which have elastic constants of $K_{1}=6.7$ pN, $K_{2}=3.4$ pN,
$K_{3}=10.6$ pN, and ${\Delta\epsilon=8}$ Nowinowski-Kruszelnicki  _et al._
(2012). The zero-twist threshold voltage is $V_{T0}=0.966~{}V$. Figures 1-2
report the numerically computed $\beta$ and $\theta_{m}$. The GA is run for
maximum twist angles of
${\phi_{m}=90^{\circ},~{}67.5^{\circ},~{}45^{\circ},22.5^{\circ}}$, and
$0^{\circ}$, at a range of voltages from 0.96 V to 7 V, with resulting cost
function $\mathcal{L}_{\beta}$ values of less than $10^{-32}$, and
$\mathcal{L}_{\theta_{m}}$ between $10^{-17}$ and $10^{-26}$. For
$\beta(\phi_{m},V)$, the first measure of whether the GA is working is whether
it can match the threshold $\beta_{T}(\phi_{m})$ values given by Eq. (9). As
shown in Fig. 1, the GA imitates the correct $\beta_{T}(\phi_{m})$ at each
maximum twist angle’s threshold voltage. The data points for each $\phi_{m}$
follow a smooth decreasing trend, asymptotically approaching zero for large
voltages. We propose a semi-analytical form to fit $\beta(\phi_{m},V)$:

$\begin{split}\beta(\phi_{m},V)=&\beta_{T}(\phi_{m})\bigg{(}1-\frac{2}{\pi}\times\\\
&\arctan\left[\sum_{i=1}^{4}s_{i}V^{i}\sqrt{V-V_{T}(\phi_{m})}\right]\bigg{)},\end{split}$
(13)

where $s_{i}$ are a set of fitting parameters for all $\phi_{m}$. Equation
(13) exhibits the expected behaviour:

1\. $\beta(\phi_{m},V)$ does not exist for $V<V_{T}(\phi_{m})$.

2\. When $V=V_{T}(\phi_{m})$, then $\beta(\phi_{m},V)=\beta_{T}(\phi_{m})$.

3\. $\beta(\phi_{m},V)\rightarrow 0$ as $V\rightarrow\infty$.

The computed data set for $\beta(45^{\circ},V)$ is used to obtain the fitting
parameters $s_{1}=-7.95043$, $s_{2}=16.5784$, $s_{3}=-10.5245$, and
$s_{4}=2.35869$, yielding a coefficient of determination
${R^{2}_{45}=99.99\%}$. With these $s_{i}$, the coefficients of determination
for the other data sets are ${R^{2}_{90}=99.98\%}$, ${R^{2}_{67.5}=99.97\%}$,
${R^{2}_{22.5}=99.98\%}$, and ${R^{2}_{0}=100\%}$.

Figure 3: Fits for $\theta_{m}(\phi_{m},V)$. With the ansatz of Eq. (14), fits
are produced using the ${\phi_{m}=0^{\circ},45^{\circ},90^{\circ}}$ data sets.
Each fit is then used to produce the $\theta_{m}(\phi_{m},V)$ curves for other
$\phi_{m}$ values, as shown in each subplot. Figure 4: GA results for
${\phi_{m}=45^{\circ}}$. a. Twist $\phi(z)$, b. tilt $\theta(z)$, and c. phase
retardation $\Gamma(z)$ distributions numerically calculated for a range of
voltages between 1.0 V and 4.2 V. The $\Gamma(z)$ are computed using
$L=35~{}\mu$m, $\Delta n=0.151$, and $\lambda=632$ nm.

For the maximum tilt angle $\theta_{m}(\phi_{m},V)$, we expect
${\theta_{m}=0}$ when ${V<V_{T}}$ for each $\phi_{m}$. The bottom inset of
Fig. 2 shows that our GA reproduces the threshold voltage values of Eq. (8).
The top inset of Fig. 2 is a zoom-in to show what appears to be a common
crossing point around 1.2 V with $\theta_{m}\sim 30.9^{\circ}$. These curves
can each be fit with the form,

$\theta_{m}(\phi_{m},V)=\arctan\left[\sum_{i=1}^{4}b_{i}V^{i}\sqrt{V-V_{T}(\phi_{m})}\right],$
(14)

where $b_{i}$ are a new set of fitting parameters. This ansatz obeys the
expected behaviour:

1\. When $V=V_{T}(\phi_{m})$, then $\theta_{m}(\phi_{m},V)=0$.

2\. $\theta_{m}(\phi_{m},V)\rightarrow 90^{\circ}$ as $V\rightarrow\infty$.

However, this form fails to reproduce the crossing point around 1.2 V for the
fitted data set, despite the corresponding $R^{2}$ values being 99.99%, as can
be seen in Fig. 3. This suggests that there is an extra $\phi_{m}$-dependence
not included in Eq. (14). A hint about this dependence can be found in Ref.
Preist  _et al._ (1989), which studies the twist and tilt distributions in the
high-voltage limit ${V\gg V_{T0}}$:

$\tan^{2}(\theta_{m})\approx\left(1+\tan^{2}(\phi_{m}/2)\right)\tan^{2}(\theta_{m}^{(0)}),$
(15)

where $\theta_{m}^{(0)}$ is the maximum tilt angle for a non-twisted cell at
the same voltage. However, an analytical expression for voltages below or near
the threshold was not found in the literature. If a semi-analytical form could
be found that reproduces the GA outputs for $\theta_{m}(\phi_{m},V)$, then
along with Eq. (13) for $\beta(\phi_{m},V)$ one would no longer need to run
the GA for every $(\phi_{m},V)$ configuration, dramatically reducing
computation time.

Figure 4a-b reports an example of the numerically obtained twist and tilt
distributions for ${\phi_{m}=45^{\circ}}$ with various voltage settings. The
twist distribution is no longer linear as we increase the applied field
strength (see Fig. 4a), as was first calculated by Deuling Deuling (1974).
Nevertheless, there appears to be a voltage above the threshold up to which
$\phi(z)$ is essentially still linear. As the field strength increases
further, the distributions feature an S-like shape, becoming sharper and more
step-like. As expected, the liquid crystals start tilting in the direction of
the applied field only above the threshold voltage (see Fig. 4b) Deuling
(1974). Above the threshold, the maximum tilt angle steadily approaches
$\theta_{m}=90^{\circ}$ with increasing voltage. For very intense fields, the
distribution flattens out in the middle. The phase retardation $\Gamma(z_{j})$
within the $j^{\text{th}}$ slice in Fig. 4c is calculated from the tilt
distribution using the trapezoidal rule in favour of left or right Riemann
sums for a better estimate,

$\Gamma(z_{j})=\frac{\pi\Delta
nd}{\lambda}[\cos^{2}\theta(z_{j+1})+\cos^{2}\theta(z_{j})].$ (16)

A cell thickness of $L=35~{}\mu$m is used, with ${N=100}$ slices for
${d=L/N}$, $\Delta n=0.151$, and $\lambda=632$ nm. The twist and tilt
distributions, along with the phase retardation, have also been extracted for
different maximum twist angles $\phi_{m}$ at a range of voltages (see Fig. 5).

Figure 5: GA results for different maximum twist angles. Twist, tilt, and
phase retardation distributions at a. $V=1.061$ V, b. $V=2.0$ V, c. $V=4.0$ V.
Figure 6: Validation of the twist scaling approximation. Comparison of output
states generated via the approximated twist and tilt distributions (colored
lines), based on ${\phi(45^{\circ},V,z)}$ and $\theta(45^{\circ},V,z)$, and
the numerically calculated distributions for each $\phi_{m}$ (colored dots),
at applied voltages of $V=3.00$, 8.00, and 11.31$~{}V$. A horizontally
polarized input was considered.

Using the numerically calculated $\phi(z)$, $\theta(z)$, and Eq. (16), we can
derive the total Jones matrix of a TNLC cell for any given $\phi_{m}$ and
applied voltage $V$. In the case with no external field, the device is
modelled as a stack of $N$ linearly twisted cells of thickness $d$, with
constant phase retardation. The Jones matrix
$\mathbf{T}_{\phi_{f}}(\phi_{m},\Gamma)$ for a given
${\phi_{m}=\phi_{b}-\phi_{f}}$ and phase retardation $\Gamma$ is given by
Yariv and Yeh (1984),

$\displaystyle\mathbf{T}_{\phi_{f}}(\phi_{m},\Gamma)$ $\displaystyle=$
$\displaystyle\mathbf{R}(-\phi_{b})\mathbf{M}_{0}(\phi_{m},\Gamma)\mathbf{R}(\phi_{f}),$
(17)

where,

$\mathbf{M}_{0}(\phi_{m},\Gamma)=\begin{bmatrix}\cos X-\frac{i\Gamma}{2X}\sin
X&\frac{\phi_{m}}{X}\sin X\\\ -\frac{\phi_{m}}{X}\sin X&\cos
X+\frac{i\Gamma}{2X}\sin X\end{bmatrix},$ (18)

$X=\sqrt{\phi_{m}^{2}+(\Gamma/2)^{2}}$, $\phi_{f}=\phi(0)$, and
$\phi_{b}=\phi(L)$ are the front and back alignment angles, respectively, and
$\mathbf{R}(\cdot)$ is the rotation matrix,

$\mathbf{R}(\cdot)=\begin{bmatrix}\cos(\cdot)&\sin(\cdot)\\\
-\sin(\cdot)&\cos(\cdot)\end{bmatrix}.$ (19)

When a voltage $V$ is applied across the cell, the total Jones matrix can be
then approximated as,

$\displaystyle\mathbf{J}_{\phi_{f}}(\phi_{m},V)$ $\displaystyle=$
$\displaystyle\prod_{j=0}^{N}\mathbf{T}_{\phi(z_{j})}\left(\phi_{m}(z_{j}),\Gamma(z_{j})\right)$
$\displaystyle=$
$\displaystyle\mathbf{R}(-\phi_{b})\left[\prod_{j=0}^{N}\mathbf{M}_{0}(\phi_{m}(z_{j}),\Gamma(z_{j}))\right]\mathbf{R}(\phi_{f}),$

where $z_{j}=jd$, $\phi_{m}(z_{j})=\phi(z_{j+1})-\phi(z_{j})$.

To further save on computation time, we make the additional approximation that
the twist distributions for different $\phi_{m}$ are scaled versions of each
other, e.g.,
$\phi(\phi_{m},V,z)=(\phi_{m}/45^{\circ})\times\phi(45^{\circ},V,z)$. From the
various twist distributions shown in the left column of Fig. 5, this appears
to be a reasonable assumption. Consequently, the tilt and phase retardation
distributions are assumed to be the same for all $\phi_{m}$ at a given
voltage, i.e., ${\theta(\phi_{m},V,z)=\theta(45^{\circ},V,z)}$, and
${\Gamma(\phi_{m},V,Z)=\Gamma(45^{\circ},V,z)}$. This is particularly the case
for high voltages, with deviations expected at lower voltages (see Fig. 5).
Figure 6 compares the output polarizations plotted on the Poincaré Sphere at
different voltages for the full numerically calculated distributions for
${\phi_{m}=-90^{\circ}}$, $-67.5^{\circ}$, $-45^{\circ}$, $-22.5^{\circ}$,
$0^{\circ}$, $22.5^{\circ}$, $45^{\circ}$, $67.5^{\circ}$, $90^{\circ}$, with
the approximated distributions based on ${\phi_{m}=45^{\circ}}$. The
simulations are obtained from Eq. (Numerical Results), assuming a horizontally
polarized input state. The differences are minimal, with average state
overlaps ${(1+\mathbf{S}_{\phi_{m}}\cdot\tilde{\mathbf{S}}_{\phi_{m}})/2}$ of
over $99\%$, where $\mathbf{S}_{\phi_{m}}$ and $\tilde{\mathbf{S}}_{\phi_{m}}$
are the output Stokes vector resulting from the numerically calculated and
approximated distributions, respectively, for a given $\phi_{m}$.

Figure 7: Fabricated sample. Image of DP(0,1/2) between crossed polarizers
under a microscope, illuminated with white light. The topological pattern on
each glass plate is also shown. The ${q=1/2}$ pattern is discretized into 16
slices to explore maximum twist angles $\phi_{m}$ from $-90^{\circ}$ to
$90^{\circ}$. Figure 8: Experimental results. Comparison between numerical
predictions and experimentally reconstructed Stokes parameters for a.
DP(0,1/2), at $V_{pp}=6.00$, 9.00, and 12.00 $V$, and b. DP(0,1/2), at
$V_{pp}=3.40$, 5.60, and 8.30 $V$.

## Experimental Results

The predictions of our numerical routine are tested with TNLC plates with a
spatially varying maximum twist angle (see Ref. Sit  _et al._ (2023) for
details). In the following, we briefly review the fabrication technique of
these devices. Two glass plates with a thin layer of conductive indium tin
oxide (ITO) are coated with an azobenzene-based dye. The dye molecules are
photoaligned when exposed to linearly polarized light at a wavelength within
the peak of the absorption spectrum. In our realization, we separately
patterned the glass plates such that the front plate features a uniform
alignment, and the back plate is patterned with a ${q=1/2}$ topology which has
been discretized into 16 slices. Finally, the two plates are sealed and the
nematic liquid crystals (6CHBT) penetrate into the sample by capillarity,
exhibiting a full range of maximum twist angles $\phi_{m}$ from $-90^{\circ}$
to $90^{\circ}$. Figure 7 shows the fabricated sample between crossed
polarizers under a microscope, illuminated with white light. We refer to these
inhomogeneous non-symmetrically patterned devices as dual-plates (DPs), as
they exhibit a different behaviour depending on the plate orientation. For the
DP used here, the configuration in which light passes through the uniform
pattern first and exits through the $q=1/2$ pattern is denoted as DP(0,1/2),
and vice versa as DP(1/2,0).

The action of each configuration on polarized light is characterized via
polarization tomography, to reconstruct the output polarization distribution
when different voltages are applied across the cell. Figure 8a-b compares the
experimentally reconstructed Stokes parameters with the predicted outputs
obtained from the approximated total Jones matrices of DP(0,1/2) and
DP(1/2,0), respectively, at different voltages. In our experiment, a
sinusoidal waveform with 4 kHz frequency was used. An excellent agreement is
observed in all realizations, with average overlaps of $93\pm 2\%$, $98.0\pm
0.3\%$, and $92.7\pm 0.4\%$ for DP(0,1/2), at voltages ${V_{pp}=6.00}$, $9.00$
and $12.00$ V, respectively (see Fig. 8a), where $V_{pp}$ is the peak-to-peak
voltage, and the average is computed over the outputs within each of the 16
slices. For DP(1/2,0), we obtain $99.7\pm 0.1\%$, $95\pm 1\%$, and $99.2\pm
0.3\%$, at ${V_{pp}=3.40}$, $5.60$ and $8.20$ V, respectively (see Fig. 8b).
Deviations from numerical predictions are mainly ascribed to fabrication
defects and environmental temperature fluctuations, which can slightly change
the liquid-crystal intrinsic birefringence and elastic constants. These
results certify that our numerical routines are suitable for predicting the
optical action of individual devices with high accuracy.

## Conclusion

We have demonstrated a novel robust approach to the determination of twist and
tilt distributions of liquid-crystal cells in the presence of an external
field. Our method directly tackles the integral form of the Euler-Lagrange
equations, thereby avoiding the necessity of trial solutions. The complexity
of the equation system is fragmented within subsequent genetic routines, each
of which uses the outputs of the previous one to converge to the optimal
solutions. Our method has been validated both numerically and experimentally
on two configurations of spatially varying TNLC cells, where the optimization
runs over multiple transverse positions. This scheme can provide a useful tool
for the experimental characterization of the next generations of dual-devices,
such as dual-lenses and gratings. At the same time, it will be interesting to
explore machine-learning approaches to extract the liquid-crystal director
distributions in real-time Sigaki  _et al._ (2020); Piven  _et al._ (2024).

## Acknowledgements

A.S. acknowledges the financial support of the Vanier graduate scholarship of
the NSERC. This work was supported by the Ontario’s Early Researcher Award
(ERA), Canada Research Chairs (CRC) and Natural Sciences and Engineering
Research Council of Canada (NSERC).

## References

  * Leslie (1970) F. M. Leslie, Molecular Crystals and Liquid Crystals 12, 57 (1970).
  * Deuling (1974) H. J. Deuling, Molecular Crystals and Liquid Crystals 27, 81 (1974).
  * De Gennes and Prost (1993) P. G. De Gennes and J. Prost,  _The Physics of Liquid Crystals_ , 2nd ed. (Oxford: Clarendon Press, 1993).
  * Preist  _et al._ (1989) T. W. Preist, K. R. Welford, and J. R. Sambles, Liquid Crystals 4, 103 (1989).
  * Yang and Wu (2014) D.-K. Yang and S.-T. Wu,  _Fundamentals of liquid crystal devices_ (John Wiley & Sons, 2014).
  * Moser  _et al._ (2019) S. Moser, M. Ritsch-Marte, and G. Thalhammer, Optics Express 27, 25046 (2019).
  * Wang  _et al._ (2004) B. Wang, X. Wang, and P. J. Bos, Journal of the Optical Society of America A 21, 1066 (2004).
  * Wang  _et al._ (2005) X. Wang, B. Wang, P. J. Bos, J. E. Anderson, J. J. Pouch, and F. A. Miranda, Journal of the Optical Society of America A 22, 346 (2005).
  * Peng  _et al._ (2016) F. Peng, Y. Huang, F. Gou, M. Hu, J. Li, Z. An, and S.-T. Wu, Optical Materials Express 6, 717 (2016).
  * Lee  _et al._ (2017) Y.-H. Lee, D. Franklin, F. Gou, G. Liu, F. Peng, D. Chanda, and S.-T. Wu, Scientific Reports 7, 16260 (2017).
  * Guirado  _et al._ (2022) R. Guirado, G. Perez-Palomino, M. Ferreras, E. Carrasco, and M. Caño-García, IEEE Transactions on Antennas and Propagation 70, 11847 (2022).
  * Gorkunov  _et al._ (2022) M. V. Gorkunov, A. V. Mamonova, I. V. Kasyanova, A. A. Ezhov, V. V. Artemov, I. V. Simdyankin, and A. R. Geivandov, Nanophotonics 11, 3901 (2022).
  * Sova  _et al._ (2023) O. Sova, T. Sluckin, S. Kaur, H. Gleeson, and V. Reshetnyak, Optical Materials Express 13, 1640 (2023).
  * Yang and Collings (2020) S. Yang and P. J. Collings, Crystals 10, 1041 (2020).
  * Holland (1992) J. H. Holland,  _Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control and Artificial Intelligence_ (Cambridge: MIT Press, 1992).
  * Yao (1993) X. Yao, Microprocessing and Microprogramming 38, 707 (1993).
  * Sit  _et al._ (2023) A. Sit, F. Di Colandrea, A. D’Errico, and E. Karimi, arXiv:2311.04773 (2023).
  * Yamauchi (2005) M. Yamauchi, Applied optics 44, 4484 (2005).
  * Deuling (1972) H. J. Deuling, Molecular Crystals and Liquid Crystals 19, 123 (1972).
  * Iri  _et al._ (1970) M. Iri, S. Moriguti, and Y. Takasawa, RIMS Kokyuroku Kyoto University 91, 82 (1970).
  * Davis and Rabinowitz (2014) P. Davis and P. Rabinowitz,  _Methods of Numerical Integration_ (San Diego: Academic Press, 2014).
  * Goldberg and Deb (1991) D. E. Goldberg and K. A. Deb, Foundations of Genetic Algorithms 1, 69 (1991).
  * Eshelman and Schaffer (1993) L. J. Eshelman and J. D. Schaffer, Foundations of Genetic Algorithms 2, 187 (1993).
  * Kramer (2017) O. Kramer,  _Genetic Algorithm Essentials_ (Springer, 2017).
  * Nowinowski-Kruszelnicki  _et al._ (2012) E. Nowinowski-Kruszelnicki, J. Kedzierski, Z. Raszewski, L. Jaroszewicz, M. A. Kojdecki, W. Piecek, P. Perkowski, M. Olifierczuk, E. Miszczyk, K. Ogrodnik, and P. Morawiak, Opto-electronics Review 20, 255 (2012).
  * Yariv and Yeh (1984) A. Yariv and P. Yeh,  _Optical Waves in Crystals: Propagation and Control of Laser Radiation_ (New York: John Wiley I& Sons, 1984).
  * Sigaki  _et al._ (2020) H. Y. Sigaki, E. K. Lenzi, R. S. Zola, M. Perc, and H. V. Ribeiro, Scientific Reports 10, 7664 (2020).
  * Piven  _et al._ (2024) A. Piven, D. D. Darmoroz, E. Skorb, and T. Orlova, Soft Matter (2024).