# The effects of cavitation position on the velocity of a laser-induced
microjet extracted using explainable artificial intelligence

Daichi Igarashi  Jingzu Yee  Yuto Yokoyama  Hiroaki Kusuno  Yoshiyuki Tagawa
Authors to whom correspondence should be addressed<EMAIL_ADDRESS>Also at the Institute of Global Innovation Research, Tokyo University of
Agriculture and Technology, Koganei, Tokyo 184-8588, Japan  Department of
Mechanical Systems Engineering, Tokyo University of Agriculture and
Technology, Koganei, Tokyo 184-8588, Japan

###### Abstract

The control of the velocity of a high-speed laser-induced microjet is crucial
in applications such as needle-free injection. Previous studies have indicated
that the jet velocity is heavily influenced by the volumes of secondary
cavitation bubbles generated through laser absorption. However, there has been
a lack of investigation of the relationship between the positions of
cavitation bubbles and the jet velocity. In this study, we investigate the
effects of cavitation bubbles on the jet velocity of laser-induced microjets
extracted using explainable artificial intelligence (XAI). An XAI is used to
classify the jet velocity from images of cavitation bubbles and to extract
features from the images through visualization of the classification process.
For this purpose, we run 1000 experiments and collect the corresponding
images. The XAI model, which is a feedforward neural network (FNN), is trained
to classify the jet velocity from the images of cavitation bubbles. After
achieving a high classification accuracy, we analyze the classification
process of the FNN. The predictions of the FNN, when considering the
cavitation positions, show a higher correlation with the jet velocity than the
results considering only cavitation volumes. Further investigation suggested
that cavitation that occurs closer to the laser focus position has a higher
acceleration effect. These results suggest that the velocity of a high-speed
microjet is also affected by the cavitation position.

††preprint: AIP/123-QED

##  I Introduction

Currently, injection needles are one of the most common devices used for drug
delivery. However, the use of needles for drug delivery can cause problems
such as trypanophobia Nir  _et al._ (2003); McLenon and Rogers (2019); Deacon
and Abramowitz (2006) and the spread of infectious diseases.Kermode (2004);
Weiss and McMichael (2004); res (2003) Therefore, a needle-free injector,
which injects liquid as a high-speed microjet into the body, has been
developed in recent years. Mitragotri (2005); Stachowiak  _et al._ (2009);
Schoppink and Rivas (2022); Daly  _et al._ (2020) Commercialized needle-free
injectors, which commonly use a diffused jet, are reported to cause
significant pain during injection.Harding  _et al._ (2002); Giudice and
Campbell (2006); Mitragotri (2006); Arora  _et al._ (2007) To solve this
problem, a laser-induced microjet with a focused tip is expected to be useful
because of the lower stress induced at the injection site and the higher
efficiency of liquid transportation.Tagawa  _et al._ (2012, 2013); Miyazaki
_et al._ (2021); Kiyama  _et al._ (2019) To develop an efficacious needle-free
injector using a focused microjet, it is necessary to control the jet
velocity.

Many studies have been carried out with the aim of understanding the mechanism
that affects jet velocity.Franco-Gómez  _et al._ (2021); Van Hoeve  _et al._
(2010); Scroggs and Settles (1996); Choo and Kang (2007); Thoroddsen  _et al._
(2009); Kamamoto, Onuki, and Tagawa (2021) Several of these studies have
reported that the total volumes of secondary cavitation bubbles in the liquid
affects the jet velocity, and have proposed a linear model to explain the
relationship.Kiyama  _et al._ (2016); Hayasaka, Kiyama, and Tagawa (2017);
Daily  _et al._ (2014); Ishikawa, Nishida, and Tagawa (2022) According to a
study by Kiyama et al.,Kiyama  _et al._ (2016) this process can be explained
as follows. When a laser is absorbed by a liquid in a microtube, compression
waves are generated, which accelerate the jet. These compression waves are
reflected at the gas-liquid interface as expansion waves, which decelerate the
jet. When cavitation bubbles occur at the positions where the pressure in the
liquid falls below the saturation vapor pressure, the cavitation bubbles
inhibit the propagation of the expansion waves. Therefore, the greater the
amount of cavitation, the more easily the propagation of the expansion wave is
inhibited, and the faster the jet velocity. However, plots in the same paper
of the jet velocity versus the total volumes of cavitation bubbles still
showed a non-negligible amount of variance. After reviewing these previous
experimental data, we have found that multiple cavitation bubbles appear at
different positions. Therefore, we argue that the cavitation bubbles’
positions, as well as the total volumes of the cavitation bubbles, may affect
the jet velocity.

To investigate the influence of both the volumes and positions of cavitation
bubbles on the jet velocity, we propose the use of artificial intelligence
(AI), which has recently been the focus of a great deal of attention from the
fluid mechanics community,Brunton, Noack, and Koumoutsakos (2020); Colvert,
Alsalman, and Kanso (2018); Ling, Kurzawski, and Templeton (2016); Kutz (2017)
even in the field of cavitation bubbles.Folden and Aschmoneit (2023) Feature
extraction has been studied by fluid researchers in various fields, using
methods such as dynamic mode decomposition (DMD),Williams, Kevrekidis, and
Rowley (2015); Jovanović, Schmid, and Nichols (2014) principal component
analysis (PCA),Abdi and Williams (2010) and proper orthogonal decomposition
(POD).Kerschen  _et al._ (2005) In addition, many fluid studies have been
performed using convolutional neural networks (CNNs) for problems dealing with
images.Vennemann and Rösgen (2020); Sekar  _et al._ (2019); Liu  _et al._
(2020); Cai  _et al._ (2019); Jin  _et al._ (2018) However, for most AI
approaches, the underlying reasoning that leads to a specific decision is
unknown or not properly understood. In contrast, explainable AI (XAI),Arrieta
_et al._ (2020); Adadi and Berrada (2018); Gunning (2017); Gunning  _et al._
(2019); Yee  _et al._ (2023) which does not have the problems of
explainability and interpretability, represents a powerful tool that can aid
in the understanding of physical phenomena. Therefore, the objective of this
study is to clarify the effects of bubble volume and position on jet velocity
through feature extraction using a large amount of experimental data and XAI.
In a study by Yee et al.,Yee, Yamanaka, and Tagawa (2022) weight visualization
with XAI was performed to extract the features that differentiate splashing
and non-splashing drops when a droplet impacts a solid surface. We adopt the
same method here to extract the positions of the cavitation bubbles that
influence the jet velocity.

The structure of this paper is as follows. In Sec. II, the experimental setup
and the collection of experimental data are described. We also describe the
relationship between the amount of cavitation bubbles and the jet velocity
here. Next, in Sec. III, the XAI implementation and the methods of feature
extraction are explained. We also show results from the training and testing
of the XAI implementation. In Sec. IV, we explain the prediction and feature
extraction results of the FNN and discuss the results. We analyze the effects
of the positions of the cavitation bubbles on the jet velocity. Finally, the
conclusions and outlook of this study are presented in Sec. V.

##  II Experiments

This section describes the methodology of the experiments and the relation
between the amount of secondary cavitation bubbles and the jet velocity.

###  II.1 Experimental setup

The experimental setup is shown in Fig. 1. It is based on a study by Tagawa et
al. Tagawa  _et al._ (2012) A green pulsed laser generated by a laser
generator (Nd: YAG laser Nano SPVI, Litron Laser Co., Ltd., wavelength: 532
nm) was passed through a mirror, a half mirror (Dielectric plate half-mirror,
OptoSigma Co., transmission: 50%), and an objective lens (MPLn10x, Olympus Co.
Ltd., Japan, magnification: 10x, N.A. value: 0.25). We used magenta ink
(GI-30M, Canon), which has high green laser absorbance, as the working liquid.
The laser energy was measured with an energy meter (ES111C, Thorlabs, energy
range: 10 µJ–150 mJ). The image sequences of the laser-induced microjet were
filmed with a high-speed camera (FASTCAM SA-X, Photron Co., Ltd.). The
temporal resolution of the camera was 100 000 f.p.s. The lens was attached
with a high-pass filter (SCF-50S-56O, Sigma Koki Co., transmission limit
wavelength: 560 nm) to shade the 532 nm laser light. A light source
(SLG-150V-CW-MN, REVOX Inc.) was used as the backlight for the high-speed
camera. The laser generator and the high-speed camera were synchronized using
a delay generator (Model 575 Pulse/Delay generator, BNC Co.). In this study,
we acquired 1000 image sequences of the laser-induced microjet.

According to the study by Tagawa et al.,Tagawa  _et al._ (2012) the parameters
that affect the jet velocity are the laser energy $E$ [µJ], the inner diameter
of the microtube $d$ [mm], the contact angle between the meniscus interface
and the inner wall $\theta$ [∘], and the distance from the meniscus interface
to the laser irradiation position $H$ [mm], as shown in Fig. 1. In this study,
the ranges of these parameters are $E=$ 500–600 µJ, $d=$ 0.5 mm, $\theta=$
40.3–59.4∘, and $H=$ 1.84–2.09 mm, respectively. In addition, the same tube
was used throughout the experiments.

Figure 1:  A schematic diagram of the experimental setup for the laser-induced
microjet.

###  II.2 Experimental observations

Fig. 2 shows an example of the image sequence taken from the experiments
described above. We set the time of the frame immediately before the laser
illumination as $t=0$ µs. A laser-induced microjet was generated after the
laser was focused on the liquid ($t=$ 20–300 µs). During this process, we
observed an expanding bubble at the laser focusing point ($t=$ 10–300 µs).
This bubble, known as a vapor bubble, was created when the liquid was
vaporized by the heat of the laser. Akhatov  _et al._ (2001) After a certain
period of time, the vapor bubble began to shrink ($t\approx 200$ µs) and
finally collapsed ($t\approx 300$ µs).

In addition to the vapor bubble, we also observed small secondary bubbles in
the images at $t=$ 10–40 µs in Fig. 2. As shown in the figure, these secondary
bubbles generated at $t\approx$ 10 µs began to shrink at $t\approx$ 30 µs and
eventually disappeared at $t\approx$ 50 µs. Such bubbles, known as cavitation
bubbles, are caused by a rapid phase change from liquid to gas due to the
decrease in the pressure of the liquid.Caupin and Herbert (2006) Previous
studies have indicated that the jet velocity is related to these cavitation
bubbles. Kiyama  _et al._ (2016); Hayasaka, Kiyama, and Tagawa (2017); Daily
_et al._ (2014); Ishikawa, Nishida, and Tagawa (2022) Therefore, we
investigated the images in the frame immediately before the laser illumination
($t=0$ µs), the frame immediately after the laser illumination ($t=10$ µs),
and the succeeding frame ($t=20$ µs). At $t=0$ µs, there were no cavitation
bubbles. In the frames at $t=10$ µs, cavitation bubbles started to form. In
the frame at $t=20$ µs, the cavitation bubbles reached their maximum sizes in
most cases, as shown in Fig. 2.

Figure 2:  An image sequence of the laser-induced microjet is obtained from
the experiment. The vapor bubble is observed at $t=$ 10–250 µs, the cavitation
bubbles are observed at $t=$ 10–40 µs, and the microjet is observed at $t=$
20–300 µs.

###  II.3 The effects of the laser energy and amount of cavitation on the jet
velocity

The relationship between the jet velocity $U_{\mathrm{jet}}$ and the laser
energy $E$ for the 1000 data is shown in Fig. 3(a). As shown in the figure,
there is significant variation in the jet velocity for a fixed laser energy.
By comparing images with the same $E$ but different $U_{\textrm{jet}}$, we
found some images with large amounts of cavitation bubbles, some with small
amounts of cavitation bubbles, and others with no cavitation bubbles. Kiyama
et al.Kiyama  _et al._ (2016) pointed out that the jet velocity increases with
the total volumes of cavitation bubbles in the liquid. Therefore, we
investigated the relationship between the jet velocity
$U_{\mathrm{\mathrm{jet}}}$ and the total projected area of cavitation bubbles
$S_{\mathrm{\mathrm{cav}}}$ of the experimental data when the cavitation
bubbles reached their maximum sizes, i.e., $t=20$ µs. The results are shown in
Fig. 3(b). We calculated the cavitation area $S_{\mathrm{\mathrm{cav}}}$ from
the number of pixels occupied by cavitation bubbles in the binarized input
images. It was observed that the two roughly correlate with each other, with
the correlation coefficient of the jet velocity and the total area of
cavitation bubbles being 0.817. However, there is still a non-negligible
variance. A possible contributing factor to this issue could be the lack of
consideration of the positions of the cavitation bubbles.

