# DLDNN: Deterministic Lateral Displacement Design Automation by Neural
Networks

Farzad Vatandoust
Department of Mechanical Engineering
Iran University of Science and Technology
<EMAIL_ADDRESS>
&Hoseyn A. Amiri
Department of Biomechanics
Iran University of Science and Technology
<EMAIL_ADDRESS>
&Sima Mas-hafi
Department of Mechanical Engineering
Kharazmi University
<EMAIL_ADDRESS>
###### Abstract

Size-based separation of bioparticles/cells is crucial to a variety of
biomedical processing steps for applications such as exosomes and DNA
isolation. Design and improvement of such microfluidic devices is a challenge
to best answer the demand for producing homogeneous end-result for study and
use. Deterministic lateral displacement (DLD) exploits a similar principle
that has drawn extensive attention over years. However, the lack of predictive
understanding of the particle trajectory and its induced mode makes designing
a DLD device an iterative procedure. Therefore, this paper investigates a fast
versatile design automation platform to address this issue. To do so,
convolutional and artificial neural networks were employed to learn velocity
fields and critical diameters of a wide range of DLD configurations. Later,
these networks were combined with a multi-objective evolutionary algorithm to
construct the automation tool. After ensuring the accuracy of the neural
networks, the developed tool was tested for 12 critical conditions. Reaching
the imposed conditions, the automation components performed reliably with
errors of less than 4%. Moreover, this tool is generalizable to other field-
based problems and since the neural network is an integral part of this
method, it enables transfer learning for similar physics. All the codes
generated and used in this study alongside the pre-trained neural network
models are available on https://github.com/HoseynAAmiri/DLDNN.

_K_ eywords Microfluidics $\cdot$ Deterministic Lateral Displacement $\cdot$
Convolutional Neural Network $\cdot$ Design Automation $\cdot$ Multi-objective
Optimization $\cdot$ Particle Separation

##  1 Introduction

The need for biomolecular analysis was the primary incentive for microfluidics
development, In recent years, cell separation studies have drawn significant
attention to microfluidic devices [1] which could be categorized based on
particle properties such as size and other attributes [2]. Deterministic
lateral displacement (DLD), a size-based particle/cell separation method, is a
crucial technique in the fields of medical science and biology for rapid
diagnosis and testing. Its low cost and easy operation has created a growing
interest in developing DLD devices.

In 2004, the first DLD device sample was proposed by Huang et al. [3]. They
used bifurcation of laminar flow around arrays of obstacles to induce
different lateral displacements to the particles. Interestingly, it was shown
that particles larger than a certain critical diameter ($D_{c}$) follow a
continuous lateral displacement path (bumped mode) while smaller particles
continue to move along their entry point towards the end of the channel
(zigzag mode). Later in 2006, Inglis et al. [4] derived an analytical
approximation of the parameter $D_{c}$ to facilitate the design of DLD and
Davis et al. [5] modified the formula for improving its accuracy.

Further advancement of DLD devices aimed at increasing the throughput,
enhancing the separation quality, and reducing clogging. Zeming et al. [6] and
Kim et al. [7] investigated the effect of the gap size parameter on separation
quality and throughput. In addition to array arrangements, researchers also
employed different pillar shapes. Wang et al. [8] evaluated the influence of
triangular pillar shapes while Zeming et al. [9, 10] discussed various shapes
like rectangular, L-shaped, and I-shaped. Furthermore, Hyun et al. [11]
performed topology optimization on the shape of the pillars in the DLD device,
and their method allowed to increase the gap between particles, thus reducing
the clogging. Finally, Dincau et al. [12] tested the performance of DLD
devices in various Reynolds numbers ($10<Re<60$). The results suggested that
the $D_{c}$ decreases with an increase in Reynolds number due to the changes
in the flow field and wakes behind pillars. Therefore, the performance of the
DLD device can be tuned by adjusting the flow rate parameter. Overall, these
works thoroughly established the direct impact of the geometrical parameters
and Reynolds number on $D_{c}$, throughput, and separation quality. Various
effective parameters make the engineering of a DLD device highly iterative and
resource-intensive; hence, a design automation platform is required to
alleviate this problem.

The automation platforms typically include a surrogate model of the actual
physics for fast prediction and an algorithm for finding the best device setup
[13, 14]. Optimization algorithms are well-established in various problems for
myriad applications [15, 16, 17, 18, 19]. Thus, the main challenge of
automation is constructing a reliable surrogate model. Recently, machine
learning methods have displayed great potential in detecting patterns and
learning PDE equations [20, 21, 22, 23] as well as extracting or
reconstructing flow field features [24]. Jin et al. [25] utilized a fusion
convolutional neural network (CNN) to predict the velocity field around a
cylinder by setting the pressure around the cylinder as an input. Lee et al.
[26] used CNN and generative adversarial neural networks (GANs) with and
without imposing conservation of mass and momentum to generate a flow field
around a cylinder. The conclusions state that all networks are able to predict
flow in near future. However, in the aspect of long-term generalization and
prediction, the network with imposed physical law have better performance. In
the field of the physics-informed neural network (PINN), Raissi et al. [27,
28, 29] have dedicated many studies to encoding the governing equations into a
neural network. In 2019, they introduced the concept of PINN and applied
Navier-Stokes equations to a standard neural network to predict flow over a
cylinder [27]. More information about these types of neural networks can be
found in Raissi’s review paper [29]. Saker et al. [30] combined CNN and ANN to
generate flow over the airfoil. They used CNN to extract geometrical features
from the image of an airfoil and then fed it into an ANN to predict the flow
past a cylinder. In addition to the flow field, there are other studies
involving the prediction of fields of any kind such as heat flux [31]. As seen
from the literature, machine learning methods are capable of imitating the
behavior of underlying fluid flow physics. However, to the best of our
knowledge, there are not many studies investigating the DLD flow field
prediction, and the use of machine learning was limited to the area of
particle detection by CNN to facilitate the experimental process [32].

Therefore, this study explores the application of deep learning methods for
surrogate modeling of DLD devices and combines them with a multi-objective
algorithm creating a design automation platform. To do so, circular pillars
were chosen as a base for the design, thereafter the data-set consisting of x
and y components of velocity field was generated by numerical simulation for a
wide range of DLD device conditions. Followed by that, a CNN was developed and
trained to reconstruct these fields around the pillars. The network
demonstrated good accuracy in predicting the field and later was used as a
data augmentation tool for training a fully-connected neural network (FCNN) to
directly predict $D_{c}$. These components alongside non-dominant sorting
genetic algorithm type 3 (NSGA3) constructed the design automation platform.
The developed tool provides many features from achieving the particular design
parameters for a user-specified desired performance, predicting the device
working range after fabrication, and modeling the particle behaviors for any
given device length. Overall, this approach can facilitate the design process
while saving time and resources by avoiding repetitive trial and error
processes to reach the desired design.

##  2 Methodology

###  2.1 Data preparation

A DLD array of circular posts was built in the CFD solver to generate a
comprehensive data-set for neural network training. To do so, the
incompressible Navier-Stokes and continuity equations were solved since the
Reynolds number lies within the laminar flow regime. Figure 1a illustrates a
general schematic of the DLD design and the role of critical diameter in the
separation of green ($D_{g}>D_{c}$) and red ($D_{r}<D_{c}$) particles. The
dashed segment demonstrates the excerpt used for solving the periodic flow
model which holds all the geometrical parameters that influence the critical
diameter. These geometrical parameters consist of R, N, $G_{x}$, and $G_{y}$.
Where, R is pillars radius, N is the pillars arrays period number
($N=1/\alpha$) defining their periodic repetition, and $G=G_{x}=G_{y}$ are the
identical horizontal and vertical gaps between pillars. Then, in order to feed
the velocity fields to the proposed CNN, the decomposed fields in the
Cartesian coordinate system were mapped to squared fields of 128$\times$128
structured data as drawn in Figure 1b. To achieve this mapping, the domain was
normalized by $L=2R+G$, followed by a shear mapping in the y direction as,

$\displaystyle\begin{pmatrix}x^{\prime}\\\
y^{\prime}\end{pmatrix}=\begin{pmatrix}1&0\\\
-\frac{1}{N}&1\end{pmatrix}\begin{pmatrix}\frac{1}{L}&0\\\
0&\frac{1}{L}\end{pmatrix}\begin{pmatrix}x\\\ y\end{pmatrix}.$ (1)

###  2.2 Field generative CNN

Field-generative neural networks provide a great alternative for CFD solvers
with lower computational costs and the capacity of transfer learning for
similar applications. Here, a convolutional network was used for
reconstructing velocity fields as illustrated in Figure 1c. The neural network
architecture used in this work, inspired by Dosovitskiy et al. [33], generates
data as a great asset in DLD cases since the existing methods can be time-
consuming and resource intensive. To capture a wide range of operational
conditions, three inputs were defined namely, Re, N, and f where,

$f=\frac{2R}{2R+G}.$ (2)

It is important to take the ratio of G and R instead of their values because
the latter does not necessarily represent a unique shape. The inputs were fed
into a separate set of dense layers for higher dimensional representation and
after concatenation, they went through two sequential dense layers. Next, they
were connected to four upsampling layers combined with convolutional. It
should be noted that since both the x and y components of the velocity field
are required, the mentioned sub-structure is used in parallel in order to
reconstruct both fields.

Figure 1: a) DLD device principle, b) data preparations, and c) CNN
architecture.

###  2.3 Particle tracing and critical diameter extraction

Particle tracing and critical diameter extraction are crucial parts of the
performance evaluation of DLD devices. Consequently, particles’ behavior in
the flow field and their interaction with obstacles were simulated based on
the following models. First, the massless particle tracing scheme was
utilized, because the particles and the flow have negligible inertia. Then,
particles positions were approximated from velocity fields with Runge-Kutta
method. Finally, for simulating the particle-pillar interaction, the wall
distance function was numerically calculated to measure the minimum spatial
distance from pillars. In case of contact (when the distance between particle
and pillar is equal to the particle radius), the reflect function modifies the
particle velocity ($V_{p}$) as,

$\vec{n}\cdot\vec{V}_{p,reflected}=-\vec{n}\cdot\vec{V}_{p},$ (3)
$\vec{t}\cdot\vec{V}_{p,reflected}=\vec{t}\cdot\vec{V}_{p},$ (4)

where the $\vec{n}$ and $\vec{t}$ are the unit vectors in the normal and
tangential directions of particle and pillar contact point, respectively (see
Appendix A). Once the particle’s accurate pathlines are obtained, the critical
diameter was extracted via the Newtonian root finding method. In this
approach, first, two particles with the highest and lowest diameters possible
are simulated. The particle with the bumped mode receives a positive value and
the particle with the zigzag mode is assigned a negative value. Thus, the
Newtonian method can narrow the diameter interval to detect the root which is
the critical diameter (see Appendix B).

###  2.4 Design automation

The design automation algorithms convert the user-defined specification to
geometry properties and flow rate to readily meet the desired criteria. To
boost the function evaluation speed, one can utilize neural networks as the
surrogate model instead of the whole physics/experimentation. This alternative
facilitates the design process of a DLD device to reach the best possible
configuration, especially when combined with an evolutionary algorithm.

As presented in Figure 2, design automation inputs are the diameters of
particles to be separated ($D_{1}$ and $D_{2}$), desired constraints to be
applied on f, N, Re ($C_{f}$, $C_{N}$, and $C_{Re}$), and finally the trade-
off between flexibility and stability of the design ($\phi$). Flexibility
indicates the ability of the proposed design to cover a wide range of critical
diameters while stability specifies its immunity to change due to minor flow
rate fluctuations (see Appendix C). These two terms were included to consider
the possible device behavior after fabrication since it can be altered by
adjusting the Reynolds number.

After defining the inputs, a multi-objective optimization algorithm coupled
with a pre-trained FCNN tunes the f, N, Re, and G to reach the optimized
design. In this design automation tool, an NSGA3 with five reference
directions and population size of 260 was chosen and implemented in Pymoo
library [34] in Python 3. The pre-trained FCNN was used for direct critical
diameter prediction (direct neural network) since extracting critical diameter
from CNN outputs slows down the optimization process. Finally, the tuned
parameters in addition to bandwidth (BW) are extracted and presented in the
outputs. BW is the difference between the maximum and the minimum critical
diameter for that design serving as an index of flexibility. Note that the
parameter G is used in optimization to nondimensionalize the critical diameter
matching the developed neural network’s $D_{c}$.