Figure 3:  (a) The relationship between the jet velocity
$U_{\mathrm{\mathrm{jet}}}$ and laser energy $E$ in the experiment, where the
blue and red areas represent the data of the fastest 25% and the slowest 25%
of the jet velocity, respectively. (b) The relationship between the jet
velocity $U_{\mathrm{\mathrm{jet}}}$ and the total projected area of the
cavitation bubbles $S_{\mathrm{\mathrm{cav}}}$ at $t=20$ µs.

##  III Explainable artificial intelligence (XAI): feedforward neural network
(FNN)

To investigate the effects of the cavitation position on the jet velocity, we
used explainable artificial intelligence (XAI) to classify the jet velocity
from images with cavitation bubbles and extract the features of the images
through visualization of the classification process. In this section, we
describe the implementation of the XAI and the visualization of the
classification process.

Figure 4:  The architecture of the feedforward neural network for classifying
the “Fast” and “Slow” jet velocity classes in this study.

###  III.1 Implementation of the FNN

For the XAI, we used the architecture of a feedforward neural
network(FNN),Rosenblatt (1958) as shown in Fig. 4. An FNN is a simple neural
network in which the decision-making process can be explained through
visualization of the classification process.Yee, Yamanaka, and Tagawa (2022)
According to a study by Hornik and Stinchcombe,Hornik, Stinchcombe, and White
(1989) an FNN with three or more layers can approximate a continuous function
with arbitrary precision (universal approximation theorem). Therefore, we used
a three-layered FNN. In addition, we adopted the method of Yee et al.Yee,
Yamanaka, and Tagawa (2022) to address the problem as a binary classification
problem to predict whether the jet velocity was “Fast” or “Slow”, and
visualized the image features of each classification.

The captured images were preprocessed before being input into the FNN. First,
we cropped the image of all experimental data to focus on the phenomena in the
microtube. The cropping range for each image was 391$\times$111 pix in the
area enclosed by the blue, green, and red dashed boxes in Fig. 2. Next, we
normalized the brightness values of the images. To feed an image into the FNN,
each two-dimensional matrix of the input image was flattened into a one-
dimensional vector. The number of neurons in the input layer equals the total
number of pixels in an input image. The input vectors of the hidden layer are
defined as $\mathbf{q_{\mathrm{1}}}$:

$\mathbf{q_{\mathrm{1}}=W_{\mathrm{1}}x+b_{\mathrm{1}}}.$ (1)

Here, $\mathbf{W_{\mathrm{1}}}$ and $\mathbf{b_{\mathrm{1}}}$, respectively,
represent the weight matrix and bias vector that connect the input and hidden
layers. In addition, $\mathbf{x}$ denotes the flattened one-dimensional vector
of an input image.

We set the number of neurons in the hidden layer to five, which is the optimal
number obtained through parameter tuning. In the hidden layer, we used the
leaky rectified linear unit (Leaky ReLU) function Maas, Hannun, and Ng (2013)
as the activation function. By representing the output vectors of the hidden
layer with $\mathbf{s_{\mathrm{1}}}$, the equation in the hidden layer could
be expressed as

$\mathbf{s_{\mathrm{1}}}=\mathrm{Leaky\,ReLU}(\mathbf{q_{\mathrm{1}}})=\begin{cases}\mathbf{q_{\mathrm{1}}}&\textrm{if}\;\mathbf{q_{\mathrm{1}}}\geq
0\\\ \alpha\mathbf{q_{\mathrm{1}}}&\textrm{if}\;\mathbf{q_{\mathrm{1}}}<0,\\\
\end{cases}$ (2)

where $\alpha$ is set to 0.01, which is the optimal value as determined by
hyperparameter tuning.

By defining the input vectors of the output layer as
$\mathbf{q_{\mathrm{out}}}$, the equation that connects the hidden and output
layers can be expressed as

$\mathbf{q_{\mathrm{out}}=W_{\mathrm{out}}s_{\mathrm{1}}+b_{\mathrm{out}}}.$
(3)

Here, $\mathbf{W_{\mathrm{out}}}$ and $\mathbf{b_{\mathrm{out}}}$ represent
the weight matrix and bias vector that connect the hidden and output layers.

For the classification problem, out of the 1000 data, we labeled the 250 with
the highest velocities ($U_{\rm{jet}}>23.8$ m/s) as “Fast” and the 250 with
the lowest velocities ($U_{\rm{jet}}<17.6$ m/s) as “Slow”, as shown in Fig.
3(a). Note that the blue and red areas in the figure represent the data of the
fastest 25% and the slowest 25% of the jet velocity, respectively. Thus, the
number of neurons in the output layer was two, i.e., the number of classes.

In the output layer, we used the sigmoid function as the activation function.
By defining the output vectors of the output layer as
$\mathbf{y_{\mathrm{pred}}}$, the equation in the output layer could be
expressed as

$\mathbf{y}_{\mathrm{pred}}=\mathrm{\sigma}(\mathbf{q_{\mathrm{out}}})=\frac{1}{1+\exp(-\mathbf{q_{\mathrm{out}}})}.$
(4)

Here, we define the two elements of the $\mathbf{y}_{\mathrm{pred}}$ vector,
which give the probabilities of an image being predicted to be “Fast” and
“Slow”, as $y_{\mathrm{pred1}}$ and $y_{\mathrm{pred2}}$, respectively. A
prediction of “Fast” is made when $y_{\mathrm{pred1}}>y_{\mathrm{pred2}}$ and
a prediction of “Slow” is made when $y_{\mathrm{pred1}}<y_{\mathrm{pred2}}$.

###  III.2 Training and testing

Figure 5:  An illustration of the cross-validation: the combinations of
training, validation, and test data used in this study (where comb. means
combination).

We evaluated the generalization performance of the FNN using cross-
validation.Géron (2019) As explained in Sec. III.1, we used 250 “Fast” data
and 250 “Slow” data as the training data for the classification problem. Note
that the remaining 500 data were used as the interpolation test data in Sec.
IV. As shown in Fig. 5, we divided the 500 training data into five parts of
100. For a given data combination, one of the five parts was used for testing,
while the remaining four parts were used for training and validation.
Therefore, there were five data combinations in total. We set the number of
“Fast” and “Slow” data to be equal in both the training and testing datasets
for each combination.

We trained the neural network described in Sec. III.1 with momentum stochastic
gradient descent (SGD)Li  _et al._ (2014) as the optimizer. We performed mini-
batch training with the number of data in each batch set to 32. The validation
data used was a randomly selected 20% of all data excluding test data. The
evaluation of the validation losses for each epoch is shown in Fig. 6. Here,
the number of epochs represents how many times an entire training dataset is
fed to the FNN. The figure shows that the error decreases with each epoch for
all combinations. Therefore, the training validated the generalization
performance. Note that in order to avoid overfitting the training data, we
used early stopping,Prechelt (2002) a method that terminates the training
process when the error increases over that of the previous epoch.

After training, we investigated the test accuracy of the classification for
the jet velocity from the trained FNN using images at different times. Tables
1, 2, and 3 show the accuracy when using images at $t=0$ µs (the blue dashed
box in Fig. 2), $t=10$ µs (the green dashed box in Fig. 2), and $t=20$ µs (the
red dashed box in Fig. 2), respectively. Each table shows the test accuracy of
the “Fast”, “Slow”, and “Total” predictions of each combination and the
average accuracy of all combinations.

Figure 6:  An evaluation of the validation losses for every epoch of the training using each combination. Note that Combination 5 continued to train until about 2500 epochs. Table 1:  The test accuracy of the “Fast” and “Slow” predictions in each combination for $t=0$ µs. Combination | Fast | Slow | Total
---|---|---|---
1 | 50/50 (100%) | 0/50 (0%) | 50/100 (50%)
2 | 50/50 (100%) | 0/50 (0%) | 50/100 (50%)
3 | 50/50 (100%) | 0/50 (0%) | 50/100 (50%)
4 | 47/50 (94%) | 11/50 (22%) | 58/100 (58%)
5 | 47/50 (94%) | 21/50 (42%) | 68/100 (68%)
Average | 97.6% | 12.8% | 55.2%

At $t=0$ µs, Table 1 shows that most combinations have a total accuracy of
50%. This is the minimum prediction accuracy because there are only two
classification classes and the number of “Fast” and “Slow” test data is equal.
Therefore, this accuracy can be obtained just by classifying all the test data
as one of the two classes. From above, with regard to the accuracy using the
images at $t=0$ µs, the FNN predicted almost all data as “Fast”, which means
the jet velocity was not accurately predicted. This is because there is not
much difference between the images at $t=0$ µs.

Next, we compared the test results at $t=10$ µs (Table 2) and $t=20$ µs (Table
3). As shown in Table 2, at $t=10$ µs, the accuracy of the “Fast” and “Slow”
cases is higher than 70.0% for all five combinations. Although it is higher
than that at $t=0$ µs, there are still inaccurate predictions in some cases.
However, as shown in Table 3, these wrongly predicted data are correctly
predicted when images at $t=20$ µs are used for training. Therefore, we
inspected the “Fast” data that were wrongly predicted using images at $t=10$
µs but correctly predicted using images at $t=20$ µs. The investigation shows
that for most of these data, cavitation bubbles appeared in the images at
$t=20$ µs but did not in the images at $t=10$ µs. This suggests that the
presence of cavitation bubbles affects the test accuracy. Finally, the
averaged total accuracy was 99.2%, indicating that the FNN extracted a
straightforward feature that affects the jet velocity. Note that the increase
in the accuracy of the “Fast” data at $t=20$ µs improves the accuracy of the
“Slow” data because of the nature of a binary classification problem. In this
study, to understand the effects of the cavitation bubbles, we investigate the
images at $t=20$ µs when the volumes of the cavitation bubbles and the
averaged total accuracy are maximal among $t=0$ µs to $t=150$ µs.

Table 2:  The test accuracy of the “Fast” and “Slow” predictions in each combination for $t=10$ µs. Combination | Fast | Slow | Total
---|---|---|---
1 | 38/50 (76%) | 47/50 (94%) | 85/100 (85%)
2 | 49/50 (98%) | 46/50 (92%) | 95/100 (95%)
3 | 50/50 (100%) | 41/50 (82%) | 91/100 (91%)
4 | 50/50 (100%) | 42/50 (84%) | 92/100 (92%)
5 | 50/50 (100%) | 35/50 (70%) | 85/100 (85%)
Average | 94.8% | 84.4% | 89.6%
Table 3:  The test accuracy of the “Fast” and “Slow” predictions in each combination for $t=20$ µs. Combination | Fast | Slow | Total
---|---|---|---
1 | 49/50 (98%) | 50/50 (100%) | 99/100 (99%)
2 | 50/50 (100%) | 50/50 (100%) | 100/100 (100%)
3 | 50/50 (100%) | 50/50 (100%) | 100/100 (100%)
4 | 50/50 (100%) | 49/50 (99%) | 99/100 (99%)
5 | 48/50 (96%) | 50/50 (100%) | 98/100 (98%)
Average | 98.8% | 99.6% | 99.2%
Figure 7:  Weight visualization is performed by reshaping each row vector of
the weight matrix $\bf{W}_{1}$ into the shape of the input images: the orange
arrow represents a two-dimensional to one-dimensional reshaping and the green
arrow represents a one-dimensional to two-dimensional reshaping.

###  III.3 Weight visualization

In this subsection, we explain the method for weight visualization and the
meaning of positive and negative values in the visualized weight.

We describe the visualization method of the weight matrices
$\mathbf{W_{\mathrm{1}}}$ for investigating the extracted features.Yee,
Yamanaka, and Tagawa (2022) Since weights represent the importance of an input
signal, visualization of the weight matrices enables us to investigate the
important features of the input images.

An illustration of the weight visualization method used in this study is shown
in Fig. 7. As shown by the orange arrow in the figure, the two-dimensional
matrix of the input image is flattened into a one-dimensional vector and
connected to neurons in the input layer. This one-dimensional matrix
$\mathbf{x}$ is input into Eq. 1.

As explained in Sec. III.1, we set the number of neurons in the hidden layer
to five in this study, so the number of rows in the weight matrix
$\mathbf{W_{\mathrm{1}}}$ is five. Note that we refer to the five neurons as
$n1$, $n2$, $n3$, $n4$, and $n5$. As shown by the green arrow in the figure,
the weight vector that corresponds to each of these neurons is reshaped into a
two-dimensional matrix, which has the same size as an input image and is
visualized as a colormap. The colormaps of the reshaped weight vectors of the
neurons trained using Combination 5 are shown in Fig. 8. Note that we used
Combination 5 for the explanation because the extracted features in this case
are the clearest among all combinations. The results are similar for all
combinations.