Furthermore, the pre-trained CNN was attached to the end of the automation
process to provide the ability to simulate the particle trajectory in the
optimum design for a given number of periods. This feature helps visualize the
particle behavior at the full length of the device.

Figure 2: The scheme of design automation platform.

##  3 Results and discussion

###  3.1 CNN training and data generation

The CNN was developed for fast velocity field generation which accelerates
particle tracing. In this section, the network architecture, performance, and
utility for data augmentation were analyzed. For this purpose, a sufficient
amount of data was selected to cover a wide range of DLD device
configurations. In order to choose the domain limits, a combination of
reported working conditions from the literature was used [2]. The range for
parameters f, N, and Re are as follows. f was chosen from 0.25 to 0.75 with
step size of 0.02 which resulted in 26 data points. N was chosen to be {3, 4,
5, 6, 7, 8, 9, 10} and the Reynolds number was selected to be {0.01, 0.1, 1,
2.5, 5, 7.5, 10, 15, 20, 25}. In total, 2288 data points were created and
split 80%-20% for train-development sets. Furthermore, another data-set was
created to test the network, consisting of data points that were not given to
the network in the development process. This data-set was build from a
combination of f between 0.25 and 0.75 with a step size of 0.06, N of {3, 4,
5, 6}, and Re of {0.05, 1.5, 6.5, 8.5, 12.5, 18.5}, which resulted in a data-
set with 239 data points.

Next, a network architecture was designed to effectively understand the
relationship between the data. As can be seen in Table 1, nine architectures
were tested starting with a base architecture (CNN1). At first, the width of
the dense layers was increased (CNN1-3) which caused no significant
improvement despite increasing the network trainable parameters dramatically.
Then, the width and depth of convolutional layers were examined (CNN4-9) where
CNN9 exhibited the best performance and its detailed structure is presented in
Table 2. Moreover, other hyper-parameters were set to constant as follows,
batch size of 64, a learning rate of $2e-3$ for the first 100 epochs, and
$2e-4$ for the other 100 epochs. The chosen network was also gone through
another 100 epochs with a learning rate of $2e-5$ and reached the validation
loss of $5.1e-6$.

Table 1: CNN architectures study. The number of convolutional layers reported do not include the two constant layers with 64 filters in the architecture. Type | Architecture | Loss | Val-Loss | #Params
---|---|---|---|---
CNN1 | 2Dense256, 3Conv32 | $8.04e-6$ | $8.81e-6$ | 8,764,162
CNN2 | 2Dense512, 3Conv32 | $8.05e-6$ | $8.62e-6$ | 17,571,586
CNN3 | 2Dense1024, 3Conv32 | $1.30e-5$ | $1.27e-5$ | 35,972,866
CNN4 | 2Dense256, 3Conv64 | $6.93e-6$ | $8.08e-6$ | 8,948,674
CNN5 | 2Dense256, 3Conv128 | $1.73e-5$ | $1.70e-5$ | 9,538,882
CNN6 | 2Dense256, 4Conv64 | $3.04e-5$ | $1.42e-5$ | 9,022,530
CNN7 | 2Dense256, 2Conv64 | $8.22e-6$ | $8.67e-6$ | 8,874,818
CNN8 | 2Dense256, 2Conv128 | $6.95e-6$ | $7.73e-6$ | 9,243,714
CNN9 | 2Dense256, 2Conv256 | $8.47e-6$ | $7.62e-6$ | 10,423,874
Table 2: The detailed architecture of the nominated CNN9. Layer | Kernel/Stride/Pad | Output size
---|---|---
Dense 1-1 | … | 16
Dense 1-2 | … | 16
Dense 1-3 | … | 16
Concate 1-4 | … | 48
Dense 2 | … | 256
Dense 3 | … | 256
Dense 4 | … | 16384
Reshape 5-1 | … | $16\times 16\times 64$
Conv2D 5-2 |  $3\times 3$/1/same | $16\times 16\times 64$
UpSampling 6-1 | $2\times 2$ | $32\times 32\times 64$
Conv2D 6-2 |  $3\times 3$/1/same | $32\times 32\times 256$
UpSampling 7-1 | $2\times 2$ | $64\times 64\times 256$
Conv2D 7-2 |  $3\times 3$/1/same | $64\times 64\times 256$
UpSampling 8-1 | $2\times 2$ | $128\times 128\times 256$
Conv2D 8-3 |  $3\times 3$/1/same | $128\times 128\times 64$
Conv2D 8-4 |  $3\times 3$/1/same | $128\times 128\times 1$

Finally, Figure 3 depicts the performance evaluation of the chosen CNN
structure to predict the critical diameter. The accuracy of the network in
development and test set is illustrated in Figure 3a, where its mean squared
errors ($MSE$) were $2.7e-3$ and $2.6e-3$, respectively. Furthermore, The
network predictions’ $MSE$ and standard deviation ($SD$) are also demonstrated
based on Re, N, and f on the development set to detect any trends in the
errors (Figure 3b). In sum, $MSE$ and $SD$ are less than 3% and 5%,
respectively, which is an indication of the satisfactory performance of the
network. However, the error trends convey that it will be intensified with the
increase in N and f. Yet, there is no strong correlation between the error and
Re. Additionally, Figure 3c, depicts the velocity fields with the most
critical diameter prediction error from the test set ($f=0.38$, $N=8$, and
$Re=15$). By comparing the original field from the numerical study (Ground
Truth) and the neural network-generated field (Prediction), it is clear that
these fields match well. Moreover, the discrepancies between the x and y
components of velocity fields (Difference) quantitatively prove the latter
observation. The results reveal that the nominated CNN has a high accuracy
($error<1\%$) in the field generation leaving most of the critical diameter
prediction errors in the process of critical diameter computation.

Figure 3: The nominated CNN performance evaluation. a) Critical diameter
prediction accuracy in the development and test sets. b) Critical diameter
prediction errors based on $f$, $N$, and $Re$. c) CNN fields analysis at the
highest predicted $D_{c}$ error.

Then, the trained CNN was used to generate new data (data augmentation) to
facilitate the training of the direct neural network. This network helped with
the augmentation of 2288 to 10400 data points by enhancing the resolution of
Re in the new data-set by decreasing the step size to 0.499.

###  3.2 Direct neural network training

Here, the process of finding the best architecture of the FCNN for direct
prediction of critical diameter alongside its training procedure is presented.
For this particular problem, a systematic search was carried out where
networks with a maximum of 10 hidden layers (HL) and nodes (N) of up to 128
were tested. The results are given in Table 3 after training with 1000 epochs
and a learning rate of $1e-4$. The results imply that the 8 HL and 128 N
architecture display superior performance. Therefore, this structure was used
in the optimization algorithm.

Table 3: FCNN architectural study. The numbers reported here are $MSE$ values of the network’s prediction on the set. HL N | 10 | 16 | 32 | 64 | 128
---|---|---|---|---|---
4 | $1.76e-4$ | $1.04e-4$ | $7.11e-5$ | $4.81e-5$ | $4.33e-5$
5 | $1.88e-4$ | $3.54e-4$ | $7.21e-5$ | $5.63e-5$ | $4.07e-5$
6 | $6.08e-4$ | $1.10e-4$ | $6.19e-5$ | $5.67e-5$ | $4.35e-5$
8 | $1.73e-4$ | $9.16e-5$ | $6.48e-5$ | $4.55e-5$ | $3.82e-5$
10 | $1.47e-4$ | $1.02e-4$ | $5.96e-5$ | $4.82e-5$ | $4.73e-5$

The convergence process and performance of the chosen network are depicted in
Figure 4. It is shown that the network’s $MSE$ on the training set (training
loss) and the development set (validation loss) decreased per episode, Figure
4a. Considering the proximity and low values of training and validation losses
($2.8e-5$ and $3.8e-5$, respectively), it can be concluded that there were no
avoidable bias and over-fitting problems. In addition, the network predictions
of critical diameter and the actual critical diameter were compared on the
entire original and generated data-sets. The network prediction has satisfying
performance in both data-sets with the $MSE$ of $4.6e-3$ on the 2288 data-set
and $4.0e-3$ on the 10400 data-set, Figure 4b.

Figure 4: Direct neural network performance evaluation. a) The convergence
process of training and validation losses. b) The prediction accuracy in the
2288 and 10400 data-sets.

###  3.3 Design automation

After ensuring the accuracy of the direct neural network, the next step was
examining the robustness and versatility of the design automation platform
consisting of the network and evolutionary algorithm. In order to do this, a
set of conditions were used as shown in Table 4. In this table, two conditions
of the maximum stability ($\phi=0$) and the maximum flexibility ($\phi=1$)
were used in each of which f, N, or Re were set to either maximum or minimum.
In addition, for all the test subjects $D_{1}$ and $D_{2}$ were kept at 5 and
8 $\mu m$, respectively. The output of these tests indicates that the proposed
NSGA3 could apply the constraints and find the design parameters efficiently.
For example, in case 1, the f reached its minimum value with BW of 2.27 $\mu
m$, while its corresponding case with the maximum flexibility (case 7)
increased the BW up to 2.61 $\mu m$. This trend for BW occurred similarly
where the cases with $\phi=1$ always had bigger BWs. Moreover, as the last
stage of checking the method’s accuracy, the proposed configurations were also
simulated in the original software for comparison with the direct neural
network’s prediction. The relative differences are shown in Table 4 as errors
(E) being acceptably less than 4%.

Furthermore, to visually elaborate the automation outputs, Figure 5a depicts
the range of $D_{c}$ based on $Re$. The results presented in these plots are
aligned with the finding from Table 4 and it can be seen that the cases of
maximum flexibility have a higher slope than the cases of maximum stability.
In addition, the post-processing of case with minimum $N$ and maximum
flexibility (case 9) is demonstrated in Figure 5b for 10 periods. In this
figure, the particle trajectory plots illustrate how particle with bumped and
zigzag modes would behave within the full length of the device and the
recurrence maps provide the exact lateral positions of the particles at the
input and output of each period.

Table 4: Design automation results for chosen samples. Inputs |  | Outputs |  |  |
---|---|---|---|---|---
No. | $\phi$ | f | N | Re | f | N | Re | G | $D_{c}$ | BW | E%
1 | 0 | Min | - | - | 0.25 | 9.00 | 1.28 | 15.79 | 6.51 | 2.27 | 3.20
2 | 0 | Max | - | - | 0.75 | 10.00 | 21.91 | 19.07 | 6.49 | 3.41 | 2.37
3 | 0 | - | Min | - | 0.51 | 3.00 | 0.13 | 8.80 | 6.49 | 3.41 | 0.82
4 | 0 | - | Max | - | 0.46 | 10.00 | 5.56 | 19.03 | 6.50 | 0.90 | 1.13
5 | 0 | - | - | Min | 0.41 | 10.00 | 0.01 | 21.62 | 6.50 | 1.63 | 1.73
6 | 0 | - | - | Max | 0.47 | 8.00 | 25.00 | 17.55 | 6.50 | 1.38 | 0.54
7 | 1 | Min | - | - | 0.25 | 10.00 | 18.75 | 20.45 | 6.49 | 2.61 | 0.56
8 | 1 | Max | - | - | 0.75 | 4.00 | 0.02 | 16.94 | 6.50 | 10.44 | 3.98
9 | 1 | - | Min | - | 0.75 | 4.00 | 4.29 | 14.28 | 6.50 | 8.81 | 2.20
10 | 1 | - | Max | - | 0.75 | 10.00 | 18.99 | 21.51 | 6.50 | 3.85 | 1.24
11 | 1 | - | - | Min | 0.75 | 4.00 | 0.03 | 16.93 | 6.50 | 10.44 | 3.89
12 | 1 | - | - | Max | 0.74 | 5.00 | 21.66 | 8.85 | 6.50 | 4.50 | 2.40
Figure 5: Design automation results. a) Optimized devices’ performance plots
showing the relationship between $D_{c}$ and $Re$. b) An example of $D_{1}=5$
$\mu m$ and $D_{2}=8$ $\mu m$ particle trajectories and their corresponding
recurrence maps.

##  4 Conclusion