Figure 8:  Weight visualization for the five neurons ($n1$–$n5$).

In Fig. 8, red, green, and blue indicate positive, zero, and negative values,
respectively. Although there are almost all zero values in $n1$, $n2$, and
$n3$, the distributions of the positive and negative values of $n4$ and $n5$
resemble the bubbles in the input images. Notably, the distributions of $n4$
and $n5$ are opposite to one another, with the distributions of the positive
values in $n4$ to those of the negative values in $n5$. This is because the
elements of $\bf{W}_{\textrm{out}}$ that correspond to $n4$ and $n5$, i.e.,
${W}_{\textrm{out},4}$ and ${W}_{\textrm{out},5}$, have different signs, with
${W}_{\textrm{out},4}$ having a negative value when ${W}_{\textrm{out},4}$ has
a positive value. Moreover, the distribution in $n4$ is more obvious than that
in $n5$. Thus, the discussion below is focused on $n4$.

In the image of the weight matrix for $n4$ shown in Fig. 8, positive elements
indicate “Fast” features and negative elements indicate “Slow” features. This
is explained by analyzing the calculation process as follows.

As shown in Fig. 4, there are two neurons in the output layer. We express
values for “Slow” using the superscript “s” and “Fast” using the superscript
“f”, i.e., $\mathbf{y_{\mathrm{pred}}}=y_{\mathrm{pred}}^{\mathrm{s}}$ for
“Slow” and $\mathbf{y_{\mathrm{pred}}}=y_{\mathrm{pred}}^{\mathrm{f}}$ for
“Fast” neurons.

When “Slow” is predicted, for the output layer,
$y_{\mathrm{pred}}^{\mathrm{s}}$ is larger than
$y_{\mathrm{pred}}^{\mathrm{f}}$, i.e., $y_{\mathrm{pred}}^{\mathrm{s}}$ is
closer to 1. Therefore, $q_{\mathrm{out}}^{\mathrm{s}}$ has to be larger than
$q_{\mathrm{out}}^{\mathrm{f}}$. From Eq. 3, $q_{\mathrm{out}}^{\mathrm{s}}$
can be expressed as

$q_{\mathrm{out}}^{\mathrm{s}}=W_{\mathrm{out,1}}^{\mathrm{s}}s_{\mathrm{1,1}}+W_{\mathrm{out,2}}^{\mathrm{s}}s_{\mathrm{1,2}}+\dots+W_{\mathrm{out,5}}^{\mathrm{s}}s_{\mathrm{1,5}}+b_{\mathrm{out}}^{\mathrm{s}}.$
(5)

In Eq. 5, since the $W_{\mathrm{out,4}}^{\mathrm{s}}$ obtained from the
training process has a positive value, $s_{\mathrm{1,4}}$ has to be as large
as possible to obtain the “Slow” prediction. Moreover, the relationship
between $s_{\mathrm{1,4}}$ and $q_{\mathrm{1,4}}$ is determined by the Leaky
ReLU function as shown in Eq. 2. Therefore, $q_{\mathrm{1,4}}$ has to be as
large as possible. From Eq. 1, $q_{\mathrm{1,4}}$ is expressed as

$q_{\mathrm{1,4}}=\mathbf{W_{\mathrm{1,4}}x}+b_{\mathrm{1,4}},$ (6)

where $\mathbf{W_{\mathrm{1,4}}}$ is the fourth row of the weight matrix
$\mathbf{W_{\mathrm{1}}}$ and $b_{\mathrm{1,4}}$ is the corresponding bias
element, which has an insignificant value in our case. From above,
$\mathbf{W_{\mathrm{1,4}}x}$ has to be as large as possible when “Slow” is
predicted.

Since shadowgraphy was used to capture the image sequences of the laser-
induced microjet, the pixels in an image covered by a vapor bubble or
cavitation bubbles are dark, and thus have low normalized intensity values
($x_{i,j}$, the elements of $\mathbf{x}$, are close to 0). On the other hand,
due to the backlight, the pixels not covered by bubbles are bright and have
high normalized intensity values ($x_{i,j}$ are close to 1). Therefore, the
elements of $\mathbf{W_{\mathrm{1,4}}}$ that correspond to the area covered by
the bubbles in the “Slow” images but not that in the “Fast” images have
negative values, so these values become zero when multiplied by a “Slow” image
to increase the value of $\mathbf{W_{\mathrm{1,4}}x}$. In contrast, for a
“Fast” image, these values would remain as negative values, thereby reducing
the value of $\mathbf{W_{\mathrm{1,4}}x}$. Therefore, the elements with
negative values, which are indicated by the blue pixels in Fig. 8, show the
positions of the bubbles in the “Slow” images. On the other hand, the elements
of $\mathbf{W_{\mathrm{1,4}}}$ that correspond to the area covered by the
bubbles in the “Fast” images but not in the “Slow” images have to be positive
values so that these values would become zero when multiplied by a “Fast”
image to reduce the value of $\mathbf{W_{\mathrm{1,4}}x}$. However, for a
“Slow” image, these values would remain as negative values, thereby increasing
the value of $\mathbf{W_{\mathrm{1,4}}x}$. These elements with positive
values, which are indicated by the red pixels in Fig. 8, show the positions of
the bubbles in the “Fast” images.

##  IV Results and discussion

In this section, we demonstrate the relationship between the jet velocity
obtained from the experiments and that predicted by the FNN in Sec. IV.1.
Next, in Sec. IV.2, we explain the influence of a vapor bubble and cavitation
bubbles on the jet velocity. Finally, we report the extracted features related
to the distribution of the cavitation positions obtained from the weight
visualization in Sec. IV.3.

###  IV.1 The relationship between the experimental and predicted jet velocity

In this subsection, we investigate and discuss the relationship between the
jet velocity and the predicted values $\mathbf{q_{\mathrm{out}}}$ obtained
from the FNN.

As described in Sec. III.1, the classification of “Fast” and “Slow” is based
on the vector $\mathbf{y_{\mathrm{pred}}}$, which is normalized by the sigmoid
function, as shown in Eq. 4. The sigmoid function saturates extreme negative
and positive values to zero and one, respectively. Therefore, we used the
vector of $\mathbf{q_{\mathrm{out}}}$, which is before normalization, in the
investigation of the relationship between the experimental values and the
predicted values. The relationship between the experimental jet velocity
$U_{\mathrm{\mathrm{jet}}}$ and the $q_{\mathrm{out}}^{\mathrm{f}}$ is shown
in Fig. 9. Note that $\mathbf{q_{\mathrm{out}}}$ has two values, but since
both of them show similar trends, we use the values that indicate “Fast”.

We also tested data other than the “Fast” and “Slow” data of the model used in
this study as interpolation test data and obtained
$q_{\mathrm{out}}^{\mathrm{f}}$ for each data point. Then, the relationship
between the jet velocity and the predicted values for the interpolation test
data was plotted as a black plot in Fig. 9. Remarkably, the correlation
coefficient is 0.912, which is significantly higher than that obtained using
the total area of the cavitation bubbles shown in Fig. 3(b). This result
indicates that the FNN can predict the jet velocity more accurately, even for
interpolation test data, by considering not only the total area of cavitation
bubbles but also other features. Therefore, we hereafter extract the features
that the FNN observes to understand the important factors that influence the
jet velocity.

Figure 9:  The experimental jet velocity $U_{\mathrm{\mathrm{jet}}}$ and the
predicted values $q_{\mathrm{out}}^{\mathrm{f}}$ for the FNN. The blue and red
areas represent the data of the fastest 25% (“Fast”) and the slowest 25%
(“Slow”) of the jet velocity, respectively.

###  IV.2 The influence of vapor bubbles vs. cavitation

Figure 10:  (a) The experimental jet velocity $U_{\mathrm{\mathrm{jet}}}$ vs.
predicted values $q_{\mathrm{out,vap}}^{\mathrm{f}}$ for the FNN trained with
images of vapor bubbles and (b) the experimental jet velocity
$U_{\mathrm{jet}}$ vs. predicted values $q_{\mathrm{out,cav}}^{\mathrm{f}}$
for the FNN trained with images of cavitation bubbles.

In this subsection, we compare the effects of a vapor bubble and cavitation
bubbles on the jet velocity by comparing the relationship between the
experimental jet velocity and the predicted values from the FNN. We compare
the training of the FNN using cropped images of only a vapor bubble and only
cavitation bubbles, respectively. Fig. 10(a) shows the relationship between
the jet velocity $U_{\mathrm{\mathrm{jet}}}$ and the predicted values
$q_{\mathrm{out,vap}}^{\mathrm{f}}$ when training only with a vapor bubble. In
addition, Fig. 10(b) shows the relationship between the jet velocity
$U_{\mathrm{\mathrm{jet}}}$ and the predicted values
$q_{\mathrm{out,cav}}^{\mathrm{f}}$ when training only with cavitation
bubbles. The correlation coefficients are 0.688 and 0.883, respectively, and
the averaged total accuracies of the classification are 90.0% and 98.4%,
respectively. These results suggest that cavitation bubbles at $t=20$ µs have
a greater impact on the prediction of the jet velocity than vapor bubbles.

###  IV.3 Weight visualization: distribution of the cavitation position

Figure 11:  (a) Input images of cases with similar sizes of cavitation bubbles
occurring at different positions: (i) $S_{\mathrm{cav}}=0.0302$ mm2,
$U_{\mathrm{jet}}=26.3$ m/s, (ii) $S_{\mathrm{cav}}=0.0294$ mm2,
$U_{\mathrm{jet}}=17.9$ m/s and (b) the binarized image of (a)(ii) with the
blue dot for centroid of the cavitation bubbles. Figure 12:  (a) The
distribution of the centroids of the cavitation bubbles, (b) the colormap of
the weight matrix $\mathbf{W_{\mathrm{1,4}}}$ for classification of the jet
velocity, and (c) an example image of (i) “Slow” and (ii) “Fast” images, which
include cavitation bubbles.

In this subsection, we describe an investigation of the feature extraction
results and the effects of the cavitation position on the jet velocity. For
this purpose, we investigate the weight matrix $\mathbf{W_{\mathrm{1}}}$,
which determines the predicted values $\mathbf{q_{\mathrm{out}}}$.

To demonstrate the effects of the cavitation positions, we compare cases with
cavitation bubbles of similar sizes occurring at different places in Figs.
11(a): (i) $S_{\mathrm{cav}}=0.0302$ mm2 and (ii) $S_{\mathrm{cav}}=0.0294$
mm2. Remarkably, the jet velocity for the case of Fig. 11(a)(i) is
$U_{\mathrm{jet}}=26.3$ m/s ($q_{\mathrm{out,cav}}^{\mathrm{f}}=0.332$), which
is higher than the value of $U_{\mathrm{jet}}=17.9$ m/s
($q_{\mathrm{out,cav}}^{\mathrm{f}}=-1.70$) in Fig. 11(a)(ii). These cases
show that the jet velocity greatly depends on the cavitation positions. To
further investigate the effects of the cavitation positions, the centroids of
the cavitation bubbles are obtained from the binarized image as shown in Fig.
11(b). The centroids of the cavitation bubbles for all data (1000 data) are
shown in Fig. 12(a). From this figure, we identified regions where the
cavitation is often located. Regions where cavitation bubbles are frequently
generated are marked as Regions [1], [2], [3], and [4] in the figure. Out of
the 1000 data, cavitation bubbles are observed in 931 (93.1%) data, including
250 out of the 250 (100%) “Fast” data, 188 out of the 250 (75.2%) “Slow” data,
and 493 out of the 500 (98.6%) data that are neither “Fast” nor “Slow”. For
comparison with the extracted features, Fig. 12(b) shows an image of the
weight matrix $\mathbf{W_{\mathrm{1,4}}}$ obtained using the method described
in Sec. III.3. Note that positive weights (red region) indicate “Fast”
features and negative weights (blue region) indicate “Slow” features, as
explained in Sec. III.3. Comparing Figs. 12(a) and (b), we verify the “Fast”
weight in Regions [1]–[3]. Note that the “Fast” weight is also shown around
the vapor bubble. This demonstrates that the larger the vapor bubble, the
faster the jet velocity, which is consistent with a finding of the previous
study by Kawamoto et al.Kawamoto  _et al._ (2016)

Although the “Fast” weight can be seen in Regions [1]–[3], “Fast” features are
not seen in Region [4] in Fig. 12(b). To investigate this region, we show an
example of a “Slow” input image, although the cavitation bubbles are generated
in Region [4] in Fig. 12(c)(i). In fact, cavitation bubbles occurred in Region
[4] in 152 out of the 250 (60.8%) “Slow” images. This shows that there are
cavitation positions that are not related to the increase in the jet velocity.
Note that, as can be seen from the example of a “Fast” input image indicating
the generation of cavitation bubbles at all defined regions shown in Fig.
12(c)(ii), we can also observe the bubble generation at Region [4] even in the
“Fast” image. Therefore, Figs. 12(c)(i) and (ii) indicate that the bubble
generations in Regions [1]–[3] affect the increase in jet velocity, while the
generation at Region [4] is hardly related to the jet velocity.

Here, the weight diagram shown in Fig. 12(b) indicates the overall evaluation
of the frequency of cavitation generation in each region and the cavitation
area. The more frequent the cavitation generation in each region, the greater
the value of the weight element at that position. Moreover, the larger the
cavitation area, the higher the possibility that the weight of the region is
large. In order to distinguish between the cavitation frequency and the
cavitation area at each position, the area of each cavitation is represented
by the scale of the colorbar and plotted at the centroid in Fig. 13. As shown
in the subfigures, we divide the data into 4 groups, each 250 data, in order
of ascending velocity. The slowest 25% are Group 1 (“Slow”), the next slowest
25% are Group 2, the second fastest 25% are Group 3, and the fastest 25% are
Group 4 (“Fast”). Fig. 13 shows that the cavitation area increases with the
jet velocity. This is consistent with what was described in the previous study
of Kiyama et al. Kiyama  _et al._ (2016)

Fig. 13 also shows that the higher the jet velocity, the more cavitations
occur in Regions [1] and [2]. We show the number of cavitation generations in
each region and the velocity range in Table 4. In Regions [1] and [2], the
number of cavitation bubbles increases with increasing velocity. On the other
hand, the number does not necessarily increase in Regions [3] and [4]. This
suggests that the effect of the velocity increase is more significant when
cavitation occurs closer to the laser focus position, where it is farther from
the gas-liquid interface. Note that, in Region [3], cavitation bubbles with a
large area are generated at high velocities, while in Region [4], small
cavitation bubbles are generated regardless of the velocity range, so the
weight is red in Region [3], and there is almost no weight in Region [4], as
shown in Fig. 12(b).

Here, we discuss the possible mechanisms of how the cavitation position
affects the jet acceleration. Kiyama et al.Kiyama  _et al._ (2016) reported
that after the propagation of the compression wave, the fluid pressure becomes
positive, while after the propagation of the expansion wave, it becomes
negative. Moreover, the compression waves are converted into expansion waves
and vice versa at the gas-liquid interface. When the waves are in close
proximity to the gas-liquid interface, the time interval between the
propagation of compression and expansion waves is minimal. In contrast, the
farther away it is from the gas-liquid interface, the longer the time interval
is. Therefore, the greater the distance from the gas-liquid interface, the
longer the time the fluid experiences the tensile stress, which means a
greater pressure impulse. As reported by Hayasaka $et\ al$.Hayasaka, Kiyama,
and Tagawa (2017), the larger the pressure impulse experienced by the liquid,
the larger the jet velocity. Consequently, the cavitation positions farther
from the gas-liquid interface, where the fluid experiences the larger pressure
impulse, makes a more significant contribution to the velocity increase. For
further understanding, we need to observe the propagation behavior of a shock
wave with cavitation bubbles. For this purpose, we suggest the use of the
background-oriented schlieren (BOS) technique Venkatakrishnan and Meier
(2004); Yamamoto  _et al._ (2022) to visualize the shock waves in the
microtube.

Figure 13:  The distribution of the centroids of the cavitation bubbles with a colormap of the cavitation area $\mathbf{S_{\mathrm{cav}}}$ for 250 data at each velocity. Groups 1,2,3 and 4 are in order of decreasing velocity. Table 4:  The number of cavitation bubbles in each region and for each velocity range. | Region 1 | Region 2 | Region 3 | Region 4 | Total
---|---|---|---|---|---
Group 1 | 3 (0.865%) | 4 (1.15%) | 177 (51.0%) | 163 (47.0%) | 347
Group 2 | 39 (8.52%) | 17 (3.71%) | 262 (57.2%) | 140 (30.6%) | 458
Group 3 | 76 (13.8%) | 30 (5.44%) | 307 (55.7%) | 138 (25.0%) | 551
Group 4 | 227 (36.7%) | 50 (8.09%) | 285 (46.1%) | 56 (9.06%) | 618

##  V Conclusions and outlook

In this study, we have conducted 1,000 experiments to obtain corresponding
images and jet velocities. We have investigated the relationship between the
amount of cavitation bubbles and the jet velocity and found that there was
still non-negligible variation.

To investigate the reason for this, we utilized explainable artificial
intelligence (XAI), namely a feedforward neural network (FNN), to make
predictions and classify the jet velocity from images of cavitation bubbles.
The results showed that at $t=20$ µs, when cavitation was most likely to
occur, the averaged total accuracy was as high as 99.2%. This confirmed that
the FNN extracted clear features related to the jet velocity.

In order to clarify the extracted features, we visualized the calculation
process in the FNN. The results suggested that the prediction values obtained
by training only cavitation bubbles had a higher correlation with the jet
velocity than those using only vapor bubbles. This demonstrates that
cavitation bubbles have a greater impact on predicting the jet velocity than
vapor bubbles. In addition, we verified that the correlation coefficient
between the experimental jet velocity and the predicted values obtained from
training only cavitation bubbles was higher than that between the jet velocity
and the area of the cavitation bubbles. This indicates that the FNN extracted
features other than the area of cavitation.

Therefore, we visualized the “weights” that represent the important signals in
the input images to investigate the extracted features. As a result, we
clarified that there are cavitation positions that do or do not affect the jet
velocity. Further investigation suggested that cavitation that occurs closer
to the laser focus position has a higher acceleration effect. We argue that
this is because in close proximity to the laser focus position, the greater
the pressure impulse.

We now summarize three possible contributions that could be developed from the
findings of this study. The first is the development of the jet velocity
prediction method. We successfully predicted the jet velocity from the images
of jet generation using machine learning, indicating that machine learning can
be used as an application for predicting jet velocities. The second
contribution is the development of XAI used in this study. We successfully
extracted features from a large amount of data acquired from experiments,
verifying the usefulness of XAI as an aid in the interpretation of fluid
experiment images. We believe that XAI is a highly applicable method for
visualizing important features, especially in further experimental fluid
systems with many parameters. Finally, the third contribution is the
usefulness of this approach not only for high-speed microjets but also for
various phenomena involving cavitation bubbles. For this, we need to further
consider the physical meaning of the critical cavitation position. Cavitation
bubbles and the jet velocity are related to the propagation of shock waves.
Therefore, we need to observe the propagation behavior of a shock wave with
cavitation bubbles. For this purpose, a possible method to visualize the shock
waves in a microtube is the background-oriented schlieren (BOS)
techniqueVenkatakrishnan and Meier (2004); Yamamoto  _et al._ (2022).

###### Acknowledgements.

This work was funded by the Japan Society for the Promotion of Science (Grant
Nos. 20H00223, 20H00222, and 20K20972) and the Japan Science and Technology
Agency PRESTO (Grant No. JPMJPR21O5). The authors would also like to thank Dr.
Masaharu Kameda (Professor, Tokyo University of Agriculture and Technology)
and Dr. Akinori Yamanaka (Professor, Tokyo University of Agriculture and
Technology) for helpful discussions and comments.

## References