This work presents a design automation platform for DLD devices by
incorporating the deep neural network’s ability to understand their underlying
physics. To achieve this goal, first, the DLD physics was numerically
simulated whereafter extracting the data, they were modified for a CNN
architecture. Furthermore, using a CNN enabled the fast prediction of the flow
fields to be readily used for particle trajectory simulation and data
augmentation; however, they were not utilized directly for critical diameter
prediction. Instead, a fully-connected neural network was trained to
understand the direct relationship between the device properties and its
critical diameter. Finally, the neural network’s fast predictability was
embedded into an NSGA3 multi-objective optimization to create the design
automation tool. The main purpose of the automation platform was to take the
particle diameters and suggest a design and working condition for their
separation. Later, more features were allowed to be entered by the user to
have more control over the suggested design and the behavior of the device
after the fabrication.

Overall, each design automation component exhibits an acceptable level of
accuracy ensuring the reliability of the tool to suggest suitable design
parameters. The convolutional neural network which was trained on a small
data-set reached the validation loss of $5.1e-6$ and the critical diameter
extracted from this predicted field had the mean squared error of $2.6e-3$.
Moreover, the fully connected neural network that was used for the direct
prediction of critical diameter achieved $MSE=4.6e-3$. Knowing the desirable
performance of each component, the automation platform was tested for 12
critical cases. In the last stage, the predicted critical diameters of each
case were revalidated by the original software as their maximum error did not
exceed 4%. This method presents itself as highly accurate where the stem of
the errors is likely due to the simplifying particle trajectory simulation
that is used for critical diameter extraction. Moreover, this approach is
functional and generalizable to other physics and applications that require
field prediction and particle trajectory simulations. In the case of DLD
devices, there is a possibility of transferring learning to other cases with
different pillar shapes and patterns. Future work can be enhancing the
simulation technique, employing experimental data, and including other
prominent shapes in DLD devices.

## Appendix

###  A Particle-pillar contact simulation

Wall function was used to simulate particle interactions with obstacles. It
stores the shortest distance from any position to pillars in a discretized
periodic domain. Then, this information serves as an efficient mapping between
the spatial positions of particles to their distance from pillars, since it is
only calculated once for each geometry. Moreover, the normal and tangential
unit vector fields for velocity modification in the contact point can be
calculated by extracting the wall function gradient and its perpendicular
direction, respectively. Figure A.1 shows an example of the wall distance
function and the normal vector field.

Figure A.1: a) Wall distance function contour and b) wall normal vector field.

###  B Critical diameter extraction with Newtonian root finding method

The Newtonian root finding method was used for extracting the critical
diameter from particle trajectory simulation. Here, the detailed explanation
of this method is provided in the Algorithm 1 where it gets the desired
tolerance as input and derives the critical diameter.

Algorithm 1 Finding critical diameter with Newtonian method.

Input $Tolerance$
Output $d_{c}$

$g\leftarrow 1-f$ $\triangleright$ Calculate $g$

$d_{1}\leftarrow 0.1\times g$ $\triangleright$ Lowest diameter

$d_{2}\leftarrow 0.95\times g$ $\triangleright$ Highest diameter

$d_{1,value}\leftarrow$Particle Tracing$(d_{1})$ $\triangleright$ Perform
particle tracing with $d_{1}$

$d_{2,value}\leftarrow$Particle Tracing$(d_{2})$ $\triangleright$ Perform
particle tracing with $d_{2}$

if $d_{1,value}=-1\And d_{2,value}=1$ then

while $(d_{2}-d_{1})\geq Tolerance$ do

$d=0.5\times(d_{1}+d_{2})$ $\triangleright$ Halve the diameter difference

$d_{value}\leftarrow$Particle Tracing$(d)$ $\triangleright$ Perform particle
tracing with $d$

if $d_{value}=1$ then

$d_{2}\leftarrow d$ $\triangleright$ Update $d_{2}$

else

$d_{1}\leftarrow d$ $\triangleright$ Update $d_{1}$

end if

$d_{c}\leftarrow d$ $\triangleright$ Obtain critical diameter

end while

else

$d_{c}\leftarrow None$ $\triangleright$ No critical diameter was found for
separation

end if

###  C Critical diameter trends in the original data-set

Figure C.1 depicts the relationship between all the effective parameters and
critical diameter in the original data-set in four f values to provide a
better understanding of the DLD physics and criteria that should be used in
automation step. The overall trends of the critical diameter changes with an
increase of f: at low f values, in almost all N values there is a peak around
Re=5, followed by the critical diameter value decreases. Nevertheless, for N=3
there are no critical diameters at all. However, the peak vanishes as the f
value grows and it only shows a positive slope. Moreover, at N=3 some area
appears to have a critical diameter around a low Reynolds number and this area
expands as f value increases. Furthermore, the area with a slope near zero
offers the best stability since the critical diameter would not change much,
and separation will occur regardless of the fluctuation in Re. On the other
hand, it appears that higher f and lower N have a steeper slope and provide a
better chance of flexible design.

Figure C.1: The effect of all influential parameters on critical diameter.

## References

  * [1] J. McGrath, M. Jimenez and H. Bridle  “Deterministic lateral displacement for particle separation: a review”  In _Lab Chip_ 14 The Royal Society of Chemistry, 2014, pp. 4139–4158
  * [2] Hao Tang et al.  “Geometric structure design of passive label-free microfluidic systems for biological micro-object separation”  In _Microsystems & Nanoengineering_ 8.1 Nature Publishing Group, 2022, pp. 1–28
  * [3] Lotien Richard Huang, Edward C Cox, Robert H Austin and James C Sturm  “Continuous particle separation through deterministic lateral displacement”  In _Science_ 304.5673 American Association for the Advancement of Science, 2004, pp. 987–990
  * [4] David W Inglis, John A Davis, Robert H Austin and James C Sturm  “Critical particle size for fractionation by deterministic lateral displacement”  In _Lab on a Chip_ 6.5 Royal Society of Chemistry, 2006, pp. 655–658
  * [5] John A Davis et al.  “Deterministic hydrodynamics: taking blood apart”  In _Proceedings of the National Academy of Sciences_ 103.40 National Acad Sciences, 2006, pp. 14779–14784
  * [6] Kerwin Kwek Zeming, Thoriq Salafi, Chia-Hung Chen and Yong Zhang  “Asymmetrical deterministic lateral displacement gaps for dual functions of enhanced separation and throughput of red blood cells”  In _Scientific reports_ 6.1 Nature Publishing Group, 2016, pp. 1–10
  * [7] Sung-Cheol Kim et al.  “Broken flow symmetry explains the dynamics of small particles in deterministic lateral displacement arrays”  In _Proceedings of the National Academy of Sciences_ 114.26 National Acad Sciences, 2017, pp. E5034–E5041
  * [8] Yanjuan Wang, Junsheng Wang, Yanghao Wu and Jianhong Dong  “The automatic and high-throughput purification and enrichment of microalgae cells using deterministic lateral displacement arrays with different post shapes”  In _Journal of Chemical Technology & Biotechnology_ 96.8 Wiley Online Library, 2021, pp. 2228–2237
  * [9] Kerwin Kwek Zeming, Shashi Ranjan and Yong Zhang  “Rotational separation of non-spherical bioparticles using I-shaped pillar arrays in a microfluidic device”  In _Nature communications_ 4.1 Nature Publishing Group, 2013, pp. 1–8
  * [10] Kerwin Kwek Zeming et al.  “Microfluidic label-free bioprocessing of human reticulocytes from erythroid culture”  In _Lab on a Chip_ 20.18 Royal Society of Chemistry, 2020, pp. 3445–3460
  * [11] Ji-chul Hyun, Jaeyub Hyun, Semyung Wang and Sung Yang  “Improved pillar shape for deterministic lateral displacement separation method to maintain separation efficiency over a long period of time”  In _Separation and Purification Technology_ 172 Elsevier, 2017, pp. 258–267
  * [12] Brian M Dincau et al.  “Deterministic lateral displacement (DLD) in the high Reynolds number regime: High-throughput and dynamic separation characteristics”  In _Microfluidics and Nanofluidics_ 22.6 Springer, 2018, pp. 1–8
  * [13] Ali Lashkaripour et al.  “Machine learning enables design automation of microfluidic flow-focusing droplet generation”  In _Nature communications_ 12.1 Nature Publishing Group, 2021, pp. 1–14
  * [14] David McIntyre, Ali Lashkaripour, Polly Fordyce and Douglas Densmore  “Machine learning for microfluidic design and control”  In _Lab on a Chip_ Royal Society of Chemistry, 2022
  * [15] Qinghai Bai  “Analysis of particle swarm optimization algorithm”  In _Computer and information science_ 3.1 Canadian Center of ScienceEducation, 2010, pp. 180
  * [16] EM Arun Sankar, Mohammad Shahab and Raghunathan Rengaswamy  “Spacing optimization for active droplet sorting in microfluidic networks using genetic algorithm”  In _Industrial & Engineering Chemistry Research_ 60.4 ACS Publications, 2021, pp. 1699–1708
  * [17] Xiao-jun Liu, Hong Yi and Zhong-hua Ni  “Application of ant colony optimization algorithm in process planning optimization”  In _Journal of Intelligent Manufacturing_ 24.1 Springer, 2013, pp. 1–13
  * [18] Xijun Huang, Chuanpei Xu and Long Zhang  “An efficient algorithm for optimizing the test path of digital microfluidic biochips”  In _Journal of Electronic Testing_ 36.2 Springer, 2020, pp. 205–218
  * [19] Yusliza Yusoff, Mohd Salihin Ngadiman and Azlan Mohd Zain  “Overview of NSGA-II for optimizing machining process parameters”  In _Procedia Engineering_ 15 Elsevier, 2011, pp. 3978–3983
  * [20] Michael I Jordan and Tom M Mitchell  “Machine learning: Trends, perspectives, and prospects”  In _Science_ 349.6245 American Association for the Advancement of Science, 2015, pp. 255–260
  * [21] Maziar Raissi and George Em Karniadakis  “Hidden physics models: Machine learning of nonlinear partial differential equations”  In _Journal of Computational Physics_ 357 Elsevier, 2018, pp. 125–141
  * [22] Justin Sirignano and Konstantinos Spiliopoulos  “DGM: A deep learning algorithm for solving partial differential equations”  In _Journal of computational physics_ 375 Elsevier, 2018, pp. 1339–1364
  * [23] Jiequn Han, Arnulf Jentzen and Weinan E  “Solving high-dimensional partial differential equations using deep learning”  In _Proceedings of the National Academy of Sciences_ 115.34 National Acad Sciences, 2018, pp. 8505–8510
  * [24] Yunfei Li, Juntao Chang, Chen Kong and Wen Bao  “Recent progress of machine learning in flow modeling and active flow control”  In _Chinese Journal of Aeronautics_ Elsevier, 2021
  * [25] Xiaowei Jin, Peng Cheng, Wen-Li Chen and Hui Li  “Prediction model of velocity field around circular cylinder over various Reynolds numbers by fusion convolutional neural networks based on pressure on the cylinder”  In _Physics of Fluids_ 30.4 AIP Publishing LLC, 2018, pp. 047105
  * [26] Sangseung Lee and Donghyun You  “Data-driven prediction of unsteady flow over a circular cylinder using deep learning”  In _Journal of Fluid Mechanics_ 879 Cambridge University Press, 2019, pp. 217–254
  * [27] Maziar Raissi, Paris Perdikaris and George E Karniadakis  “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations”  In _Journal of Computational physics_ 378 Elsevier, 2019, pp. 686–707
  * [28] Maziar Raissi, Zhicheng Wang, Michael S Triantafyllou and George Em Karniadakis  “Deep learning of vortex-induced vibrations”  In _Journal of Fluid Mechanics_ 861 Cambridge University Press, 2019, pp. 119–137
  * [29] Salvatore Cuomo et al.  “Scientific Machine Learning through Physics-Informed Neural Networks: Where we are and What’s next”  In _arXiv preprint arXiv:2201.05624_ , 2022
  * [30] Vinothkumar Sekar, Qinghua Jiang, Chang Shu and Boo Cheong Khoo  “Fast flow field prediction over airfoils using deep learning approach”  In _Physics of Fluids_ 31.5 AIP Publishing LLC, 2019, pp. 057103
  * [31] Junhyuk Kim and Changhoon Lee  “Prediction of turbulent heat transfer using convolutional neural networks”  In _Journal of Fluid Mechanics_ 882 Cambridge University Press, 2020
  * [32] Eric Gioe, Mohammed Raihan Uddin, Jong-Hoon Kim and Xiaolin Chen  “Deterministic Lateral Displacement (DLD) Analysis Tool Utilizing Machine Learning towards High-Throughput Separation”  In _Micromachines_ 13.5 MDPI, 2022, pp. 661
  * [33] Alexey Dosovitskiy, Jost Tobias Springenberg, Maxim Tatarchenko and Thomas Brox  “Learning to generate chairs, tables and cars with convolutional networks”  In _IEEE transactions on pattern analysis and machine intelligence_ 39.4 IEEE, 2016, pp. 692–705
  * [34] J. Blank and K. Deb  “pymoo: Multi-Objective Optimization in Python”  In _IEEE Access_ 8, 2020, pp. 89497–89509