## References

  * Nir  _et al._ (2003) Y. Nir, A. Paz, E. Sabo, and I. Potasman, “Fear of injections in young adults: prevalence and associations,” The American Journal of Tropical Medicine and Hygiene 68, 341–344 (2003).
  * McLenon and Rogers (2019) J. McLenon and M. A. Rogers, “The fear of needles: A systematic review and meta-analysis,” Journal of Advanced Nursing 75, 30–42 (2019).
  * Deacon and Abramowitz (2006) B. Deacon and J. Abramowitz, “Fear of needles and vasovagal reactions among phlebotomy patients,” Journal of Anxiety Disorders 20, 946–960 (2006).
  * Kermode (2004) M. Kermode, “Unsafe injections in low-income country health settings: need for injection safety promotion to prevent the spread of blood-borne viruses,” Health Promotion International 19, 95–103 (2004).
  * Weiss and McMichael (2004) R. A. Weiss and A. J. McMichael, “Social and environmental risk factors in the emergence of infectious diseases,” Nature Medicine 10, S70–S76 (2004).
  * res (2003) Centers for Disease Control and Prevention.: Incorporating HIV prevention into the medical care of persons living with HIV. Recommendations of CDC, the Health Resources and Services Administration, the National Institutes of Health, and the HIV Medicine Association of the Infectious Diseases Society of America. Morbidity and Mortality Weekly Report. 52, 1–24 (2003).
  * Mitragotri (2005) S. Mitragotri, “Immunization without needles,” Nature Reviews Immunology 5, 905–916 (2005).
  * Stachowiak  _et al._ (2009) J. C. Stachowiak, T. H. Li, A. Arora, S. Mitragotri, and D. A. Fletcher, “Dynamic control of needle-free jet injection,” Journal of Controlled Release 135, 104–112 (2009).
  * Schoppink and Rivas (2022) J. Schoppink and D. F. Rivas, “Jet injectors: Perspectives for small volume delivery with lasers,” Advanced Drug Delivery Reviews , 114109 (2022).
  * Daly  _et al._ (2020) C. Daly, N. A. Molodecky, M. Sreevatsava, A. D. Belayneh, S. A. Chandio, J. Partridge, A. Shaikh, M. Laghari, J. Agbor, R. M. Safdar, _et al._ , “Needle-free injectors for mass administration of fractional dose inactivated poliovirus vaccine in karachi, pakistan: a survey of caregiver and vaccinator acceptability,” Vaccine 38, 1893–1898 (2020).
  * Harding  _et al._ (2002) L. Harding, A. Adeniyi, R. Everson, S. Barker, D. Ralph, and A. Baranowski, “Comparison of a needle-free high-pressure injection system with needle-tipped injection of intracavernosal alprostadil for erectile dysfunction,” International Journal of Impotence Research 14, 498–501 (2002).
  * Giudice and Campbell (2006) E. L. Giudice and J. D. Campbell, “Needle-free vaccine delivery,” Advanced Drug Delivery Reviews 58, 68–89 (2006).
  * Mitragotri (2006) S. Mitragotri, “Current status and future prospects of needle-free liquid jet injectors,” Nature Reviews Drug Discovery 5, 543–548 (2006).
  * Arora  _et al._ (2007) A. Arora, I. Hakim, J. Baxter, R. Rathnasingham, R. Srinivasan, D. A. Fletcher, and S. Mitragotri, “Needle-free delivery of macromolecules across the skin by nanoliter-volume pulsed microjets,” Proceedings of the National Academy of Sciences 104, 4255–4260 (2007).
  * Tagawa  _et al._ (2012) Y. Tagawa, N. Oudalov, C. W. Visser, I. R. Peters, D. van der Meer, C. Sun, A. Prosperetti, and D. Lohse, “Highly focused supersonic microjets,” Physical Review X 2, 031002 (2012).
  * Tagawa  _et al._ (2013) Y. Tagawa, N. Oudalov, A. El Ghalbzouri, C. Sun, and D. Lohse, “Needle-free injection into skin and soft matter with highly focused microjets,” Lab on a Chip 13, 1357–1363 (2013).
  * Miyazaki  _et al._ (2021) Y. Miyazaki, M. Usawa, S. Kawai, J. Yee, M. Muto, and Y. Tagawa, “Dynamic mechanical interaction between injection liquid and human tissue simulant induced by needle-free injection of a highly focused microjet,” Scientific Reports 11, 1–10 (2021).
  * Kiyama  _et al._ (2019) A. Kiyama, N. Endo, S. Kawamoto, C. Katsuta, K. Oida, A. Tanaka, and Y. Tagawa, “Visualization of penetration of a high-speed focused microjet into gel and animal skin,” Journal of Visualization 22, 449–457 (2019).
  * Franco-Gómez  _et al._ (2021) A. Franco-Gómez, H. Onuki, Y. Yokoyama, Y. Nagatsu, and Y. Tagawa, “Effect of liquid elasticity on the behaviour of high-speed focused jets,” Experiments in Fluids 62, 1–15 (2021).
  * Van Hoeve  _et al._ (2010) W. Van Hoeve, S. Gekle, J. H. Snoeijer, M. Versluis, M. P. Brenner, and D. Lohse, “Breakup of diminutive Rayleigh jets,” Physics of Fluids 22, 122003 (2010).
  * Scroggs and Settles (1996) S. Scroggs and G. Settles, “An experimental study of supersonic microjets,” Experiments in Fluids 21, 401–409 (1996).
  * Choo and Kang (2007) Y. Choo and B. Kang, “The effect of jet velocity profile on the characteristics of thickness and velocity of the liquid sheet formed by two impinging jets,” Physics of Fluids 19, 112101 (2007).
  * Thoroddsen  _et al._ (2009) S. T. Thoroddsen, K. Takehara, T. Etoh, and C.-D. Ohl, “Spray and microjets produced by focusing a laser pulse into a hemispherical drop,” Physics of Fluids 21, 112101 (2009).
  * Kamamoto, Onuki, and Tagawa (2021) K. Kamamoto, H. Onuki, and Y. Tagawa, “Drop-on-demand painting of highly viscous liquids,” Flow 1, E6 (2021).
  * Kiyama  _et al._ (2016) A. Kiyama, Y. Tagawa, K. Ando, and M. Kameda, “Effects of a water hammer and cavitation on jet formation in a test tube,” Journal of Fluid Mechanics 787, 224–236 (2016).
  * Hayasaka, Kiyama, and Tagawa (2017) K. Hayasaka, A. Kiyama, and Y. Tagawa, “Effects of pressure impulse and peak pressure of a shockwave on microjet velocity in a microchannel,” Microfluidics and Nanofluidics 21, 166 (2017).
  * Daily  _et al._ (2014) J. Daily, J. Pendlebury, K. Langley, R. Hurd, S. Thomson, and T. Truscott, “Catastrophic cracking courtesy of quiescent cavitation,” Physics of Fluids 26, 091107 (2014).
  * Ishikawa, Nishida, and Tagawa (2022) T. Ishikawa, H. Nishida, and Y. Tagawa, “Numerical investigation on influence of number of bubbles on laser-induced microjet,” Water 14, 3707 (2022).
  * Brunton, Noack, and Koumoutsakos (2020) S. L. Brunton, B. R. Noack, and P. Koumoutsakos, “Machine learning for fluid mechanics,” Annual Review of Fluid Mechanics 52, 477–508 (2020).
  * Colvert, Alsalman, and Kanso (2018) B. Colvert, M. Alsalman, and E. Kanso, “Classifying vortex wakes using neural networks,” Bioinspiration & Biomimetics 13, 025003 (2018).
  * Ling, Kurzawski, and Templeton (2016) J. Ling, A. Kurzawski, and J. Templeton, “Reynolds averaged turbulence modelling using deep neural networks with embedded invariance,” Journal of Fluid Mechanics 807, 155–166 (2016).
  * Kutz (2017) J. N. Kutz, “Deep learning in fluid dynamics,” Journal of Fluid Mechanics 814, 1–4 (2017).
  * Folden and Aschmoneit (2023) T. S. Folden and F. J. Aschmoneit, “A classification and review of cavitation models with an emphasis on physical aspects of cavitation,” Physics of Fluids 35 (2023).
  * Williams, Kevrekidis, and Rowley (2015) M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, “A data–driven approximation of the koopman operator: Extending dynamic mode decomposition,” Journal of Nonlinear Science 25, 1307–1346 (2015).
  * Jovanović, Schmid, and Nichols (2014) M. R. Jovanović, P. J. Schmid, and J. W. Nichols, “Sparsity-promoting dynamic mode decomposition,” Physics of Fluids 26, 024103 (2014).
  * Abdi and Williams (2010) H. Abdi and L. J. Williams, “Principal component analysis,” Wiley Interdisciplinary Reviews: Computational Statistics 2, 433–459 (2010).
  * Kerschen  _et al._ (2005) G. Kerschen, J.-C. Golinval, A. F. Vakakis, and L. A. Bergman, “The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: an overview,” Nonlinear Dynamics 41, 147–169 (2005).
  * Vennemann and Rösgen (2020) B. Vennemann and T. Rösgen, “A dynamic masking technique for particle image velocimetry using convolutional autoencoders,” Experiments in Fluids 61, 1–11 (2020).
  * Sekar  _et al._ (2019) V. Sekar, Q. Jiang, C. Shu, and B. C. Khoo, “Fast flow field prediction over airfoils using deep learning approach,” Physics of Fluids 31, 057103 (2019).
  * Liu  _et al._ (2020) B. Liu, J. Tang, H. Huang, and X.-Y. Lu, “Deep learning methods for super-resolution reconstruction of turbulent flows,” Physics of Fluids 32, 025105 (2020).
  * Cai  _et al._ (2019) S. Cai, S. Zhou, C. Xu, and Q. Gao, “Dense motion estimation of particle images via a convolutional neural network,” Experiments in Fluids 60, 1–16 (2019).
  * Jin  _et al._ (2018) X. Jin, P. Cheng, W.-L. Chen, and H. Li, “Prediction model of velocity field around circular cylinder over various Reynolds numbers by fusion convolutional neural networks based on pressure on the cylinder,” Physics of Fluids 30, 047105 (2018).
  * Arrieta  _et al._ (2020) A. B. Arrieta, N. Díaz-Rodríguez, J. Del Ser, A. Bennetot, S. Tabik, A. Barbado, S. García, S. Gil-López, D. Molina, R. Benjamins, _et al._ , “Explainable artificial intelligence (XAI): Concepts, taxonomies, opportunities and challenges toward responsible AI,” Information Fusion 58, 82–115 (2020).
  * Adadi and Berrada (2018) A. Adadi and M. Berrada, “Peeking inside the black-box: a survey on explainable artificial intelligence (XAI),” IEEE Access 6, 52138–52160 (2018).
  * Gunning (2017) D. Gunning, “Explainable artificial intelligence (XAI),” Defense Advanced Research Projects Agency (DARPA) 2, 1 (2017).
  * Gunning  _et al._ (2019) D. Gunning, M. Stefik, J. Choi, T. Miller, S. Stumpf, and G.-Z. Yang, “XAI—explainable artificial intelligence,” Science Robotics 4, eaay7120 (2019).
  * Yee  _et al._ (2023) J. Yee, D. Igarashi, S. Miyatake, and Y. Tagawa, “Prediction of the morphological evolution of a splashing drop using an encoder–decoder,” Machine Learning: Science and Technology 4, 025002 (2023).
  * Yee, Yamanaka, and Tagawa (2022) J. Yee, A. Yamanaka, and Y. Tagawa, “Image features of a splashing drop on a solid surface extracted using a feedforward neural network,” Physics of Fluids 34, 013317 (2022).
  * Akhatov  _et al._ (2001) I. Akhatov, O. Lindau, A. Topolnikov, R. Mettin, N. Vakhitova, and W. Lauterborn, “Collapse and rebound of a laser-induced cavitation bubble,” Physics of Fluids 13, 2805–2819 (2001).
  * Caupin and Herbert (2006) F. Caupin and E. Herbert, “Cavitation in water: a review,” Comptes Rendus Physique 7, 1000–1017 (2006).
  * Rosenblatt (1958) F. Rosenblatt, “The perceptron: a probabilistic model for information storage and organization in the brain.” Psychological Review 65, 386 (1958).
  * Hornik, Stinchcombe, and White (1989) K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural Networks 2, 359–366 (1989).
  * Maas, Hannun, and Ng (2013) A. L. Maas, A. Y. Hannun, and A. Y. Ng, “Rectifier nonlinearities improve neural network acoustic models,” Proc. ICML (2013).
  * Géron (2019) A. Géron, “Hands-on machine learning with scikit-learn, keras, and tensorflow: Unsupervised learning techniques. O’Reilly Media, Incorporated,” O’Reilly Media, Incorporated (2019).
  * Li  _et al._ (2014) M. Li, T. Zhang, Y. Chen, and A. J. Smola, “Efficient mini-batch training for stochastic optimization,” Proceedings of the 20th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining , 661–670 (2014).
  * Prechelt (2002) L. Prechelt, “Early stopping-but when?” in  _Neural Networks: Tricks of the trade_ (Springer, 2002) pp. 55–69.
  * Kawamoto  _et al._ (2016) S. Kawamoto, K. Hayasaka, Y. Noguchi, and Y. Tagawa, “Volume of a laser-induced microjet,” Transactions of the JSME (in Japanese) 82, 16–00094 (2016).
  * Venkatakrishnan and Meier (2004) L. Venkatakrishnan and G. Meier, “Density measurements using the background oriented schlieren technique,” Experiments in Fluids 37, 237–247 (2004).
  * Yamamoto  _et al._ (2022) S. Yamamoto, T. Shimazaki, A. Franco-Gómez, S. Ichihara, J. Yee, and Y. Tagawa, “Contactless pressure measurement of an underwater shock wave in a microtube using a high-resolution background-oriented schlieren technique,” Experiments in Fluids 63, 142 (2022).


#  Rotating Bose gas dynamically enters the lowest Landau level

Vaibhav Sharma<EMAIL_ADDRESS>Erich J Mueller<EMAIL_ADDRESS>Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, New
York

###### Abstract

Motivated by recent experiments, we model the dynamics of a condensed Bose gas
in a rotating anisotropic trap, where the equations of motion are analogous to
those of charged particles in a magnetic field. As the rotation rate is ramped
from zero to the trapping frequency, the condensate stretches along one
direction and is squeezed along another, becoming long and thin. When the trap
anisotropy is slowly switched off on a particular timescale, the condensate is
left in the lowest Landau level. We use a time dependent variational approach
to quantify these dynamics and give intuitive arguments about the structure of
the condensate wavefunction. This preparation of a lowest Landau level
condensate can be an important first step in realizing bosonic analogs of
quantum Hall states.

##  I Introduction

For many years the ultracold atom community has been attempting to use neutral
atoms to explore the physics of 2D electron gases in magnetic fields. The most
exciting possibility is realizing fractional quantum Hall physics and
observing non-abelian phases of matter [1]. A recent experiment at MIT [2]
made an important step towards achieving this goal. They found that when
rotating a Bose gas in an anisotropic trap, the resulting non-equilibrium
dynamics can drive the atoms into the lowest Landau level, a pre-requisite for
observing the quantum Hall effect. Motivated by this observation, we model the
dynamics of a rotating Bose gas in an anisotropic trap using a time dependent
variational wavefunction. We characterize the protocols which drive the atoms
into the lowest Landau level.

There are a number of techniques for using neutral atoms to emulate the
properties of charged particles in a magnetic field: These include optical
dressing [3], shaking [4] and rotating [5]. Here we focus on rotation, where
the connection with electronic systems can be understood classically: In the
rotating frame, the Coriolis force, $F_{\rm c}=2m\,\vec{v}\times\vec{\Omega}$
has the same form as the Lorentz force $F_{\rm L}=q\,\vec{v}\times\vec{B}$ on
a charged particle $q$ in a magnetic field $\vec{B}$. Here $\vec{\Omega}$ is a
vector whose magnitude is the rotation frequency, and direction is the axis of
rotation (which we take to be $\hat{z}$). One can formally identify
$q\vec{B}=2m\,\vec{\Omega}$. The rotating system also experiences an outward
centrifugal force, which reduces the strength of the in-plane harmonic
trapping potential. If the trapping frequency equals the rotation frequency,
$\omega=\Omega$, the effective trapping frequency vanishes, and in the
rotating frame the neutral atoms have equations of motion which are identical
to electrons in a uniform magnetic field. The quantum mechanical eigenstates
correspond to highly degenerate Landau levels, whose energies are separated by
$2\hbar\Omega$. The resulting behavior is quite rich; for example, rotating
superfluids display vortex lattices similar to those of type-2 superconductors
in magnetic fields [6, 7, 8] and also exhibit a rotational Meissner effect
[9].