# Maximum Energy Growth Rate in Dilute Quantum Gases

Ran Qi<EMAIL_ADDRESS>Department of Physics, Renmin University of China,
Beijing, 100872, P. R. China  Zheyu Shi  Key State Laboratory of Precision
Spectroscopy, East China Normal University, Shanghai 200062, China  Hui Zhai
<EMAIL_ADDRESS>Institute for Advanced Study, Tsinghua University,
Beijing 100084, China

###### Abstract

In this letter we study how fast the energy density of a quantum gas can
increase in time, when the inter-atomic interaction characterized by the
$s$-wave scattering length $a_{\text{s}}$ is increased from zero with
arbitrary time dependence. We show that, at short time, the energy density can
at most increase as $\sqrt{t}$, which can be achieved when the time dependence
of $a_{\text{s}}$ is also proportional to $\sqrt{t}$, and especially, a
universal maximum energy growth rate can be reached when $a_{\text{s}}$ varies
as $2\sqrt{\hbar t/(\pi m)}$. If $a_{\text{s}}$ varies faster or slower than
$\sqrt{t}$, it is respectively proximate to the quench process and the
adiabatic process, and both result in a slower energy growth rate. These
results are obtained by analyzing the short time dynamics of the short-range
behavior of the many-body wave function characterized by the contact, and are
also confirmed by numerical solving an example of interacting bosons with
time-dependent Bogoliubov theory. These results can also be verified
experimentally in ultracold atomic gases.

The ability of tuning interactions between particles is a major advantage of
ultracold atomic systems FR1 ; FR2 . Especially, by ultilizing magnetic and
optical tools, the interaction strength between atoms, usually characterized
by the $s$-wave scattering length $a_{\text{s}}$, can be tuned in a time scale
much faster than the many-body relaxation time. This has led to a number of
interesting ultracold atomic experiments reported in recent years, such as
universal quench dynamics observed by quenching interaction to unitarity Jin ;
cambridge , and coherent excitation of the Higgs model in superfluid Fermi
gases and the Bogoliubov quasi-particles in Bose condensate by periodically
modulating interactions Higgs ; Cheng1 ; Cheng2 ; Cheng3 . These experimental
progresses are also accompanied by lots of theoretical interests on studying
non-equilibrium dynamics driven by time-dependent interactions theory1 ;
theory2 ; theory3 ; theory4 ; theory5 ; theory6 ; theory7 ; theory8 ; theory9
; theory10 ; theory11 ; theory12 ; theory13 ; theory14 ; theory15 ; theory16 ;
theory17 ; theory18 ; theory19 ; theory20 ; theory21 ; theory22 .

Motivated by these progresses, here we address a fundamental issue whether
there is a universal upper limit for the energy increasing rate. To be
concrete, suppose that we start with a non-interacting quantum gas with
$a_{\text{s}}=0$ and then change $a_{\text{s}}$ in time, and suppose that
$a_{\text{s}}(t)$ can be controlled in any function form, the question is
whether there is an upper bound for the rate of how fast the total energy can
increase as a function of time. In this letter we show that there does exist
such a universal rate limit, as far as the initial growth rate is concerned.
This result is quite counter-intuitive, because normally the interaction
energy increases as the interaction strength increases. Thus, intuitively, one
would think that a faster increasing of interaction strength should result in
a faster increasing of interaction energy, and consequently, a faster
increasing of the total energy. Since we consider that $a_{\text{s}}$ can be
increased as fast as one wants, it seems to indicate that there should not be
such a bound.

However, our results show that this intuition is not correct. Before
presenting rigorous mathematical statement, we first emphasize that our result
is closely tied to a key quantity of ultracold atomic gases called the contact
contact1 ; contact2 ; contact3 ; contact4 ; contact5 ; contact6 ; contact7 ;
bose_contact ; contact_pwave . It is now well known that, for quantum gases
with zero-range interactions, one can define contact $\mathcal{C}$ through the
short-range behavior of the many-body wave function when any two atoms are
brought close to each other, or equivalently, through the high-momentum tail
of the momentum distribution. It has been shown that the total energy of a
quantum gas is directly related to the contact contact1 ; contact2 ; contact3
; contact4 ; contact5 ; contact6 ; contact7 ; bose_contact ; contact_pwave .

To gain an intuitive understanding of our results, let us first consider two
limits. On one limit, the fastest change of the interaction strength is the
quench process, during which $a_{\text{s}}$ instantaneously jumps from zero to
any non-zero value. However, it can be shown that the contact does not change
and retains zero right after the quench contact2 , and therefore, the total
energy also does not change after the quench. This means that the fastest
change of interaction actually does not result in a fast change of the total
energy, and in contrast, the interaction energy does not change at all. On the
opposite limit, we can consider an adiabatic varying of the interaction
strength, during which the interaction energy does vary in time but it varies
adiabatically with sufficiently slow rate. The physical pictures in these two
limits motivate us to expect a universal maximum growth rate driven by
intermediate speed of varying the interaction strength.

General Expression for the Contact Growth. Here we consider a uniform Bose gas
or spin-1/2 Fermi gas starting from any non-interacting state
$|\Psi_{0}\rangle$ at $t=0$, and then the $s$-wave scattering length
$a_{\text{s}}(t)$ can vary with arbitrary time dependence. Below we use $n$
and $n_{\sigma}$ to denote the densities of bosons and fermions with
spin-$\sigma$ ($\sigma=\uparrow,\downarrow$), respectively, and $\hat{\psi}$
and $\hat{\psi}_{\sigma}$ to denote boson operator and fermion operator with
spin-$\sigma$, respectively. One of the main results of this work states as
follows:

In the short-time limit, the dynamics of the contact $\mathcal{C}(t)$ is given
by

$\mathcal{C}(t)=g_{2}(0)|\eta(t)|^{2}.$ (1)

Here $g_{2}({\bf r})$ is defined as
$\langle\Psi_{0}|\hat{\psi}^{{\dagger}}(\frac{\mathbf{r}}{2})\hat{\psi}^{{\dagger}}(-\frac{\mathbf{r}}{2})\hat{\psi}(-\frac{\mathbf{r}}{2})\hat{\psi}(\frac{\mathbf{r}}{2})|\Psi_{0}\rangle$
for bosons and
$\langle\Psi_{0}|\hat{\psi}^{{\dagger}}_{\uparrow}(\frac{\mathbf{r}}{2})\hat{\psi}^{{\dagger}}_{\downarrow}(-\frac{\mathbf{r}}{2})\hat{\psi}_{\downarrow}(-\frac{\mathbf{r}}{2})\hat{\psi}_{\uparrow}(\frac{\mathbf{r}}{2})|\Psi_{0}\rangle$
for fermions, and $g_{2}(0)$ means $g_{2}({\bf r})$ evaluated at ${\bf r}=0$.
Especially, if $|\Psi_{0}\rangle$ is the non-interacting ground state, then
$g_{2}(\mathbf{r})=n^{2}$ or $n_{\uparrow}n_{\downarrow}$ for the Bose or the
spin-$1/2$ Fermi gas. The key result is that the function $\eta(t)$ obeys the
following integral equation

$\displaystyle\left[\hat{L}+\frac{1}{4\pi a_{\text{s}}(t)}\right]\eta(t)=-1,$
(2)

where $\hat{L}$ is a linear operator acting on $\eta(t)$ as

$\displaystyle\hat{L}\eta(t)=$
$\displaystyle\left(\frac{m}{\hbar}\right)^{\frac{1}{2}}\frac{1}{8\pi^{3/2}\sqrt{i}}\times$
$\displaystyle\lim_{\epsilon\rightarrow
0^{+}}\left[\int_{0}^{t-\epsilon}\frac{\eta(\tau)}{(t-\tau)^{\frac{3}{2}}}d\tau-\frac{2\eta(t)}{\sqrt{\epsilon}}\right].$
(3)

This result is motivated by solving the two-body problem, which satisfies the
following Schr$\ddot{\text{o}}$dinger equation in the relative coordinate
${\bf r}$ frame as

$\displaystyle i\hbar\frac{\partial\psi}{\partial
t}=-\frac{\hbar^{2}\nabla^{2}\psi}{m}+\frac{4\pi\hbar^{2}a_{\text{s}}(t)}{m}\delta(\mathbf{r})\frac{\partial}{\partial
r}r\psi.$ (4)

Starting from an initial state $\psi({\bf r})=1/\sqrt{V}$ ($V$ is the total
volume of the system), the time evolution of the wave function always obeys
the following asymptotic form in the short-range $r\rightarrow 0$ limit, that
is

$\displaystyle\psi(\mathbf{r},t)=\frac{\eta(t)}{4\pi\sqrt{V}}\left[\frac{1}{r}-\frac{1}{a_{\text{s}}(t)}\right]+O(r),$
(5)

and it can be shown that $\eta(t)$ satisfies Eq. 28 supple . Generalizing this
result from the two-body problem to the many-body problem utilizes the short-
time expansion and is quite straightforward, which yields Eq. 1 supple . Here
we note that for the two-body problem, $\eta(t)$ satisfies Eq. 28 for all time
scales, but for the many-body problem, the result is only valid for the short-
time scale. Here short-time is defined as the time scale much shorter than the
typical many-body time scale $t_{n}=\hbar/E_{\text{n}}$, where
$E_{\text{n}}=\hbar^{2}k_{\text{n}}^{2}/(2m)$ and $k_{n}=(6\pi^{2}n)^{1/3}$
(with $n$ replaced by $n_{\sigma}$ for fermions). In other word, in such short
time scale, the short-range behavior of many-body wave function is still
dominated by the two-body physics.

Figure 1: (a1-a2) The time-dependence of the scattering length
$a_{\text{s}}(t)$ (in unit of $l_{0}$) with different power-law functions of
Eq. 6. (a1) $\alpha=\frac{1}{4},~{}\frac{1}{2},~{}\frac{3}{4}$, and $\beta$ is
fixed at $\beta=1$. (a2) $\alpha$ is fixed at $\alpha=\frac{1}{2}$ and
$\beta=\frac{1}{2},~{}1.128,~{}1.596$ and $4$. (b1-b2) The short time behavior
of the contact $\mathcal{C}$ (in unit of $g_{2}(0)l_{0}^{2}$) with
$a_{\text{s}}(t)$ plotted in (a1) and (a2), respectively. (c1-c2) The time-
dependence of the energy density change $\delta\mathcal{E}$ (in unit of
$g_{2}(0)l_{0}\hbar^{2}/m$) with $a_{\text{s}}(t)$ plotted in (a1) and (a2),
respectively.

Contact Growth Rate. Here, without loss of generality, we consider that
$a_{\text{s}}(t)$ grows from zero to a positive value in a power-law function
as

$a_{\text{s}}(t)=\sqrt{2}\beta l_{0}\left(\frac{t}{t_{0}}\right)^{\alpha},$
(6)

where $l_{0}$ is an arbitrary length unit and $t_{0}$ is the time units, and
$l_{0}$ and $t_{0}$ are both related to the same energy unit as
$\hbar/t_{0}=\hbar^{2}/(2ml^{2}_{0})$. $\alpha,\beta$ are two constants
describing the power and the coefficient, respectively and a factor $\sqrt{2}$
is introduced just for the later convenience. The operator $\hat{L}$ has an
important property that