It is meaningless to talk of the rotation frequency of a rotationally
symmetric potential. Thus, following the experimental approach in Ref. [2], we
include a small rotating elliptical deformation which stirs the cloud. In the
rotating frame this potential takes the form $V_{\varepsilon}=-\varepsilon
m\omega^{2}xy$, and can be described in terms of an electric field
$q\vec{E}=-\vec{\nabla}V_{\varepsilon}$. In the steady state, the bosons will
experience a Hall drift with a velocity proportional to
$\vec{E}\times\vec{B}$, like classical charged particles in crossed electric
and magnetic fields. This drift “squeezes” the gas into a long snake-like
shape, as shown in Fig. 1. Remarkably, one can interpret the atomic state as
being confined to a distorted Landau Level.

If the trap deformation is turned off at the appropriate time, the
wavefunction of each atom ends up in the traditional lowest Landau level. The
turn-off must occur adiabatically on a timescale which is moderately larger
than $1/(\varepsilon\omega)$, which is the characteristic dynamical time. If
the turn-off is too slow, however, non-linear terms drive the system back out
of the lowest Landau level. This can also be interpreted in terms of
transitions between levels which cross during the dynamics, and requires that
the deformation is turned off on timescale which is small compared to
$8/(\varepsilon\omega)\ln(8/\varepsilon)$.

It is worth noting that the number of single-particle states in the Landau
level is proportional to the area of the system, and it represents a large
phase space. One cannot equate being in the lowest Landau level with being in
a quantum Hall state. Nonetheless, driving the system into the lowest Landau
level is an important first step in realizing a cold atom quantum Hall effect.
Furthermore, the dynamics in the lowest Landau level are quite rich, and are a
useful area of study on their own [10, 11].

There have been numerous other works exploring rotating Bose gases in order to
access the quantum Hall regime [1]. The lowest Landau level was nearly
achieved by Schweikhard et al. [12] by using an evaporation technique. In the
few-body context, Gemelke et al. [13] found a driving protocol which not only
produced lowest-Landau level wavefunctions, but even analogs of the Laughlin
state. Stirring protocols similar to Ref. [2] have been previously used, but
the lack of in-situ imaging meant they could not study the physics explored
here [14, 15].

Our paper is organized as follows. In Sec. II, we introduce the model. In Sec.
III, we write down our time dependent wavefunction ansatz to model the
dynamics. In Sec. IV, we use the time dependent variational principle to
derive the equations of motion for the wavefunction parameters. In Sec. V, we
solve for the dynamics of the system and interpret the results.

##  II Model

We consider bosons at zero temperature in a harmonic trap with frequencies
$(\omega_{x},\omega_{y},\omega_{z})$. For $\omega_{x},\omega_{y}=\omega$ and
$\omega_{z}\gg\omega$, the z-motion is frozen and one effectively has a two
dimensional Bose gas.

We add a rotating deformation potential given by $V_{r}(t)=-\varepsilon
m\omega^{2}X(t)Y(t)$ where, $X(t)=x\cos{\phi(t)}+y\sin{\phi(t)}$ and
$Y(t)=-x\sin{\phi(t)}+y\cos{\phi(t)}$. We assume that the rotation rate
$\Omega(t)=d\phi/dt$ is slowly varying. The strength of the deformation is
given by the dimensionless variable $\varepsilon$. Note, if $\phi$ is phase
shifted by $\pi/4$, then $V_{r}(t)\to\varepsilon
m\omega^{2}(X(t)^{2}-Y(t)^{2})/2$, which is a form used in other papers [2,
10].

In the mean field regime, the condensate wavefunction $\psi$ obeys the 2D
Gross-Pitaevskii equation (GPE),

$i\hbar\frac{d\psi}{dt}=\left(\frac{-\hbar^{2}\nabla^{2}}{2m}+V_{r}(t)+\frac{m\omega^{2}(x^{2}+y^{2})}{2}+g|\psi|^{2}\right)\psi$
(1)

Here $\nabla^{2}$ is the two dimensional laplacian and $g$ denotes the
interaction strength. We approximate the interactions to be contact
interactions such that $g=4\pi\hbar^{2}a_{s}/(md_{\perp})$. $a_{s}$ is the
s-wave scattering length and $d_{\perp}\sim\sqrt{\hbar/m\omega_{z}}$ is the
size of the wavefunction in the z-direction [2].

We can transform into the rotating frame by taking $\psi\to
e^{-i\phi(t)L_{z}/\hbar}\psi$ where the angular momentum operator along z is
given by $L_{z}=i\hbar(x\partial_{y}-y\partial_{x})$. The rotating deformation
becomes time-independent and the GPE reads,

$i\hbar\frac{d\psi}{dt}=\left(\frac{-\hbar^{2}\nabla^{2}}{2m}+V_{\rm
eff}(x,y)-\Omega(t)L_{z}+g|\psi|^{2}\right)\psi$ (2)

where,

$V_{\rm eff}(x,y)=\frac{m\omega^{2}(x^{2}+y^{2})}{2}-\varepsilon
m\omega^{2}xy.$ (3)

We further adimensionalize by scaling all length scales by
$\sqrt{\hbar/(2m\omega)}$ and time by $1/\omega$. We also scale $\psi$ by
$1/\sqrt{N}$ such that $\int d^{2}r_{\perp}|\psi|^{2}=1$, where $N$ is the
total number of particles. The resulting adimensionalized equation is

$i\frac{d\psi}{dt}=\left(-\nabla^{2}+\frac{x^{2}}{4}+\frac{y^{2}}{4}-\frac{\varepsilon
xy}{2}-r(t)L_{z}+\tilde{g}N|\psi|^{2}\right)\psi$ (4)

Here $L_{z}=i(x\partial_{y}-y\partial_{x})$, $r(t)=\Omega(t)/\omega$ and
$\tilde{g}=8\pi a_{s}/d_{\perp}$. From this point, all equations would be
expressed in terms of these dimensionless space and time variables.

The only free parameters in this model are the interaction strength,
$\tilde{g}N$, the strength of the deformation, $\varepsilon$, and the function
$r(t)$ which specifies the rotation rate. We choose $\tilde{g}N=1000$ and
$\varepsilon=0.125$ in accordance with the experiment in [2]. The rotation
rate, $r(t)$ is ramped up from 0 to 1 in a time $\tau_{r}$. For $t<\tau_{r}$
we take $r(t)=\sin{\frac{\pi t}{2\tau_{r}}}$, while for longer times we take
$r(t)=1$.

Note that in the non-interacting limit ($\tilde{g}N=0$), if $\varepsilon\to
0$, the GPE can be written as,

$i\frac{d\psi}{dt}=\left((i\vec{\nabla}-q\vec{A})^{2}+\frac{(1-r(t)^{2})(x^{2}+y^{2})}{4}\right)\psi$
(5)

This is reminiscent of the Schrodinger equation for a particle with charge $q$
in the presence of an external vector potential $\vec{A}$ and a residual
harmonic trap [16, 17]. Here $q\vec{A}=r(t)(-y/2,x/2)$ is the effective vector
potential in the symmetric gauge. When $r(t)\to 1$, we recover the dynamics of
a charged particle purely in the presence of a magnetic field. As discussed in
Sec. I, this is the case when the centrifugal force vanishes and the coriolis
force mimics the Lorentz force. The eigenstates are Landau level wavefunctions
where the particles undergo cyclotron orbits. Each Landau level has a massive
degeneracy due to the different possible in-plane center of mass positions.

##  III Variational Wavefunction Ansatz

We use a time dependent variational wavefunction approach to approximate the
solution of Eq. (4). At all times we use a normalized gaussian ansatz for our
wavefunction,

$\begin{split}\psi=\frac{1}{\sqrt{\pi
l_{x}(t)l_{y}(t)}}\exp\left(-\frac{(x-\alpha(t)y)^{2}}{2l_{x}(t)^{2}}-\frac{y^{2}}{2l_{y}(t)^{2}}\right)\times\\\
\exp(i(\phi_{x}(t)x^{2}+\phi_{y}(t)y^{2}+\phi_{xy}(t)xy)).\end{split}$ (6)

The dynamics are encoded in six time dependent parameters: $l_{x}(t),l_{y}(t)$
represent the size of the condensate while $\phi_{x}(t),\phi_{y}(t)$ quantify
particle currents in the x and y directions respectively. The $\phi_{xy}(t)$
term represents a quadrupolar flow. In particular, it serves as a vortex-free
irrotational approximation to solid body rotation of the gas [5]. The
orientation of the principle axes of the cloud is determined by $\alpha$. This
form is convenient because during the dynamics, the principle axes of the
cloud are almost aligned with $x$ and $y$ axes. In that limit, $\alpha$
corresponds to the slight tilt from the equipotentials of the deformation term
in Eq. 4.

Gaussian ansatzes of this form have been used to obtain dynamical solutions of
the Gross-Pitaevskii equation [18, 19]. It should be a good model as long as
there are no dynamical instabilites to vortex formation [20]. We compared our
variational calculations to numerical integration of the GPE using finite
differences, and find good agreement. We note that, with a fixed grid, the
finite difference approach is limited in the length of time it can describe:
Eventually the cloud becomes larger than the simulation volume. Furthermore
the finite difference approach is numerically expensive compared to the
variational approach.

##  IV Variational Wavefunction dynamics

According to the standard formulation of the time dependent variational
principle [18], we calculate the action

$S=\int\\!dt\int\\!dx\,dy\,\left(i\psi^{*}\frac{d\psi}{dt}-\psi^{*}H(t)\psi\right)$
(7)

as a functional of the parameters
$l_{x}(t),l_{y}(t),\alpha(t),\phi_{x}(t),\phi_{y}(t),$ and $\phi_{xy}(t)$. We
take

$H(t)=-\nabla^{2}+\frac{x^{2}}{4}+\frac{y^{2}}{4}-\frac{\varepsilon
xy}{2}-r(t)L_{z}+\frac{1}{2}\tilde{g}N|\psi|^{2}.$ (8)

Note the factor of 2 in the interaction term compared to the right hand side
of Eq. (4). If we extremize Eq. 7 with respect to $\psi^{*}$, we recover the
Gross-Pittaevskii equation. As described in Appendix A, our variational
equations of motion are found by making Eq. (7) stationary with respect to the
variational parameters. We obtain six coupled ordinary differential equations
corresponding to the six parameters.

###  IV.1 Initial Conditions

For our initial conditions, we find the variational ground state in the
stationary trap. We set the rotation rate to zero, $r(0)=0$, and minimize

$E=\int dx\,dy\,\langle H(0)\rangle=\int dx\,dy\,\psi^{*}H(0)\psi$ (9)

under the constraint that the total number of particles is fixed. For the bulk
of our results, we use the experimentally relevant values $\tilde{g}N=1000$
and $\epsilon=0.125$, in which case $l_{x}(0)=5.02$, $l_{y}(0)=5.06$,
$\alpha(0)=0.125$, $\phi_{x}(0)=\phi_{y}(0)=\phi_{xy}(0)=0$. Note, the cloud
is almost cylindrically symmetric, with $l_{x}$ and $l_{y}$ differing by less
than a percent.

##  V Dynamics

We use the NDSolve package in Mathematica to numerically evolve the equations
of motion from Appendix A. This uses a robust general purpose algorithm with
dynamically adjusted step size. In the experiment, the anisotropy $\epsilon$
is time-independent, and in section V.1 we analyze that case. The anisotropic
term mixes Landau levels, and we can only drive the system into the
traditional Lowest Landau level by slowly turning off $\epsilon$ after ramping
up $r(t)$. Thus in section V.2, we add this ramp-down step.

###  V.1 Constant $\varepsilon$

Figure 1 shows the spatial evolution of the condensate density ($|\psi|^{2}$)
as a function of time for the experimental parameters. The condensate expands
along the y-direction while getting squeezed in the x-direction.

As illustrated in Fig. 2(a), the dynamics has two regimes: $t<\tau_{r}$, where
the rotation rate $r(t)$ is being swept, and $t>\tau_{r}$ where the rotation
rate is being held constant. The behavior in the latter regime is particularly
simple: $l_{x}$ falls exponentially towards 1, $l_{y}$ grows exponentially,
and $\phi_{x},\phi_{y},\phi_{xy},$ and $\alpha$ all approach constants.