$\displaystyle\hat{L}t^{\alpha}=\left(\frac{m}{\hbar}\right)^{\frac{1}{2}}B(\alpha)t^{\alpha-\frac{1}{2}},$
(7)

where $B(\alpha)$ is a constant given by
$B(\alpha)=i^{3/2}\Gamma(\alpha+1)/(4\pi\Gamma(\alpha+1/2))$. That is to say,
suppose $\eta(t)$ is a power-law function in $t$, when $\hat{L}$ acts on
$\eta(t)$, it lowers the power of $\eta(t)$ by $1/2$. This property plays a
crucial role in the following analysis because it means whether $\alpha$ in
Eq. 6 is greater or smaller than $1/2$ makes significant difference.

Case I: $\alpha>1/2$. In this case, the $1/(4\pi a_{\text{s}})$ term dominates
Eq. 28, and thus, to the leading order of $t$, $\eta(t)$ and $\mathcal{C}(t)$
are given by

$\eta(t)=-4\pi a_{\text{s}}(t);\ \
\mathcal{C}(t)=16\pi^{2}a^{2}_{\text{s}}(t)g_{2}(0).$ (8)

This is consistent with the adiabatic regime where the physical quantities
only depend on the instantaneous scattering length at time $t$.

Case II: $\alpha<1/2$. In this case, the $\hat{L}$ term dominates Eq. 28, and
thus, to the leading order of $t$, $\eta(t)$ and $\mathcal{C}(t)$ are given by

$\eta(t)=-\left(\frac{\hbar}{m}\right)^{\frac{1}{2}}\frac{1}{B(1/2)}\sqrt{t};\
\ \mathcal{C}(t)=\frac{\hbar}{m}\frac{g_{2}(0)}{B(1/2)^{2}}t,$ (9)

where $B(1/2)=-1/(8\sqrt{i\pi})$. Surprisingly, in this case this result shows
that the growth of contact at the short-time is independent of parameters
$l_{0}$, $t_{0}$, $\alpha$ and $\beta$ in Eq. 6. That is to say, it is
independent of how fast $a_{\text{s}}$ varies in time. Even if $l_{0}$ or
$\beta$ is infinitely large, or $\alpha$ is infinitesimally small, and then
$a_{\text{s}}(t)$ initially grows infinitely fast, the contact always grows
linearly in time with a constant rate. This means that as long as
$\alpha<1/2$, the short-range physics at the short time is the same as a
quench process where the scattering length instantaneously jumps to unitarity.

Case III: $\alpha=1/2$. In this case, Eq. 6 becomes

$a_{\text{s}}(t)=\beta\sqrt{\frac{\hbar t}{m}},$ (10)

By dimension analysis, it is easy to see that $l_{0}$ and $t_{0}$ cancel each
other and only the coefficient $\beta$ enters the expression. In this case,
both $\hat{L}$ term and the $1/(4\pi a_{\text{s}})$ term are equally
important. Also to the leading order of $t$, we obtain

$\eta(t)=-A(\beta)\sqrt{t};\ \ \mathcal{C}(t)=|A(\beta)|^{2}g_{2}(0)t,$ (11)

where $A(\beta)$ is also a constant given by

$A(\beta)=\left(\frac{\hbar}{m}\right)^{\frac{1}{2}}\frac{1}{B\left(\frac{1}{2}\right)+\frac{1}{4\pi\beta}}.$
(12)

As one can see from here, this is a critical case. In the Case III, by taking
$\beta\rightarrow\infty$, Eq. 11 recovers Eq. 9, consistent with the quench
limit, and by taking $\beta\rightarrow 0$, Eq. 11 recovers Eq. 8, consistent
with the adiabatic limit.

Here an important point is that $|A(\beta)|^{2}$ is not a monotonic function
in $\beta$. For a given initial state, $g_{2}(0)$ is fixed, and we can then
define the initial growth rate for contact as
$v_{\mathcal{C}}=\lim_{t\rightarrow 0}d\mathcal{C}(t)/dt$. As one can see from
Eq. 8, $v_{\mathcal{C}}=0$ for case I. And for both the Case II and the Case
III, $v_{\text{C}}$ is a constant, given by $|A(\beta)|^{2}g_{2}(0)$ for the
Case III and $|A(\beta=\infty)|^{2}g_{2}(0)$ for the Case II. It turns out
that $|A(\beta)|^{2}$ reaches its maximum at
$\beta_{\text{c1}}=2\sqrt{2/\pi}\approx 1.596$, at which
$v^{\text{max}}_{\mathcal{C}}=(\hbar/m)128\pi g_{2}(0)$.

Figure 2: Initial growth rate for contact (a) and for energy (b) as a function
of $\beta$ for $a_{\text{s}}(t)=\beta\sqrt{\hbar t/m}$. Arrows mark
$\beta_{\text{c1}}$ and $\beta_{\text{c2}}$ where the maximum contact growth
rate and the maximum energy growth rate are reached. $\upsilon_{C}$ and
$\upsilon_{E}$ are plotted in units of $g_{2}(0)\hbar/m$ and
$g_{2}(0)\sqrt{\hbar^{3}/m}$ respectively.

Energy Growth Rate. The total energy density of a uniform zero-range
interacting quantum gas can be measured through its momentum distribution
$n_{{\bf k}}$. For example, for spin-$1/2$ fermions, it is given by

$\mathcal{E}=\int\frac{d^{3}k}{(2\pi)^{3}}\epsilon_{{\bf k}}\left(n_{{\bf
k}}-\frac{2\mathcal{C}}{k^{4}}\right)+\frac{\mathcal{C}}{4\pi ma_{s}},$ (13)

where $n_{{\bf k}}=n_{{\bf k}\uparrow}+n_{{\bf k}\downarrow}$, $\epsilon_{{\bf
k}}=\hbar^{2}{\bf k}^{2}/(2m)$, and the contact $\mathcal{C}$ is related to
$n_{{\bf k}\sigma}$ through
$\mathcal{C}\equiv\lim_{k\rightarrow\infty}k^{4}n_{{\bf k}\sigma}$ contact1 .
The same expression, replacing all $\mathcal{C}$ by $\mathcal{C}/2$, also
holds for the spinless Bose gas as long as the three-body contact can be
ignored bose_contact .

On the other hand, there is a direct relation between the time evolution of
the energy and the contact. For spin-$1/2$ fermions it is given as

$\displaystyle\frac{d}{dt}\mathcal{E}(t)=\frac{\hbar^{2}\mathcal{C}(t)}{4\pi
ma_{\text{s}}^{2}(t)}\frac{da_{\text{s}}}{dt}.$ (14)

For spinless bosons, an extra $1/2$ factor should also be added in the r.h.s.
of Eq. 14. Therefore, based on the contact growth discussed above, we can
determine the energy growth.

Case I: $\alpha>1/2$. With the help of Eq. 8, one can obtain that

$\delta\mathcal{E}(t)=\frac{4\pi\hbar^{2}a_{\text{s}}(t)}{m}g_{2}(0).$ (15)

where $\delta\mathcal{E}(t)=\mathcal{E}(t)-\mathcal{E}(t=0)$. This result
again shows that the physics in this regime is consistent with adiabatic
regime where the energy is determined by the instantaneous scattering length.
Since $\alpha>1/2$, the energy increases slower than $\sqrt{t}$ at the short
time.

Case II: $\alpha<1/2$. In this regime, Eq. 9 gives rise to

$\delta\mathcal{E}(t)=\frac{16\sqrt{2}\alpha}{\beta(1-\alpha)}\frac{\hbar^{2}}{m}g_{2}(0)l_{0}\left(\frac{t}{t_{0}}\right)^{1-\alpha}.$
(16)

Since $\alpha<1/2$, the energy also increases slower than $\sqrt{t}$ at the
short time. When taking the $\alpha\rightarrow 0$ limit, or
$\beta\rightarrow\infty$ or $l_{0}\rightarrow\infty$ limit,
$\delta\mathcal{E}(t)\rightarrow 0$, and it is consistent with the fact there
is no energy change for the quench process as discussed above.

Case III: $\alpha=1/2$. In this regime, Eq. 11 yields

$\delta\mathcal{E}(t)=\sqrt{\frac{\hbar^{3}}{m}}\frac{|A(\beta)|^{2}}{4\pi\beta}g_{2}(0)\sqrt{t}.$
(17)

It is in this case that the energy growth at the short time is the fastest.
Now we can define an energy growth rate $v_{E}=\lim_{t\rightarrow
0}d\mathcal{E}(t)/d\sqrt{t}$. For case I and II, this rate is zero. In case
II, $v_{E}$ is given by
$\sqrt{\hbar^{3}/m}|A(\beta)|^{2}g_{2}(0)/(4\pi\beta)$, which reaches its
maximum at $\beta_{\text{c2}}=2/\sqrt{\pi}\approx 1.128$ with
$v^{\text{max}}_{E}=4(2+\sqrt{2})\sqrt{\pi}g_{2}(0)\sqrt{\hbar^{3}/m}\approx
24.2g_{2}(0)\sqrt{\hbar^{3}/m}$. Note that this value of $v^{\text{max}}_{E}$
applies for the spin-$1/2$ Fermi gas, and for the spinless Bose gas an extra
1/2 factor should be added.

This maximum energy growth rate is the main result of this work. We note that,
although this result is obtained by assuming power-law function of
$a_{\text{s}}(t)$ and by considering positive $a_{\text{s}}(t)$, it can be
extended to other function forms, such as including the logarithmic function
corrections, and including the situations where $a_{\text{s}}$ varies to
negative values. The results discussed above are summarized in Fig. 1 and Fig.
2. Fig. 1(a1) and (a2) show different power-law function of $a_{\text{s}}(t)$
given by Eq. 6, either with different power $\alpha$, or with different
coefficient $\beta$ and fixed $\alpha=1/2$. Fig. 1(b1) and (b2) show the
corresponding contact growth, and Fig. 1(c1) and (c2) show the corresponding
energy growth, using spinless bosons as an example. It clearly shows that a
faster increasing of $a_{\text{s}}$ does not necessarily lead to a faster
increasing of the contact and the energy density. One can see that for
different powers, $\alpha=1/2$ gives the fastest contact growth and energy
growth at the short time. And for $\alpha$ fixed at $1/2$,
$\beta=\beta_{\text{c1}}$ yields the fastest contact growth and
$\beta=\beta_{\text{c2}}$ yields the fastest energy growth, as also shown in
Fig. 2.

Figure 3: Dynamics of the total energy density of Bose gas for
$a_{\text{s}}(t)$ with different power-law functions of Eq. 6. (a) $\beta=1$
and $\alpha=\frac{1}{4},~{}\frac{1}{2},~{}\frac{3}{4}$. (b)
$\alpha=\frac{1}{2}$ and $\beta=4,~{}0.5,~{}1.128$. $\delta\mathcal{E}$ is
plotted in unit of $n^{2}l_{0}\hbar^{2}/m$ and we have set $t_{0}=t_{n}$ and
thus $l_{0}=1/k_{n}$ in the numerical calculation.

Example. The analysis above is based on the short time expansion. Here, we
consider a concrete example of spinless bosons, which can be described by the
following time-dependent Hamiltonian

$\displaystyle\hat{H}(t)=\sum_{\mathbf{k}}\epsilon_{\mathbf{k}}\hat{b}_{\mathbf{k}}^{{\dagger}}\hat{b}_{\mathbf{k}}+\frac{U(t)}{2V}\sum_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}\hat{b}_{\mathbf{k}}^{{\dagger}}\hat{b}_{\mathbf{q}-\mathbf{k}}^{{\dagger}}\hat{b}_{\mathbf{q}-\mathbf{k}^{\prime}}\hat{b}_{\mathbf{k}^{\prime}},$
(18)

where $\hat{b}_{{\bf k}}$ are boson creation operators with momentum ${\bf
k}$. $U(t)$ is related to $a_{\text{s}}(t)$ through the renormalization
relation

$\frac{1}{U(t)}=\frac{m}{4\pi\hbar^{2}a_{\text{s}}(t)}-\frac{1}{V}\sum_{\mathbf{k}}\frac{1}{2\epsilon_{\mathbf{k}}}.$
(19)

We solve this Hamiltonian by adopting the Bogoliubov-type variational ansatz
as

$\displaystyle|\Psi(t)\rangle=\mathcal{N}(t)\exp\left[g_{0}(t)\hat{b}^{\dagger}_{0}+\sum_{\mathbf{k}\neq
0}g_{\mathbf{k}}(t)\hat{b}_{\mathbf{k}}^{{\dagger}}\hat{b}_{-\mathbf{k}}^{{\dagger}}\right]|0\rangle,$
(20)

where $\mathcal{N}(t)$ is a normalization factor, $|0\rangle$ is vacuum of
particles, and $g_{0}$ and $g_{{\bf k}}$ are all variational parameters. This
approach is not restricted to the short time and has been successfully used in
the previous studies of degenerate Bose gas quenched to unitarity theory2 ;
theory5 ; theory18 . The evolution of variational parameters $g_{0}(t)$ and
$g_{\mathbf{k}}(t)$ can be obtained from the Euler-Lagrange equation for the
Lagrangian
$\mathcal{L}=\frac{1}{2}[\langle\Psi(t)|\dot{\Psi}(t)\rangle-\langle\dot{\Psi}(t)|\Psi(t)\rangle]-\langle\Psi(t)|\hat{H}(t)|\Psi(t)\rangle$,
which yields a set of differential equations for $g_{0}$ and $g_{{\bf k}}$.
Since we start with a non-interacting Bose condensate, we take $g_{0}=1$ and
$g_{{\bf k}}=0$ at $t=0$ as the initial conditions for these differential
equations. We can obtain the variational wave function by solving these
equations, and subsequently, we can determine the total energy density with
Eq. 13. The results for the total energy density are shown in Fig. 3 for
different powers and different coefficients. One can see that the short time
behaviors agree very well with that given in Fig. 1(c1) and (c2).

Summary. In summary, we have studied the energy growth rate of degenerate
quantum gas driven by increasing the $s$-wave scattering length $a_{\text{s}}$
from zero, by both analyzing the short time behavior on general situations and
numerically solving a concrete example of interacting bosons. Two main results
are summarized as follows: (i) For energy increasing as $t^{\alpha}$ at the
short time, $\alpha$ cannot be smaller than $1/2$ and $\alpha=1/2$ is achieved
when $a_{\text{s}}(t)$ varies as $\propto\sqrt{t}$. (ii) For energy increasing
as $\sqrt{t}$ at the short time, the fastest energy increasing is achieved
when $a_{\text{s}}(t)=2\sqrt{\hbar t/(\pi m)}$, with a maximum energy growth
given by $4(2+\sqrt{2})\sqrt{\pi\hbar^{3}t/m}g_{2}(0)$ for the spin-$1/2$
Fermi gas and half of that for the spinless Bose gas. This prediction can be
directly verified in cold atom experiments. We emphasize that this maximum
energy growth rate is universal, that is, it is independent of any length or
energy scale. This is because when $a_{\text{s}}$ varies as $\sqrt{t}$, the
entire many-body Schrödinger equation is invariant under a space-time scaling
transformation $t\rightarrow\lambda^{2}t$ and $r\rightarrow\lambda r$. Similar
examples of such scale invariant many-body dynamics have been studied in
Efimovian ; Efimovian2 ; Efimovian3 . Hence, this result ties together the
fastest energy growth with the scaling symmetry, and this is reminiscent of an
equilibrium analogy, where the interaction effect is the strongest at
unitarity where the system is also scale invariant.

Acknowledgment. We thank Peng Zhang for helpful discussions. The project was
supported by NSFC under Grant No. 12022405 (RQ), No. 11774426 (RQ) and No.
11734010 (HZ and RQ), Beijing Outstanding Young Scholar Program (HZ), the
National Key R and D Program of China under Grant No. 2018YFA0306501(RQ), the
Research Funds of Renmin University of China under Grant No. 19XNLG12 (RQ),
the Beijing Natural Science Foundation under Grant No. Z180013 (RQ), and
Program of Shanghai Sailing Program Grant No. 20YF1411600 (ZYS).

## References

  * (1) T. K$\ddot{\text{o}}$hler, K. G$\acute{\text{o}}$ral, and P. S. Julienne, Rev. Mod. Phys. 78, 1311 (2006).
  * (2) C. Chin, R. Grimm, P. S. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * (3) P. Makotyn, C. E. Klauss, D. L. Goldberger, E. A. Cornell, and D. S. Jin, Nat. Phys. 10, 116 (2014).
  * (4) C. Eigen, J. A. Glidden, R. Lopes, E. A. Cornell, R. P. Smith, and Z. Hadzibabic, Nature (London) 563, 221 (2018).
  * (5) A. Behrle, T. Harrison, J. Kombe, K. Gao, M. Link, J.-S. Bernier, C. Kollath and M. Köhl, Nat. Phys. 14, 781 (2018)
  * (6) L. W. Clark, A. Gaj, L. Feng and C. Chin, Nature 551, 356 (2017).
  * (7) L. Feng, J. Hu, L. W. Clark, and C. Chin, Science 363, 521 (2019).
  * (8) J. Hu, L. Feng, Z. Zhang, and C. Chin, Nat. Phys. 15, 785 (2019).
  * (9) X. Yin and L. Radzihovsky, Phys. Rev. A 88, 063611 (2013).
  * (10) A. G. Sykes, J. P. Corson, J. P. D’Incao, A. P. Koller, C. H. Greene, A. M. Rey, K. R. Hazzard, and J. L. Bohn, Phys. Rev. A 89, 021601 (2014).
  * (11) A. Rançon and K. Levin, Phys. Rev. A 90, 021602 (2014).
  * (12) B. Kain and H. Y. Ling, Phys. Rev. A 90, 063626 (2014).
  * (13) J. P. Corson and J. L. Bohn, Phys. Rev. A 91, 013616 (2015)
  * (14) F. Ancilotto, M. Rossi, L. Salasnich, and F. Toigo, Few-Body Syst. 56, 801 (2015).
  * (15) X. Yin and L. Radzihovsky, Phys. Rev. A 93, 033653 (2016).
  * (16) A. Eckardt, Rev. Mod. Phys. 89, 011004 (2017).
  * (17) V. E. Colussi, J. P. Corson, and J. P. D’Incao, Phys. Rev. Lett. 120, 100401 (2018).
  * (18) V. E. Colussi, S. Musolino, and S. J. J. M. F. Kokkelmans, Phys. Rev. A 98, 051601 (2018).
  * (19) M. Van Regemortel, H. Kurkjian, M. Wouters, and I. Carusotto, Phys. Rev. A 98, 053612 (2018).
  * (20) J. P. D’Incao, J. Wang, and V. E. Colussi, Phys. Rev. Lett. 121, 023401 (2018).
  * (21) H. Fu, L. Feng, B. M. Anderson, L. W. Clark, J. Hu, J. W. Andrade, C. Chin, and K. Levin, Phys. Rev. Lett. 121, 243001 (2018).
  * (22) T. Chen and B. Yan, Phys. Rev. A 98, 063615 (2018).
  * (23) S. Musolino, V. E. Colussi, and S. J. J. M. F. Kokkelmans, Phys. Rev. A 100, 013612 (2019).
  * (24) A. Muñoz de las Heras, M. M. Parish, F. M. Marchetti, Phys. Rev. A 99, 023623 (2019).
  * (25) Z. Wu and H. Zhai, Phys. Rev. A 99, 063624 (2019).
  * (26) C. Gao, M. Sun, P. Zhang, and H. Zhai, Phys. Rev. Lett. 124, 040403 (2020).
  * (27) M. Sun, P. Zhang, and H. Zhai, Phys. Rev. Lett. 125, 110404 (2020).
  * (28) Y.-Y. Chen, P. Zhang, W. Zheng, Z. Wu, and H. Zhai, Phys. Rev. A 102, 011301(R) (2020).
  * (29) Y. Cheng and Z. Y. Shi, arXiv:2004.12754
  * (30) C. Lv, R. Zhang, Q. Zhou, Phys. Rev. Lett. 125, 253002 (2020)
  * (31) S. Tan, Ann. Phys. (N.Y.) 323, 2952 (2008).
  * (32) S. Tan, Ann. Phys. (N.Y.) 323, 2971 (2008).
  * (33) S. Tan, Ann. Phys. (N.Y.) 323, 2987 (2008).
  * (34) M. Punk and W. Zwerger, Phys. Rev. Lett. 99, 170404 (2007)
  * (35) G. Baym, C. J. Pethick, Z. Yu and M. W. Zwierlein, Phys. Rev. Lett. 99, 190407 (2007)
  * (36) E. Braaten and L. Platter, Phys. Rev. Lett. 100, 205301 (2008).
  * (37) S. Zhang and A. J. Leggett. Phys. Rev. A 79, 023601 (2009)
  * (38) E. Braaten, D. Kang, and L. Platter, Phys. Rev. Lett. 106, 153005 (2011).
  * (39) Z. Yu, J. H. Thywissen, and S. Zhang, Phys. Rev. Lett. 115, 135304 (2015).
  * (40) See the supplementary material for the detailed proof.
  * (41) S. Deng, Z.-Y. Shi, P. Diao, Q. Yu, H. Zhai, R. Qi, and H. Wu, Science 353, 371 (2016).
  * (42) Z.-Y. Shi, R. Qi, H. Zhai, and Z. Yu, Phys. Rev. A 96, 050702(R) (2017).
  * (43) S. Deng, P. Diao, F. Li, Q. Yu, S. Yu, and H. Wu, Phys. Rev. Lett. 120, 125301 (2018).

## Supplementary material:
Maximum Energy Growth Rate in Dilute Quantum Gases

In this supplementary material, we provide the details of proof for Eq.
(1)-(3) in the main text.

##  I solution of two-body problem

In this section, we solve the time-dependent two-body problem and establish
Eq. (2) and (3). We consider the following time-dependent Schrödinger equation
in the relative coordinate ${\bf r}$ frame

$\displaystyle i\hbar\frac{\partial}{\partial
t}\psi(\mathbf{r},t)=-\frac{\hbar^{2}}{m}\nabla^{2}\psi(\mathbf{r},t)+\frac{4\pi\hbar^{2}a_{\text{s}}(t)}{m}\delta(\mathbf{r})\frac{\partial}{\partial
r}r\psi(\mathbf{r},t),$ (21)

which is Eq. (4) in the main text and we choose the initial state
$\psi(\mathbf{r})=1/\sqrt{V}$. We first define an auxiliary function $\eta(t)$
as

$\displaystyle\eta(t)=-\sqrt{V}4\pi a_{s}(t)\lim_{r\rightarrow
0}\frac{\partial}{\partial r}[r\psi(\mathbf{r},t)].$ (22)

Then Eq. (21) can be rewritten as

$\displaystyle i\hbar\frac{\partial}{\partial
t}\psi(\mathbf{r},t)=-\frac{\hbar^{2}}{m}\nabla^{2}\psi(\mathbf{r},t)-\frac{1}{\sqrt{V}}\frac{\hbar^{2}}{m}\delta(\mathbf{r})\eta(t).$
(23)

Now Eq. (23) can be solved with the standard Green’s function approach and the
solution is given as