Figure 1: Condensate density ($\rho=|\psi|^{2}$) as a function of time, for a
Bose gas in a 2D harmonic potential of frequencty $\omega$ in the presence of
a rotating quadrupolar deformation of fixed amplitude $\varepsilon$, as
defined in Eq. (3). From left to right, the rotation rate is ramped up,
starting at $\Omega=0$ at time $t=0$ and ending at $\Omega=\omega$ at time
$\tau_{r}$. The rotation frequency is then held fixed for longer times. Here
$\varepsilon=0.125$. Brighter colors represent higher normalized density
$\rho/\rho_{\rm max}$.  Figure 2: Model and Wavefunction ansatz parameters as
a function for time for constant $\varepsilon$. (a) $r(t)$ (blue) and
$\varepsilon=0.125$ (red) (b) $l_{x}(t)$ (blue) saturates at 1 while
$l_{y}(t)$ (red) grows exponentially at long times (c) $\phi_{xy}(t)$ (blue)
saturates at 0.5 and $\alpha(t)$ (red) saturates at $\varepsilon/8$ (d)
$\phi_{y}(t)$ (blue) and $\phi_{x}(t)$ (red) saturate at $\varepsilon/16$ at
long times.

This behavior can be understood in terms of an analog of the Lorentz
$\vec{E}\times\vec{B}$ drift of charged particles in electric and magnetic
fields. In this analogy, $q\vec{E}=\vec{\nabla}(\varepsilon xy/2)$, and
$\vec{B}=r(t)\hat{z}=\hat{z}$ at long times. The $\vec{E}\times\vec{B}$ drift
velocities in the x and y direction are respectively given by,
$v_{x}\sim-\varepsilon x/2$ and $v_{y}\sim\varepsilon y/2$. These equations
are consistent with the exponential squeezing along the x-axis and exponential
expansion along the y axis. This classical argument misses the fact that
$l_{x}$ saturates at 1. This feature comes from wave mechanics: The decreasing
x-width costs kinetic energy due to the $1/(2l_{x}^{2})$ term and thus there
is a quantum pressure which halts the squeezing. An increasing y-width doesn’t
cost any energy in the rotating frame due to the cancellation of the external
harmonic trap with the centrifugal force. At long times, as
$l_{y}\rightarrow\infty$, the kinetic energy, $1/(2l_{y}^{2})$ and interaction
energy $\tilde{g}N/4\pi l_{x}l_{y}$ tend to zero.

At long times, $\alpha(t)\rightarrow\varepsilon/8$ denoting a finite tilt of
the condensate away from the y axis. This wavefunction tilt can be understood
from a slightly more elaborate analysis of the classical equations of motion.
For this model, the Lorentz force equation is linear,
$md\vec{v}/dt=q\vec{E}+q\vec{v}\times\vec{B}$, and Fourier transforming it in
the time domain yields

$-\nu^{2}\begin{bmatrix}\bar{x}(\nu)\\\
\bar{y}(\nu)\end{bmatrix}=\begin{bmatrix}0&\varepsilon+2i\nu\\\
\varepsilon-2i\nu&0\end{bmatrix}\begin{bmatrix}\bar{x}(\nu)\\\
\bar{y}(\nu)\end{bmatrix}$ (10)

with $\bar{x}(\nu)=\int dt\,e^{i\nu t}x(t)$, and a similar expression for
$\bar{y}$. From Eq. (10), we find a quartic equation for the allowed
frequencies of free oscillation, $\nu$. As $\epsilon\to 0$, the solutions are
$\nu=\pm i\varepsilon/2,\pm 2$. The first solutions are purely imaginary and
represent the squeezing already discussed. The latter describe cylclotron
motion. The eigenvector of the exponentially growing solution has
$\bar{y}/\bar{x}=\varepsilon/8$, corresponding to the observed tilt. This is
analogous to inertial drift in plasma physics [21].

The modulus and phase of the condensate wavefunction can be interpreted in
terms of fluid mechanics, $\psi=\sqrt{\rho}e^{i\Phi}$, where $\rho$ is the
density and ${\vec{v}}=(h/m)\vec{\nabla}\Phi$. Hence, in dimensionless units,
our ansatz corresponds to a flow-pattern $v_{x}=2(2\phi_{x}x+\phi_{xy}y)$ and
$v_{y}=2(2\phi_{y}y+\phi_{xy}x)$. When $\phi_{xy}=0.5$, the angular momentum
of the condensate resembles that of solid body rotation of the condensate (see
Appendix B). In Fig. 2(d), one sees that the particle current terms $\phi_{y}$
and $\phi_{x}$ saturate at $\varepsilon/16$ at long times.

Further insight into the structure of the wavefunction at long times comes
from transforming into a coordinate system which is aligned with the principle
axes of the condensate. We define $\tilde{x}=(x-\alpha y)$ and
$\tilde{y}=(y+\alpha x)$, and take the parameters to have their long-time
values, $\phi_{x}=\phi_{y}=\varepsilon/16$, $\alpha=\varepsilon/8$, and
$\phi_{xy}=1/2$, $l_{x}=1$. Up to linear order in $\varepsilon$, the
variational wavefunction is then

$\psi\sim\frac{1}{\sqrt{\pi\ell_{y}}}\exp\left(-\frac{\tilde{x}^{2}}{2}+i\frac{\varepsilon\tilde{y}^{2}}{8}+i\frac{\tilde{x}\tilde{y}}{2}-\frac{\tilde{y}^{2}}{2l_{y}^{2}}\right)$
(11)

At long times $l_{y}\to\infty$, and one can ignore the last term in the
exponential. The three remaining terms have physical meanings. The first
limits the spread of the particles in the $\tilde{x}$ axis. The second
represents out-going currents in the $\tilde{y}$ axis, corresponding to the
expansion. The third is an irrotational flow pattern that approximates solid
body rotation. If $\epsilon=0$ and $l_{y}\to\infty$, this is in the lowest
Landau level, in the symmetric gauge. Transforming to the Landau gauge
corresponds to multiplying by $\exp(-i\tilde{x}\tilde{y}/2)$, which eliminates
the third term.

In our ’squeezing’ solution of the classical Lorentz force equations earlier
in the section, the rate of exponential expansion was $\varepsilon/2$. This is
consistent with $\phi_{y}=\varepsilon/16$ at long times since the outgoing
velocity in the $\tilde{y}$ direction becomes $\varepsilon\tilde{y}/2$ by
calculating gradient of the phase in Eq. (11). Due to a finite wavefunction
tilt and exponential expansion, we expect $\phi_{x}=\phi_{y}$ at long times.

Figure 3: Logarithm of the average traditional Landau level occupation,
ln$(N_{L}+1)$ (blue) and logarithm of the average distorted Landau level
occupation, ln$(N_{L}^{\prime}+1)$ (red), as a function of time for the ramp
described in Fig. 2(a). At long times, $N_{L}$ grows exponentially while
$N_{L}^{\prime}$ decreases and tends to 0.

The saturation of $l_{x}$ is consistent with being in the Lowest Landau level.
There are some subtleties however, as the eigenstates of the deformed trap are
somewhat different than those of traditional Landau levels. We can quantify
this feature by calculating the average number of landau levels which are
occupied,

$N_{L}=\int dx\,dy\,\psi^{*}\frac{(H_{L}-1)}{2}\psi$ (12)

where

$H_{L}=-\nabla^{2}+\frac{x^{2}}{4}+\frac{y^{2}}{4}-i(x\partial_{y}-y\partial_{x})$
(13)

is the symmetric gauge Hamiltonian of a charged particle in a uniform magnetic
field: $N_{L}=0$ when the wavefunction is in the lowest Landau level. We plot
ln$(N_{L}+1)$ as a function of time in Fig. 3. As the rotation rate increases,
$N_{L}$ drops. However, after the ramp is concluded
($t\varepsilon\omega/2=4$), the average Landau level number grows
exponentially. This behavior is evident from Eq. (11), for which
$N_{L}\sim\varepsilon^{2}l_{y}^{2}/16$. Since $l_{y}$ grows exponentially, so
does $N_{L}$.

Following the supplementary of [2], one can alternatively use a gauge
transformation to define distorted Landau levels which take into account the
elliptic deformation of the trap. The gauge transformation is given by
$\psi^{\prime}=G\psi$ where $G=\exp(i\varepsilon(x^{2}-y^{2})/8)$. One can
then define the distorted Landau levels and count their average occupation via

$N_{L}^{\prime}=\int
dx\,dy\,\psi^{\prime*}\frac{(H_{L}^{\prime}-1)}{2}\psi^{\prime}$ (14)

where

$H_{L}^{\prime}=-\frac{\nabla^{2}}{\gamma}+\gamma\frac{(x^{2}+y^{2})}{4}-i(x\partial_{y}-y\partial_{x}),$
(15)

with $\gamma=\sqrt{1+\varepsilon^{2}/4}$. We plot ln$(N_{L}^{\prime}+1)$ as a
function of time in Fig. 3. In terms of these distorted Landau levels,
$N_{L}^{\prime}$ drops as the rotation rate increases and tends to 0 at long
times. The condensate is thus driven into this distorted lowest Landau level.

When $l_{y}\to\infty$, the difference between Eq. (11) and an undistorted
lowest Landau level wavefunction is the phase factor
$\exp(i\varepsilon\tilde{y}^{2}/8)$, which, as already explained, quantifies
the outward expansion of the condensate. This expansion is driven by the
Lorentz $E\times B$ drift, which occurs as long as $\varepsilon\neq 0$. Thus
if we wish to end up in the undistorted lowest Landau level, we must ramp
$\varepsilon$ to zero. In the next section we include such a ramp, and analyze
the resulting dynamics.

###  V.2 Varying $\varepsilon$

We adiabatically turn off the anisotropy $\varepsilon$ after the rotation rate
$r(t)$ has been ramped up to 1. For $t<\tau_{r}$,
$\varepsilon(t)=\varepsilon_{0}=0.125$. For $\tau_{r}<t<\tau_{\varepsilon}$,
$\varepsilon(t)=\varepsilon_{0}\frac{\tau_{r}+\tau_{\varepsilon}-t}{\tau_{\varepsilon}}$.
It is held at zero for all later times. Thus $\varepsilon(t)$ reaches zero at
time $t=\tau_{r}+\tau_{\varepsilon}$.

For adiabaticity, we need $\tau_{\varepsilon}>1/\omega$. In order to enter the
lowest Landau level, we want the squeezing dynamics to continue until
$l_{y}\gg 1$ and $l_{x}\sim 1$. Thus $\tau_{\varepsilon}$ needs to be
moderately larger than $1/(\varepsilon_{0}\omega)$. In Fig. 4 and 5, we
consider the case where $\tau_{\epsilon}=8/(\varepsilon_{0}\omega)$. We only
show data up to $t=\tau_{r}+\tau_{\varepsilon}$ as all parameters are time
independent once $\varepsilon=0$. Remarkably, as shown in Fig. 4, this ramp
results in a final $N_{L}\approx 0.03$, which means that the probability of a
particle being in the lowest Landau level exceeds 97%.

For these parameters, $N_{L}$ in Fig. 4 is a non-monotonic function of time.
The exponential growth, seen in Fig. 3 competes with the expected drop due to
reducing $\varepsilon$. As described below, this competition also results in a
non-monotonic dependence of the final value of $N_{L}$ on $\tau_{\epsilon}$.

The system and wavefunction parameters as a function of time are shown in Fig.
5. Similar to the constant $\varepsilon$ case in Fig. 2(b), $l_{x}$ falls
exponentially towards 1, and $l_{y}$ grows. For moderate times, the growth is
roughly exponential, but it flattens out and saturates as
$t\to\tau_{r}+\tau_{\varepsilon}$. As before, $\phi_{xy}\to 0.5$. The
wavefunction tilt, parameterized by $\alpha$, switches directions during the
$\varepsilon$ ramp. Unlike the findings in Fig. 2(c), the final tilt depends
on details of the ramp.

The particle current terms, $\phi_{x}$ and $\phi_{y}$ take values $-\alpha/2$
and $\alpha/2$ respectively. Transforming into the tilted coordinates again
(treating $\alpha$ as small) gives the wavefunction in Eq. (11), but this time
with $\varepsilon=0$. This is a lowest Landau level wavefunction. There is,
however, a correction term of higher order in $\alpha$ in the wavefunction. It
is given by $\exp(i\alpha^{3}(\tilde{x}^{2}-\tilde{y}^{2})/2)$. Due to this
correction, the leading order correction to $N_{L}$ from Eq. (12) is given by
$\alpha^{4}l_{y}^{2}$. Thus as long as $\alpha^{4}l_{y}^{2}\ll 1$, $N_{L}\sim
0$. Since $\alpha\sim\varepsilon_{0}/8$ and $l_{y}\sim
e^{\varepsilon_{0}\omega t/2}$, this condition limits the length of the ramp
down period of $\varepsilon$ to about
$8/(\varepsilon_{0}\omega)\ln(8/\varepsilon_{0})$. When $\varepsilon$ is
switched off on this timescale, the condensate is left in the undistorted
lowest Landau level at the end.

Figure 4: Average Landau level occupation, $N_{L}$, as a function of time,
when the trap asymmetry is ramped to zero. Parameters are shown in Fig. 5.
$N_{L}=0.03$ at the end, implying particles having high probability of being
in the lowest Landau level. Figure 5: Model and Wavefunction ansatz parameters
as a function of time in the protocol where $\varepsilon$ is ramped to zero.
(a) $r(t)$ (blue) and $\varepsilon$ (red) (b) $l_{x}(t)$ (blue) and
ln$(l_{y}(t))$ (red); ramping $\varepsilon\to 0$ cuts off the exponential
growth of $l_{y}$ that was seen in Fig. 2 (c) $\phi_{xy}(t)$ (blue) again
saturates at 0.5 and $\alpha(t)$ (red) saturates at a negative value which
depends on the ramp. (d) $\phi_{y}(t)$ (blue) and $\phi_{x}(t)$ (red) saturate
at $\pm\alpha/2$. All parameters are subsequently time independent.

##  VI Summary and Outlook

We modeled the dynamics of a condensed Bose gas in a rotating anisotropic
trap. We used a time dependent variational wavefunction to calculate the time
evolution of the system. We find that the behavior at long times is well
described by the classical equations of motion of a charged particle in
crossed electric and magnetic fields. We calculate the average Landau level
occupation as a function of time. We find that when the trap anisotropy is
constant in time, the condensate gets squeezed into a long thin shape which
occupies a deformed lowest Landau level. If the trap anisotropy is switched
off at the appropriate time, the condensate enters the traditional lowest
Landau level.

Confinement to the Lowest Landau level is a necessary prerequisite for
studying analogs of quantum Hall physics with cold atoms. We emphasize,
however, that the states produced by this protocol are rather different from
the homogeneous liquids which lead to the quantum Hall effect. Nonetheless,
this is an important first step. Moreover there are a number of interesting
purely interactions-driven phenomena which can be explored starting from these
highly anisotropic lowest Landau level states. For example, the longer-time
dynamics are described by correlated hopping models with unique sets of
conserved quantities that lead to unusual transport [22]. Variants of this
protocol can also lead to interesting pattern formation [23].

##  VII Acknowledgements

We thank Martin Zwierlein and his research group for insightful discussions.
This work was supported by the NSF Grant No. PHY- 2110250.

##  Appendix A Variational equations of motion

Using time dependent variational principle, we first calculate the action
given by,

$\displaystyle S$ $\displaystyle=$
$\displaystyle\int\\!\\!dt\\!\int\\!\\!dx\,dy\,\left(i\psi^{*}\frac{d\psi}{dt}-\psi^{*}H(t)\psi\right)$
$\displaystyle=$
$\displaystyle\int\\!\\!dt\,\mathcal{F}(l_{x}(t),l_{y}(t),\alpha(t),\phi_{x}(t),\phi_{y}(t),\phi_{xy}(t))$

where the variational wavefunction ansatz $\psi$ is given by Eq. (6),

$\begin{split}\psi=\frac{1}{\sqrt{\pi
l_{x}(t)l_{y}(t)}}\exp\left(-\frac{(x-\alpha(t)y)^{2}}{2l_{x}(t)^{2}}-\frac{y^{2}}{2l_{y}(t)^{2}}\right)\times\\\
\exp(i(\phi_{x}(t)x^{2}+\phi_{y}(t)y^{2}+\phi_{xy}(t)xy))\end{split}$ (17)

The spatial integrals in the action are gaussian and can be done analytically,
yielding $\mathcal{F}=\mathcal{F}_{t}+\mathcal{F}_{p}+\mathcal{F}_{\rm
trap}+\mathcal{F}_{\rm rot}+\mathcal{F}_{\rm int}$,

$\displaystyle\mathcal{F}_{t}=$
$\displaystyle-\frac{d\phi_{x}}{dt}\left(\frac{l_{x}^{2}}{2}+\frac{\alpha^{2}l_{y}^{2}}{2}\right)-\frac{d\phi_{y}}{dt}\frac{l_{y}^{2}}{2}-\frac{d\phi_{xy}}{dt}\frac{\alpha
l_{y}^{2}}{2}$ $\displaystyle\mathcal{F}_{p}=$
$\displaystyle\frac{1+\alpha^{2}}{2l_{x}^{2}}+\frac{1}{2l_{y}^{2}}+l_{x}^{2}\left(2\phi_{x}^{2}+\frac{\phi_{xy}^{2}}{2}\right)$
$\displaystyle+l_{y}^{2}$
$\displaystyle\left(2\alpha^{2}\phi_{x}^{2}+2\phi_{y}^{2}+2\alpha\phi_{xy}(\phi_{x}+\phi_{y})+\phi_{xy}^{2}\frac{1+\alpha^{2}}{2}\right)$
$\displaystyle\mathcal{F}_{\rm rot}=$ $\displaystyle
r\left(\alpha(\phi_{x}+\phi_{y})l_{y}^{2}+\phi_{xy}\frac{l_{x}^{2}+(\alpha^{2}-1)l_{y}^{2}}{2}\right)$
$\displaystyle\mathcal{F}_{\rm trap}=$ $\displaystyle-\varepsilon\alpha
l_{y}^{2}/4+\frac{l_{x}^{2}}{8}+\frac{\alpha^{2}l_{y}^{2}}{8}+\frac{l_{y}^{2}}{8}$
$\displaystyle\mathcal{F}_{\rm int}=$ $\displaystyle\frac{\tilde{g}N}{2\pi
l_{x}l_{y}}$

We extremize the action using the Euler-Lagrange equations,

$\frac{\partial\mathcal{F}}{\partial\Phi}-\frac{d}{dt}\frac{\partial\mathcal{F}}{\partial\dot{\Phi}}=0$
(18)

where $\Phi=\\{l_{x},l_{y},\alpha,\phi_{x},\phi_{y},\phi_{xy}\\}$ and
$\dot{\Phi}$ denotes a time derivative. We thus get six coupled nonlinear
ordinary differential equations. We numerically solve these equations.

##  Appendix B Parameter $\phi_{xy}$ approximating solid body rotation

The term $\exp(i\phi_{xy}xy)$ in Eq. (17) represents an irrotational
quadrupolar flow. This term contributes to a non-zero angular momentum,
$L_{z}^{c}$ given by,

$L_{z}^{c}=\langle\psi|L_{z}|\psi\rangle=\frac{\phi_{xy}(l_{y}^{2}-l_{x}^{2})}{2}$
(19)

If the condensate were rotating as a solid body with angular velocity
$\Omega$, condensate velocity profile is $\vec{v}=\vec{r}\times\vec{\Omega}$
and the solid body angular momentum $L_{z}^{sb}$ is $mr^{2}\Omega$. In
dimensionless units, $L_{z}^{sb}$ is,

$L_{z}^{sb}=\langle\psi|\frac{r^{2}}{2}|\psi\rangle=\frac{(l_{y}^{2}+l_{x}^{2})}{4}$
(20)

The condensate approximates solid body rotation when $l_{y}\gg l_{x}$ and
$\phi_{xy}=0.5$, making $L_{z}^{c}/L_{z}^{sb}\to 1$.

## References

  * Cooper [2008] N. Cooper, Rapidly rotating atomic gases, Advances in Physics 57, 539 (2008).
  * Fletcher  _et al._ [2021] R. J. Fletcher, A. Shaffer, C. C. Wilson, P. B. Patel, Z. Yan, V. Crépel, B. Mukherjee, and M. W. Zwierlein, Geometric squeezing into the lowest landau level, Science 372, 1318 (2021).
  * Dalibard  _et al._ [2011] J. Dalibard, F. Gerbier, G. Juzeliūnas, and P. Öhberg, Colloquium: Artificial gauge potentials for neutral atoms, Rev. Mod. Phys. 83, 1523 (2011).
  * Hauke  _et al._ [2012] P. Hauke, O. Tieleman, A. Celi, C. Ölschläger, J. Simonet, J. Struck, M. Weinberg, P. Windpassinger, K. Sengstock, M. Lewenstein, and A. Eckardt, Non-abelian gauge fields and topological insulators in shaken optical lattices, Phys. Rev. Lett. 109, 145301 (2012).
  * Fetter [2009] A. L. Fetter, Rotating trapped bose-einstein condensates, Rev. Mod. Phys. 81, 647 (2009).
  * Abo-Shaeer  _et al._ [2001] J. R. Abo-Shaeer, C. Raman, J. M. Vogels, and W. Ketterle, Observation of vortex lattices in bose-einstein condensates, Science 292, 476 (2001).
  * Haljan  _et al._ [2001] P. C. Haljan, I. Coddington, P. Engels, and E. A. Cornell, Driving bose-einstein-condensate vorticity with a rotating normal cloud, Phys. Rev. Lett. 87, 210403 (2001).
  * Engels  _et al._ [2003] P. Engels, I. Coddington, P. C. Haljan, V. Schweikhard, and E. A. Cornell, Observation of long-lived vortex aggregates in rapidly rotating bose-einstein condensates, Phys. Rev. Lett. 90, 170405 (2003).
  * Hess and Fairbank [1967] G. B. Hess and W. M. Fairbank, Measurements of angular momentum in superfluid helium, Phys. Rev. Lett. 19, 216 (1967).
  * Sinha and Shlyapnikov [2005] S. Sinha and G. V. Shlyapnikov, Two-dimensional bose-einstein condensate under extreme rotation, Phys. Rev. Lett. 94, 150401 (2005).
  * Sonin [2005] E. B. Sonin, Ground state and tkachenko modes of a rapidly rotating bose-einstein condensate in the lowest-landau-level state, Phys. Rev. A 72, 021606 (2005).
  * Schweikhard  _et al._ [2004] V. Schweikhard, I. Coddington, P. Engels, V. P. Mogendorff, and E. A. Cornell, Rapidly rotating bose-einstein condensates in and near the lowest landau level, Phys. Rev. Lett. 92, 040404 (2004).
  * Gemelke  _et al._ [2010] N. Gemelke, E. Sarajlic, and S. Chu, Rotating few-body atomic systems in the fractional quantum hall regime (2010), arXiv:1007.2677 [cond-mat.quant-gas] .
  * Madison  _et al._ [2000] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard, Vortex formation in a stirred bose-einstein condensate, Phys. Rev. Lett. 84, 806 (2000).
  * Bretin  _et al._ [2004] V. Bretin, S. Stock, Y. Seurin, and J. Dalibard, Fast rotation of a bose-einstein condensate, Phys. Rev. Lett. 92, 050403 (2004).
  * Darwin [1931] C. G. Darwin, The diamagnetism of the free electron, Mathematical Proceedings of the Cambridge Philosophical Society 27, 86–90 (1931).
  * Maksym and Chakraborty [1990] P. A. Maksym and T. Chakraborty, Quantum dots in a magnetic field: Role of electron-electron interactions, Phys. Rev. Lett. 65, 108 (1990).
  * Pérez-García  _et al._ [1997] V. M. Pérez-García, H. Michinel, J. I. Cirac, M. Lewenstein, and P. Zoller, Dynamics of bose-einstein condensates: Variational solutions of the gross-pitaevskii equations, Phys. Rev. A 56, 1424 (1997).
  * Castin [2001] Y. Castin, Bose-einstein condensates in atomic gases: Simple theoretical results, in  _Coherent atomic matter waves_ , edited by R. Kaiser, C. Westbrook, and F. David (Springer Berlin Heidelberg, Berlin, Heidelberg, 2001) pp. 1–136.
  * Recati  _et al._ [2001] A. Recati, F. Zambelli, and S. Stringari, Overcritical rotation of a trapped bose-einstein condensate, Phys. Rev. Lett. 86, 377 (2001).
  * Carlqvist [1984] P. Carlqvist, Plasma drift towards a plane equipotential surface (1984).
  * Fremling  _et al._ [2018] M. Fremling, C. Repellin, J.-M. Stéphan, N. Moran, J. K. Slingerland, and M. Haque, Dynamics and level statistics of interacting fermions in the lowest landau level, New J. Phys 20, 103036 (2018).
  * Mukherjee  _et al._ [2021] B. Mukherjee, A. Shaffer, P. B. Patel, Z. Yan, C. C. Wilson, V. Crépel, R. J. Fletcher, and M. Zwierlein, Crystallization of bosonic quantum hall states (2021), arXiv:2106.11300 [cond-mat.quant-gas] .