$\displaystyle\psi(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\psi_{0}(\mathbf{r},t)+\frac{i}{\sqrt{V}}\frac{\hbar}{m}\int_{0}^{t}G_{0}^{\text{rel}}(\mathbf{r},t-\tau)\eta(\tau)\tau$
(24)

where $G_{0}^{\text{rel}}(\mathbf{r},t-\tau)$ is the noninteracting Green’s
function in the relative coordinate frame given by

$\displaystyle
G_{0}^{\text{rel}}(\mathbf{r},t-\tau)=\left[\frac{m}{i4\pi\hbar(t-\tau)}\right]^{3/2}\exp\left[i\frac{mr^{2}}{4\hbar(t-\tau)}\right],$
(25)

and $\psi_{0}(\mathbf{r},t)$ satisfies the non-interacting Schrödinger
equation

$\displaystyle i\hbar\frac{\partial}{\partial
t}\psi_{0}(\mathbf{r},t)=-\frac{\hbar^{2}}{m}\nabla^{2}\psi_{0}(\mathbf{r},t),~{}~{}~{}~{}~{}~{}\psi_{0}(\mathbf{r},t=0)=\psi(\mathbf{r})$
(26)

where $\psi(\mathbf{r})$ is the initial wave function of $\psi(\mathbf{r},t)$.

For $\psi(\mathbf{r})=1/\sqrt{V}$, we have $\psi_{0}(\mathbf{r},t)=1/\sqrt{V}$
and, based on Eq. (24), we obtain the following asymptotic expansion at
$r\rightarrow 0$

$\displaystyle\psi(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{V}}\left\\{1+\frac{\eta(t)}{4\pi
r}+\left(\frac{m}{\hbar}\right)^{\frac{1}{2}}\frac{1}{8\pi^{3/2}\sqrt{i}}\lim_{r\rightarrow
0^{+}}\left[\int_{0}^{t}\frac{\eta(\tau)\exp\left[i\frac{mr^{2}}{4\hbar(t-\tau)}\right]}{(t-\tau)^{\frac{3}{2}}}d\tau-\int_{-\infty}^{t}\frac{\eta(t)\exp\left[i\frac{mr^{2}}{4\hbar(t-\tau)}\right]}{(t-\tau)^{\frac{3}{2}}}d\tau\right]+\text{O}(r)\right\\},$
(27) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{V}}\left\\{1+\frac{\eta(t)}{4\pi
r}+\left(\frac{m}{\hbar}\right)^{\frac{1}{2}}\frac{1}{8\pi^{3/2}\sqrt{i}}\lim_{\epsilon\rightarrow
0^{+}}\left[\int_{0}^{t-\epsilon}\frac{\eta(\tau)}{(t-\tau)^{\frac{3}{2}}}d\tau-\int_{-\infty}^{t-\epsilon}\frac{\eta(t)}{(t-\tau)^{\frac{3}{2}}}d\tau\right]+\text{O}(r)\right\\},$
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{V}}\left\\{1+\frac{\eta(t)}{4\pi
r}+\left(\frac{m}{\hbar}\right)^{\frac{1}{2}}\frac{1}{8\pi^{3/2}\sqrt{i}}\lim_{\epsilon\rightarrow
0^{+}}\left[\int_{0}^{t-\epsilon}\frac{\eta(\tau)}{(t-\tau)^{\frac{3}{2}}}d\tau-\frac{2\eta(t)}{\sqrt{\epsilon}}\right]+\text{O}(r)\right\\},$
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{V}}\left[1+\frac{\eta(t)}{4\pi
r}+\hat{L}\eta(t)+\text{O}(r)\right],$

where $\hat{L}$ is defined in Eq. (3) in the main text. Then substituting Eq.
(27) into the r.h.s. of Eq. (22), we immediately obtain

$\displaystyle\left[\hat{L}+\frac{1}{4\pi a_{\text{s}}(t)}\right]\eta(t)=-1,$
(28)

which is exactly the Eq. (2) in the main text.

##  II Generalization to many-body wave function

In this section, we generalize the two-body solution to many-body wave
functions and establish Eq. (1) in the main text. To simplify notations, we
set $\hbar=m=1$ in this section. Since we are only interested in the short
time evolution of the short distance behavior of the many-body wave function,
it is convenient to divide the entire $3N$ dimensional space into two regions:
$\mathcal{D}_{\epsilon}$ and $\mathcal{I}_{\epsilon}$.
$\mathcal{D}_{\epsilon}$ denotes the region in which the distance between any
two bosons is larger than $\epsilon$, and $\mathcal{I}_{\epsilon}$ is the
complementary space to $\mathcal{D}_{\epsilon}$. In the short time limit $t\ll
t_{n}$, it is possible to choose an intermediate length scale $\epsilon$ such
that $\sqrt{t}\ll\epsilon\ll 1/k_{n}$ ($t_{n}$ and $k_{n}$ was already defined
in the main text). Such dividing of configuration space will be very useful in
the following proof of this section.

We consider the following many-body Schr$\ddot{\text{o}}$inger equation for
N-interacting bosons with a time-dependent s-wave scattering length
$a_{\text{s}}(t)$

$\displaystyle i\frac{\partial}{\partial
t}\psi(\mathbf{r}_{1},\mathbf{r}_{2},\cdots\mathbf{r}_{N};t)=-\frac{1}{2}\sum_{i=1}^{N}\nabla^{2}_{i}\psi(\mathbf{r}_{1},\mathbf{r}_{2},\cdots\mathbf{r}_{N};t)+4\pi
a_{\text{s}}(t)\sum_{i<j}\delta(\mathbf{r}_{ij})\frac{\partial}{\partial
r_{ij}}r_{ij}\psi(\mathbf{r}_{1},\mathbf{r}_{2},\cdots\mathbf{r}_{N};t)$ (29)

where $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$ and
$\psi(\mathbf{r}_{1},\mathbf{r}_{2},\cdots\mathbf{r}_{N};t)$ is the normalized
N-body wave function satisfying

$\int d^{3}\mathbf{r}_{1}\,\cdots
d^{3}\mathbf{r}_{N}|\psi(\mathbf{r}_{1},\mathbf{r}_{2},\cdots\mathbf{r}_{N};t)|^{2}=1.$
(30)

We first define an auxiliary function

$\displaystyle\Phi(\mathbf{R};\mathbb{Z};t)=-4\pi
a_{\text{s}}(t)\lim_{\mathbf{r}\rightarrow 0}\frac{\partial}{\partial
r}r\psi(\mathbf{R}+\mathbf{r}/2,\mathbf{R}-\mathbf{r}/2,\mathbf{z}_{1},\cdots\mathbf{z}_{N-2};t)$
(31)

where $\mathbb{Z}=\mathbf{z}_{1},\cdots\mathbf{z}_{N-2}$. Then based on Eq.
(29) and (31), one can show that the N-body wave function satisfies the
following asymptotic boundary condition at $r_{ij}\rightarrow 0$ with any
$i<j$:

$\displaystyle\psi(\mathbb{R};t)=\frac{1}{4\pi}\Phi(\mathbf{R}_{ij};\mathbb{R}_{\overline{ij}};t)\left(\frac{1}{r_{ij}}-\frac{1}{a_{\text{s}}(t)}\right)+O(r_{ij}),$
(32)

where $\mathbf{R}_{ij}=(\mathbf{r}_{i}+\mathbf{r}_{j})/2$ and, to simplify
notations, we have defined
$\mathbb{R}=\mathbf{r}_{1},\mathbf{r}_{2},\cdots\mathbf{r}_{N}$ and
$\mathbb{R}_{\overline{ij}}=\mathbf{r}_{1},\cdots,\mathbf{r}_{i-1},\mathbf{r}_{i+1},\cdots,\mathbf{r}_{j-1},\mathbf{r}_{j+1},\cdots\mathbf{r}_{N}$
denotes the coordinates excluding ${\bf r}_{i}$ and ${\bf r}_{j}$.

The function $\Phi(\mathbf{R};\mathbb{Z};t)$ is directly related to the
contact as

$\displaystyle\mathcal{C}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{V}16\pi^{2}N(N-1)\int
d^{3}\mathbf{R}d^{3(N-2)}\mathbb{Z}\left|\frac{1}{4\pi}\Phi(\mathbf{R};\mathbb{Z};t)\right|^{2},$
(33)

where $V$ is the total volume of system and
$d^{3(N-2)}\mathbb{Z}=d^{3}\mathbf{z}_{1}\,\cdots d^{3}\mathbf{z}_{N-2}$.

Now Eq. (29) can be rewritten as

$\displaystyle i\frac{\partial}{\partial
t}\psi(\mathbb{R};t)=-\frac{1}{2}\sum_{i=1}^{N}\nabla^{2}_{i}\psi(\mathbb{R};t)-\sum_{i<j}\delta(\mathbf{r}_{ij})\Phi\left(\frac{\mathbf{r}_{i}+\mathbf{r}_{j}}{2};\mathbb{R}_{\overline{ij}};t\right).$
(34)

This equation can still be solved by standard Green’s function approach and we
obtain the following form of solution

$\displaystyle\psi(\mathbb{R};t)=\psi_{0}(\mathbb{R};t)+\sum_{i<j}\psi_{\text{int}}(\mathbf{r}_{i},\mathbf{r}_{j};\mathbb{R}_{\overline{ij}};t),$
(35)
$\displaystyle\psi_{\text{int}}(\mathbf{a},\mathbf{b};\mathbb{Z};t)=i\int_{0}^{t}d\tau
G_{0}^{\text{rel}}(\mathbf{a}\\!-\\!\mathbf{b},t\\!-\\!\tau)\Phi_{0}^{t-\tau}\left(\frac{\mathbf{a}\\!+\\!\mathbf{b}}{2},\mathbb{Z};\tau\right),$
(36) $\displaystyle\Phi_{0}^{t-\tau}(\mathbf{R};\mathbb{Z};\tau)=\int
d^{3}\mathbf{R}^{\prime}d^{3(N-2)}\mathbb{Z}^{\prime}G_{0}^{\text{cm}}\left(\mathbf{R}-\mathbf{R}^{\prime},t\\!-\\!\tau\right)\frac{\exp\left[i\frac{(\mathbb{Z}-\mathbb{Z}^{\prime})^{2}}{2(t-\tau)}\right]}{[2\pi
i(t-\tau)]^{3(N-2)/2}}\Phi(\mathbf{R}^{\prime};\mathbb{Z}^{\prime};\tau),$
(37)

where we defined the following two Green’s functions

$G_{0}^{\text{cm}}(\mathbf{r},t)=\frac{1}{(\pi
it)^{3/2}}\exp\left(i\frac{r^{2}}{t}\right),~{}~{}~{}~{}~{}~{}G_{0}^{\text{rel}}(\mathbf{r},t)=\frac{1}{(4\pi
it)^{3/2}}\exp\left(i\frac{r^{2}}{4t}\right).$ (38)

The function $\Phi_{0}^{t-\tau}(\mathbf{R};\mathbb{Z};\tau)$ can be seen as
the free expansion with time $t-\tau$ from an ‘initial wave function’
$\Phi(\mathbf{R};\mathbb{Z};\tau)$.

The non-interacting wave function $\psi_{0}(\mathbb{R};t)$ satisfies

$\displaystyle i\frac{\partial}{\partial
t}\psi_{0}(\mathbb{R};t)=-\frac{1}{2}\sum_{i=1}^{N}\nabla^{2}_{i}\psi_{0}(\mathbb{R};t),~{}~{}~{}~{}~{}~{}~{}~{}\psi_{0}(\mathbb{R};t=0)=\psi(\mathbb{R};t=0).$
(39)

It is easy to show that
$\psi_{\text{int}}(\mathbf{a},\mathbf{b};\mathbb{Z};t)$ is only singular at
$\mathbf{a}-\mathbf{b}\rightarrow 0$, and we obtain the following asymptotic
behavior as $r_{ij}\equiv|{\bf r}_{ij}|\rightarrow 0$

$\displaystyle\psi_{\text{int}}(\mathbf{r}_{i},\mathbf{r}_{j};\mathbb{R}_{\overline{ij}};t)=\frac{\Phi(\mathbf{R}_{ij},\mathbb{R}_{\overline{ij}};t)}{4\pi
r_{ij}}+i\text{Z}\int_{0}^{t}d\tau\frac{\Phi_{0}^{t-\tau}\left(\mathbf{R}_{ij};\mathbb{R}_{\overline{ij}};\tau\right)}{[4\pi
i(t-\tau)]^{3/2}}+O(r_{ij}),$ (40)

where we have defined

$\displaystyle\text{Z}\int_{0}^{t}d\tau\frac{\Phi_{0}^{t-\tau}\left(\mathbf{R};\mathbb{Z};\tau\right)}{[4\pi
i(t-\tau)]^{3/2}}=\frac{1}{(4\pi i)^{3/2}}\lim_{\epsilon\rightarrow
0^{+}}\left[\int_{0}^{t-\epsilon}d\tau\frac{\Phi_{0}^{t-\tau}\left(\mathbf{R};\mathbb{Z};\tau\right)}{(t-\tau)^{3/2}}-2\frac{\Phi(\mathbf{R};\mathbb{Z};t)}{\sqrt{\epsilon}}\right].$
(41)

Finally, combining Eq. (31), (35) and (40) we obtain an closed integral
equation for $\Phi(\mathbf{R};\mathbb{Z};t)$

$\displaystyle-\frac{\Phi(\mathbf{R};\mathbb{Z};t)}{4\pi
a(t)}=\psi_{0}(\mathbf{R},\mathbf{R},\mathbb{Z};t)+iZ\int_{0}^{t}d\tau\frac{\Phi_{0}^{t-\tau}\left(\mathbf{R};\mathbb{Z};\tau\right)}{[4\pi
i(t-\tau)]^{3/2}}+2\sum_{i=1}^{N-2}\psi_{\text{int}}(\mathbf{z}_{i},\mathbf{R};\mathbb{\bar{Z}}_{i};t)+\sum_{i<j}^{N-2}\psi_{\text{int}}(\mathbf{z}_{i},\mathbf{z}_{j};\mathbb{\bar{Z}}_{ij};t)$
(42)

where
$\mathbb{\bar{Z}}_{i}=\mathbf{z}_{1},\cdots,\mathbf{z}_{i-1},\mathbf{R},\mathbf{z}_{i+1},\cdots,\mathbf{z}_{N-2}$
and
$\mathbb{\bar{Z}}_{ij}=\mathbf{z}_{1},\cdots,\mathbf{z}_{i-1},\mathbf{R},\mathbf{z}_{i+1},\cdots,\mathbf{z}_{j-1},\mathbf{R},\mathbf{z}_{j+1},\cdots\mathbf{z}_{N-2}$.

In the case of spinless bosons starting from a pure non-interacting BEC state,
we have $\psi(\mathbb{R},t=0)=1/V^{N/2}$ and thus
$\psi_{0}(\mathbb{R},t)=1/V^{N/2}$. Below we will show that the ansatz
$\Phi(\mathbf{R};\mathbb{Z};t)=\eta(t)/V^{N/2}$ provides an asymptotic exact
solution in the region $\mathcal{D}_{\epsilon}$. Substituting this ansatz into
Eq. (42), we obtain

$\displaystyle-\frac{\eta(t)}{4\pi a(t)}-1-\hat{L}\eta(t)$ $\displaystyle=$
$\displaystyle
2\sum_{i=1}^{N-2}\phi(\mathbf{z}_{i}-\mathbf{R},t)+\sum_{i<j}^{N-2}\phi(\mathbf{z}_{i}-\mathbf{z}_{j},t),$
(43) $\displaystyle\phi(\mathbf{r},t)$ $\displaystyle=$ $\displaystyle
i\int_{0}^{t}d\tau G_{0}^{\text{rel}}(\mathbf{r},t\\!-\\!\tau)\eta(\tau).$
(44)

On a first sight, Eq. (43) cannot be full filled because the l.h.s. only
depends on $t$ but the r.h.s. is a function of $t,~{}\mathbb{Z}$ and
$\mathbf{R}$. However, it is important to note that since
$\int\phi(\mathbf{r},t)d^{3}\mathbf{r}=i\int_{0}^{t}\eta(\tau)d\tau$ is
finite, $\phi(\mathbf{r},t)$ must decay very fast as $r\gg\sqrt{t}$. For
example, for $a_{\text{s}}(t)=\beta\sqrt{t}$, we have

$\displaystyle\phi_{\beta}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{1-\frac{2\sqrt{i}}{\sqrt{\pi}\beta}}\left[\text{Erf}\left(\frac{r}{2\sqrt{it}}\right)+\frac{2\sqrt{it}}{\sqrt{\pi}r}\exp\left(i\frac{r^{2}}{4t}\right)-1\right],$
(45)

where $\text{Erf}(z)$ is the error function in the complex plane. In this
case, we have $\phi_{\beta}(\mathbf{r},t)\rightarrow t^{3/2}/r^{3}$ as
$r\gg\sqrt{t}$. For other form of $a_{\text{s}}(t)$, $\phi(\mathbf{r},t)$
decays even faster as $r\gg\sqrt{t}$. On the other hand, in region
$\mathcal{D}_{\epsilon}$, we have $r_{ij}>\epsilon\gg\sqrt{t}$. This means the
r.h.s. of Eq. (43) vanishes in the region $\mathcal{D}_{\epsilon}$ and Eq.
(43) is reduced to Eq. (28) ( or Eq. (2) in the main text). As a result, the
ansatz solution

$\displaystyle\Phi(\mathbf{R};\mathbb{Z};t)=\frac{\eta(t)}{V^{N/2}},$ (46)

with $\eta(t)$ satisfying Eq. (28) becomes asymptotically exact in the region
$\mathcal{D}_{\epsilon}$.

Finally, substituting Eq. (28) into (33) and the contact is given as

$\displaystyle\mathcal{C}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{V}16\pi^{2}N(N-1)\int_{\mathcal{D}_{\epsilon}+\mathcal{I}_{\epsilon}}d^{3}\mathbf{R}d^{3}\mathbf{r}_{3}\,\cdots
d^{3}\mathbf{r}_{N}\left|\frac{1}{4\pi}\Phi(\mathbf{R};\mathbf{r}_{3}\,\cdots\mathbf{r}_{N};t)\right|^{2}$
(47) $\displaystyle=$
$\displaystyle\frac{1}{V}16\pi^{2}N(N-1)\int_{\mathcal{D}_{\epsilon}}d^{3}\mathbf{R}d^{3}\mathbf{r}_{3}\,\cdots
d^{3}\mathbf{r}_{N}\left|\frac{1}{4\pi}\Phi(\mathbf{R};\mathbf{r}_{3}\,\cdots\mathbf{r}_{N};t)\right|^{2}$
(48) $\displaystyle=$
$\displaystyle\frac{1}{V}16\pi^{2}N(N-1)V^{N-1}\left|\frac{1}{4\pi}\frac{\eta(t)}{V^{N/2}}\right|^{2}$
(49) $\displaystyle=$ $\displaystyle n^{2}|\eta(t)|^{2},$ (50)

where $V$ is the total volume of the system. The equality between Eq. (47) and
(48) is due to the fact that the contribution to the integral from the region
$\mathcal{I}_{\epsilon}$ vanishes in the limit $\epsilon\ll 1/k_{n}$. The
relative error caused by omitting the term
$\int_{\mathcal{I}_{\epsilon}}\cdots$ is on the order of $k_{n}\epsilon$.

It is also straightforward to generalize above analysis to an initial state of
arbitrary low energy plane wave state where

$\displaystyle\psi(\mathbf{r}_{1},\cdots\mathbf{r}_{N};t=0)$ $\displaystyle=$
$\displaystyle\frac{1}{N!V^{N/2}}\sum_{P}\exp\left(i\sum_{j=1}^{N}\mathbf{k}_{P_{j}}\cdot\mathbf{r}_{j}\right),$
(51) $\displaystyle\psi_{0}(\mathbf{r}_{1},\cdots\mathbf{r}_{N};t)$
$\displaystyle=$ $\displaystyle
e^{-iE_{0}t}\psi(\mathbf{r}_{1},\cdots\mathbf{r}_{N};0),$ (52) $\displaystyle
E_{0}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\sum_{j=1}^{N}k_{j}^{2},$
(53)

where $P$ runs over all permutations. Here low energy means that
$\max\\{k_{i}\\}\ll 1/\sqrt{t}$. In this case, one can show that

$\displaystyle\Phi(\mathbf{R};\mathbb{Z};t)=\eta(t)\psi_{0}(\mathbf{R},\mathbf{R},\mathbb{Z};t),$
(54)

also provides an asymptotically exact solution in the region
$\mathcal{D}_{\epsilon}$ and the contact is given by

$\displaystyle\mathcal{C}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{V}N(N-1)\int_{\mathcal{D}_{\epsilon}}d^{3}\mathbf{R}d^{3}\mathbf{r}_{3}\,\cdots
d^{3}\mathbf{r}_{N}\left|\psi(\mathbf{R},\mathbf{R},\mathbf{r}_{3}\,\cdots\mathbf{r}_{N};0)\right|^{2}|\eta(t)|^{2}$
(55) $\displaystyle=$ $\displaystyle\frac{1}{V}\int
d^{3}\mathbf{R}\langle\Psi_{0}|\hat{\psi}^{{\dagger}}(\mathbf{R})\hat{\psi}^{{\dagger}}(\mathbf{R})\hat{\psi}(\mathbf{R})\hat{\psi}(\mathbf{R})|\Psi_{0}\rangle|\eta(t)|^{2}$
(56) $\displaystyle=$ $\displaystyle\frac{1}{V}\int
d^{3}\mathbf{R}\langle\Psi_{0}|\hat{\psi}^{{\dagger}}(\mathbf{0})\hat{\psi}^{{\dagger}}(\mathbf{0})\hat{\psi}(\mathbf{0})\hat{\psi}(\mathbf{0})|\Psi_{0}\rangle|\eta(t)|^{2}$
(57) $\displaystyle=$ $\displaystyle g_{2}(0)|\eta(t)|^{2},$ (58)

which is Eq. (1) in the main text.

The above proof can be generalized straightforwardly to spin-1/2 fermions. In
this case, we have

$\displaystyle\psi(\mathbf{r}_{1\uparrow},\cdots\mathbf{r}_{N\uparrow};\mathbf{r}_{1\downarrow},\cdots\mathbf{r}_{M\downarrow};t=0)$
$\displaystyle=$
$\displaystyle\frac{1}{N!M!V^{(N+M)/2}}\sum_{P_{\uparrow},P_{\downarrow}}(-1)^{P_{\uparrow}}(-1)^{P_{\downarrow}}\exp\left(i\sum_{j=1}^{N}\mathbf{k}_{P_{j\uparrow}}\cdot\mathbf{r}_{j\uparrow}+i\sum_{j=1}^{M}\mathbf{k}_{P_{j\downarrow}}\cdot\mathbf{r}_{j\downarrow}\right),$
$\displaystyle\psi_{0}(\mathbf{r}_{1\uparrow},\cdots\mathbf{r}_{N\uparrow};\mathbf{r}_{1\downarrow},\cdots\mathbf{r}_{M\downarrow};t)$
$\displaystyle=$ $\displaystyle
e^{-iE_{0}t}\psi(\mathbf{r}_{1\uparrow},\cdots\mathbf{r}_{N\uparrow};\mathbf{r}_{1\downarrow},\cdots\mathbf{r}_{M\downarrow};0),$
(59) $\displaystyle E_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{j=1}^{N}k_{j\uparrow}^{2}+\frac{1}{2}\sum_{j=1}^{M}k_{j\downarrow}^{2},$
(60)

and following a similar derivation as given above for Bose gas, one can show
that the contact is given as

$\displaystyle\mathcal{C}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{V}NM\int_{\mathcal{D}_{\epsilon}}d^{3}\mathbf{R}\prod_{\sigma=\uparrow,\downarrow}d^{3}\mathbf{r}_{2\sigma}\,\cdots
d^{3}\mathbf{r}_{N\sigma}|\psi(\mathbf{R},\mathbf{r}_{2\uparrow},\cdots\mathbf{r}_{N\uparrow};\mathbf{R},\mathbf{r}_{2\downarrow},\cdots\mathbf{r}_{M\downarrow};0)|^{2}|\eta(t)|^{2}$
(61) $\displaystyle=$ $\displaystyle\frac{1}{V}\int
d^{3}\mathbf{R}\langle\Psi_{0}|\hat{\psi}^{{\dagger}}_{\uparrow}(\mathbf{R})\hat{\psi}^{{\dagger}}_{\downarrow}(\mathbf{R})\hat{\psi}_{\downarrow}(\mathbf{R})\hat{\psi}_{\uparrow}(\mathbf{R})|\Psi_{0}\rangle|\eta(t)|^{2}$
(62) $\displaystyle=$ $\displaystyle\frac{1}{V}\int
d^{3}\mathbf{R}\langle\Psi_{0}|\hat{\psi}^{{\dagger}}_{\uparrow}(\mathbf{0})\hat{\psi}_{\downarrow}^{{\dagger}}(\mathbf{0})\hat{\psi}_{\downarrow}(\mathbf{0})\hat{\psi}_{\uparrow}(\mathbf{0})|\Psi_{0}\rangle|\eta(t)|^{2}$
(63) $\displaystyle=$ $\displaystyle g_{2}(0)|\eta(t)|^{2},$ (64)

which also gives Eq. (1) in the main text. Here for spin-1/2 fermions,
$\mathcal{D}_{\epsilon}$ is defined as the region in which the distance
between any $\uparrow$-particle and $\downarrow$-particle is larger than
$\epsilon$.