# Quantum speed limit and non-Markovianity in structured environments

Maryam Hadipour  0000-0002-6573-9960 Faculty of Physics, Urmia University of
Technology, Urmia, Iran  Soroush Haseli  0000-0003-1031-4815
<EMAIL_ADDRESS>Faculty of Physics, Urmia University of Technology,
Urmia, Iran  School of Physics, Institute for Research in Fundamental Sciences
(IPM), P.O. Box 19395-5531, Tehran, Iran  Saeed Haddadi  0000-0002-1596-0763
<EMAIL_ADDRESS>Faculty of Physics, Semnan University, P.O. Box
35195-363, Semnan, Iran

###### Abstract

We investigate the relationship between quantum speed limit time and the non-
Markovianity of an atom in structured environments. We show that there exists
an inverse relation between them, which means that the non-Markovian feature
of the quantum process leads to speedup of the process. Our results might shed
light on the relationship between the speedup of quantum evolution and the
backflow of information from the environment to the system.

Quantum speed limit; non-Markovianity; quantum speedup; structured
environments

##  I Introduction

How to accelerate the evolution of quantum systems is one of the most
important questions in quantum mechanics theory. Based on this theory, quantum
speed limit (QSL) quantifies the minimum evolution time for a quantum system
to pass through a predetermined distance [1, 2, 3, 4, 5, 6, 7, 8]. It refers
to the shortest possible time required for a quantum system to evolve from an
arbitrary initial state to a target final state [9, 10, 11, 12, 13, 14, 15,
16]. QSL time and its applications have been extensively studied
experimentally and theoretically [17, 18, 19, 20, 21]. Moreover, QSL time
plays a particularly important role in the development of quantum information
protection and quantum optimal theory as well [22, 23, 24, 25]. In recent
years, there has been a particular focus on utilizing the advantages of QSL
time in quantum batteries [26, 27, 28, 29]. A deeper knowledge of QSL has led
to the expansion of the speed limit to classical systems [30, 31, 32]. There
is a strong correlation between the vanishing QSL times and the emerging
classical behavior due to the reduced uncertainty in quantum observables [33].
In Ref. [34], the application of QSL time to quantum resource theories has
been investigated. In the beginning, only isolated systems were studied in
accordance with the concept of QSL time [35, 36, 37].

Initially, Mandelstam and Tamm (MT) introduced the minimal time for evolution
from an initial to an orthogonal state for a closed quantum system as
$\tau\geqslant\frac{\pi\hbar}{2\Delta E}$, where $(\Delta
E)^{2}=\left\langle\psi\left|H^{2}\right|\psi\right\rangle-(\langle\psi|H|\psi\rangle)^{2}$
is the energy variance of the system [35]. This bound is known as MT bound.
Another bound related to the mean energy $E=\langle\psi|H|\psi\rangle$ was
then derived by Margolus and Levitin (ML) as $\tau\geq\frac{\pi\hbar}{2E}$,
where it is assumed that the ground state energy $E_{0}$ is equal to zero
[36]. This bound is known as ML bound. It is possible to implement a unified
QSL bound for the unitary evolution of closed quantum systems by combining
both ML and MT bounds as [37]

$\tau\geq\max\left\\{\frac{\pi\hbar}{2\Delta E},\frac{\pi\hbar}{2E}\right\\}.$
(1)

The ML and MT bounds have been expanded and improved using a wide range of
metrics, including trace distance, Bures angle, relative purity, and many
other metrics [38, 39, 40, 41]. The distance among generalized Bloch vectors
for both unitary and non-unitary processes can also be used to obtain tight
QSL bounds [42, 43, 44, 45, 46, 47].

It should be noted that in the real world, quantum systems are open and their
evolution is non-unitary [48]. A variety of approaches can be used to analyze
the behavior of open systems, including the master equation formalism and the
quantum trajectory method. The master equation is a central tool for studying
the dynamics of open systems. It provides a mathematical description of how
the system’s density matrix changes over time. The nature of an interaction
between a quantum system and its environment can influence its dynamics and
evolution. The system-environment interaction is most often analyzed from both
the non-Markovian and Markovian perspectives [48]. In a Markovian regime, the
evolution of the system is purely determined by its current state. In this
case, the system-environment interactions occur over a much shorter timescale
than the significant system evolution. However, in a non-Markovian regime,
system-environment interactions are longer in time than the intrinsic
evolution of the system [49, 50, 51, 52, 53, 54, 55, 56]. An important feature
of non-Markovian regimes is that the interaction behavior of the system with
its environment in the past determines the evolution of the system in the
future. Due to the unavoidable interaction between a quantum system and its
environment, the decoherence effect caused by this coupling can have a
significant effect on quantum statistics.

The ratio between the QSL time and the actual time quantifies the possibility
of speeding up the quantum evolution. QSL bounds are saturated with no
potential for speedup when the ratio equals one. However, the QSL bound is
unsaturated if the ratio is less than one, and there exists potential for
speeding up the evolution of the dynamical system. In Ref. [42], it has been
shown that the non-Markovian effect can accelerate quantum evolution [42].
Both theoretically and experimentally, non-Markovianity regulation has been
proposed as a mechanism for regulating dynamical speedup. Moreover, it has
been shown that non-Markovian effects induce unsaturated QSL bounds in open
dynamics [57, 58, 59, 60, 61, 62]. It is important to note that dynamical
speedup in this way depends very much on the feedback from the environment.
Therefore, it will be essential to investigate other sources that can generate
unsaturated QSL bounds besides non-Markovianity as a source of unsaturated QSL
bounds [63, 64, 65, 66, 67].

Motivated by this, we investigate the possibility of quantum speedup using a
convenient cavity-based engineered environment. We show how the Markovian
environments consisting of two cavities can accelerate the dynamics of an
artificial atom interacting with a pseudomode. In the considered model, the
coupling between the pseudomodes, atom-pseudomode coupling, and detuning
between atom and pseudomode can have significant impacts on QSL time and non-
Markovianity.

The work is organized as follows. First, we present a physical model and its
solution in Sec. II. Then, we discuss the QSL time and non-Markovianity for
the considered model in Sec. III. Next, we turn to the results and discussion
in Sec. IV. Finally, we provide a brief conclusion of the most important
results in Sec. V.

##  II Model and pseudomode method

The considered model consists of a two-level atom interacting with a
Lorentzian structured environment. Figure 1 shows the schematic representation
of the scheme. For this model, the Hamiltonian can be expressed as follows
[68, 69]

$\displaystyle\hat{H}=$
$\displaystyle\frac{\omega_{s}}{2}\hat{\sigma}_{z}+\omega_{1}\hat{a}_{1}^{\dagger}\hat{a}_{1}+\omega_{2}\hat{a}_{2}^{\dagger}\hat{a}_{2}$
(2)
$\displaystyle+k\left(\hat{a}_{1}^{\dagger}\hat{\sigma}_{-}+\hat{a}_{1}\hat{\sigma}_{+}\right)+J\left(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{1}\hat{a}_{2}^{\dagger}\right),$

where $\hat{\sigma}_{z}=|e\rangle\langle e|-|g\rangle\langle g|$ is the
ordinary Pauli operator in $z$-direction of atom, $\omega_{s}$ is the
transition frequency of the atom, $\sigma_{+}$($\sigma_{-}$) is raising
(lowering) operator of the atom, and $\omega_{1}$ and $\omega_{2}$ are the
frequencies of the modes inside the Markovian environments shown by cavity $1$
and cavity $2$, respectively. Let us take $\omega_{1}=\omega_{2}=\omega$ and
set $\omega_{0}=\omega+\delta$, where $\delta$ is atom-pseudomode detuning.
Besides $\hat{a}_{1}$($\hat{a}^{{\dagger}}_{1}$) and
$\hat{a}_{2}$($\hat{a}^{{\dagger}}_{2}$) are the annihilation (creation)
operators.

In what follows, we study the dynamics using the pseudomode theory [70, 71,
72, 73]. This approach relies on the relationship between the atom dynamics
and the shape of the spectral distribution of the environment. The
significance of precise approaches is of current interest due to the
experimental realization of quantum systems. As shown in Fig. 1, the
environment can be depicted by two nondegenerate pseudomodes that leak into
the Markovian environments (cavities $1$ and $2$) with dissipation rates
$\Gamma_{1}$ and $\Gamma_{2}$. Here, the two-level atom only interacts with
the first pseudomode Pm1 (the strength of the coupling $k$) which is in turn
coupled to the second pseudomode Pm2 (the strength of the coupling $J$).

The dynamics of the whole system for a Lorentzian spectral distribution can be
obtained by the following exact pseudomode master equation

$\displaystyle\dot{\rho}(t)=$ $\displaystyle-i[\hat{H},\rho(t)]$ (3)
$\displaystyle-\sum_{n=1,2}\frac{\Gamma_{n}}{2}\left[a_{n}^{\dagger}a_{n}\rho(t)-2a_{n}\rho(t)a_{n}^{\dagger}+\rho(t)a_{n}^{\dagger}a_{n}\right].$

Figure 1: Schematic representation of the model. Two-level atom interacts with
pseudomode $1$ (coupling strength $k$) which leaks into cavity $1$ via
dissipation rate $\Gamma_{1}$, and pseudomode $1$ is coupled to pseudomode $2$
with coupling strength $J$ which the latter leaks into cavity $2$ via
dissipation rate $\Gamma_{2}$.

The above master equation represents the coherent interaction between the two-
level atom and two pseudomodes in the presence of the pseudomodes decay caused
by the interaction with Markovian environments.

Here, it is assumed that the atom is initially in an excited state
$|e_{a}\rangle$ while both pseudomodes Pm${}_{\textmd{1}}$ and
Pm${}_{\textmd{2}}$ are in the ground state
$|g_{\textmd{Pm}_{1}}g_{\textmd{Pm}_{2}}\rangle$, that is, the initial state
of the overall system is assumed to be
$\rho(0)=|e_{a}g_{\textmd{Pm}_{1}}g_{\textmd{Pm}_{2}}\rangle\langle
e_{a}g_{\textmd{Pm}_{1}}g_{\textmd{Pm}_{2}}|$. Considering that there is only
one excitation in the whole system including a two-level atom and two
cavities, the state of the system at time $t$ can be written as

$\displaystyle\rho(t)=$
$\displaystyle(1-\eta(t))|\psi(t)\rangle\langle\psi(t)|$
$\displaystyle+\eta(t)|g_{a}g_{\textmd{Pm}_{1}}g_{\textmd{Pm}_{2}}\rangle\langle
g_{a}g_{\textmd{Pm}_{1}}g_{\textmd{Pm}_{2}}|,$ (4)

where $0\leq\eta(t)\leq 1$ with $\eta(0)=0$, and we have

$\displaystyle|\psi(t)\rangle=$ $\displaystyle
c_{1}(t)|g_{a}e_{\textmd{Pm}_{1}}g_{\textmd{Pm}_{2}}\rangle+c_{2}(t)|g_{a}g_{\textmd{Pm}_{1}}e_{\textmd{Pm}_{2}}\rangle$
$\displaystyle+h(t)|e_{a}g_{\textmd{Pm}_{1}}g_{\textmd{Pm}_{2}}\rangle,$ (5)

with $h(0)=1$ and $c_{1}(0)=c_{2}(0)=0$. For convenience, it is appropriate to
consider the following unnormalized state as

$\displaystyle|\widetilde{\psi}(t)\rangle\equiv$
$\displaystyle\sqrt{1-\eta(t)}|\psi(t)\rangle$ (6) $\displaystyle=$
$\displaystyle\tilde{c}_{1}(t)|g_{a}e_{\textmd{Pm}_{1}}g_{\textmd{Pm}_{2}}\rangle+\tilde{c}_{2}(t)|g_{a}g_{\textmd{Pm}_{1}}e_{\textmd{Pm}_{2}}\rangle$
$\displaystyle+\tilde{h}(t)|e_{a}g_{\textmd{Pm}_{1}}g_{\textmd{Pm}_{2}}\rangle,$

where $\widetilde{h}(t)=\sqrt{1-\eta(t)}h(t)$ is the probability amplitude
that the two-level atom is in an excited state while
$\tilde{c}_{n}(t)\equiv\sqrt{1-\eta(t)}c_{n}(t)$ (with $n=1,2$) is the
probability amplitude that the pseudomodes being in their excited states.
Using the above unnormalized state, the density matrix of the total system can
be rewritten as

$\rho(t)=|\widetilde{\psi}(t)\rangle\langle\widetilde{\psi}(t)|+\eta(t)|g_{a}g_{\textmd{Pm}_{1}}g_{\textmd{Pm}_{2}}\rangle\langle
g_{a}g_{\textmd{Pm}_{1}}g_{\textmd{Pm}_{2}}|.$ (7)

The time-dependent amplitudes of $\tilde{h}(t)$, $\tilde{c}_{1}(t)$, and
$\tilde{c}_{2}(t)$ can be obtained by solving the following coupled
differential equations

$\displaystyle i\frac{d\widetilde{h}(t)}{dt}$
$\displaystyle=(\omega+\delta)\widetilde{h}(t)+k\tilde{c}_{1}(t),$ (8)
$\displaystyle i\frac{d\tilde{c}_{1}(t)}{dt}$
$\displaystyle=\left(\omega-\frac{i}{2}\Gamma_{1}\right)\tilde{c}_{1}(t)+k\widetilde{h}(t)+J\tilde{c}_{2}(t),$
$\displaystyle i\frac{d\tilde{c}_{2}(t)}{dt}$
$\displaystyle=\left(\omega-\frac{i}{2}\Gamma_{2}\right)\tilde{c}_{2}(t)+J\tilde{c}_{1}(t).$

In the above set of coupled differential equations, the coefficients
$\widetilde{h}(t)$, $\tilde{c}_{1}(t)$, and $\tilde{c}_{2}(t)$ can be obtained
using the standard Laplace transformations. By knowing these, the time-evolved
density matrix of the two-level atom can be obtained by taking the partial
trace over two pseudomodes, which we will provide in the next section.

##  III QSL time and non-Markovianity

###  III.1 QSL time

We first deal with QSL time using the method introduced in Ref. [43] for open
quantum systems. A comprehensive bound of the QSL time for any arbitrary
initial state can be defined as

$\tau_{\textmd{QSL}}=\max\left\\{\frac{1}{\Lambda_{\tau}^{op}},\frac{1}{\Lambda_{\tau}^{tr}},\frac{1}{\Lambda_{\tau}^{hs}}\right\\}\sin^{2}[\Theta(\rho_{0},\rho_{\tau})]\operatorname{tr}\left[\rho_{0}^{2}\right],$
(9)

where
$\Theta(\rho_{0},\rho_{t})=\arccos\left(\sqrt{\frac{\operatorname{tr}[\rho_{0}\rho_{t}]}{\operatorname{tr}\left[\rho_{0}^{2}\right]}}\right)$
is a function of relative purity between initial state $\rho_{0}$ and the
state at arbitrary time $\rho_{t}$ that determined by the time-dependent non-
unitary master equation $\dot{\rho}_{t}=\mathcal{L}_{t}(\rho_{t})$ and

$\displaystyle\Lambda_{\tau}^{op}=\frac{1}{\tau}\int_{0}^{\tau}dt\|\mathcal{L}_{t}(\rho_{t})\|_{op},$
(10)
$\displaystyle\Lambda_{\tau}^{tr}=\frac{1}{\tau}\int_{0}^{\tau}dt\|\mathcal{L}_{t}(\rho_{t})\|_{tr},$
$\displaystyle\Lambda_{\tau}^{hs}=\frac{1}{\tau}\int_{0}^{\tau}dt\|\mathcal{L}_{t}(\rho_{t})\|_{hs},$

where $\|\mathcal{L}_{t}(\rho_{t})\|_{op}=\lambda_{1}$,
$\|\mathcal{L}_{t}(\rho_{t})\|_{tr}=\sum_{i}\lambda_{i}$ and
$\|\mathcal{L}_{t}(\rho_{t})\|_{hs}=\sqrt{\sum_{i}\lambda_{i}^{2}}$ are
operator norm, trace norm and Hilbert-Schmidth norm, respectively. Also,
$\lambda_{i}$’s are singular values of $\mathcal{L}_{t}(\rho_{t})$ and
$\lambda_{1}$ is the largest singular value of $\mathcal{L}_{t}(\rho_{t})$.
When the denominator of the fraction is $\Lambda_{\tau}^{op}$ and
$\Lambda_{\tau}^{tr}$, we have a generalized ML type of QSL bound for open
quantum systems, whereas when it is $\Lambda_{\tau}^{hs}$, we have MT type
bound on the QSL time. According to Ref. [43], the ML type bound based on the
operator norm provides the tightest QSL bound.

Using the recently achieved bound for QSL time and a given actual driving time
$\tau$, one can estimate the intrinsic speed of a dynamical evolution. In the
case of $\tau_{\textmd{QSL}}=\tau$, the further speedup is not possible and
the speedup is not observed. In contrast, when $\tau_{\textmd{QSL}}<\tau$,
this indicates the possibility of quantum dynamics speeding up. Based on the
assumption that the two-level atom is initially in the excited state as
mentioned in the previous section, the time-evolved density matrix of the atom
takes the following form

$\rho_{t}=|\widetilde{h}(t)|^{2}|e\rangle\langle
e|+(1-|\widetilde{h}(t)|^{2})|g\rangle\langle g|.$ (11)

So, QSL time based on the tightest bound is obtained as

$\tau_{\textmd{QSL}}=\frac{1-|\widetilde{h}(\tau)|^{2}}{\frac{1}{\tau}\int_{0}^{\tau}|\partial_{t}|\widetilde{h}(t)|^{2}|dt}.$
(12)

###  III.2 Non-Markovianity

During evolution, the backflow of information from the environment to the
system is responsible for the non-Markovian feature of the process. For
quantification, the degree of non-Markovianity can be defined as [51]

$\mathcal{N}=\max_{\rho_{1}(0),~{}\rho_{2}(0)}\int_{\sigma>0}\sigma\left[t,~{}\rho_{1}(0),~{}\rho_{2}(0)\right]dt,$
(13)

in which

$\sigma\left[t,~{}\rho_{1}(0),~{}\rho_{2}(0)\right]=\frac{d}{dt}D\left[\rho_{1}(t),~{}\rho_{2}(t)\right],$
(14)

where
$D\left[\rho_{1}(t),~{}\rho_{2}(t)\right]=\frac{1}{2}\operatorname{Tr}\left|\rho_{1}(t)-\rho_{2}(t)\right|$
is the trace distance of a pair of states with the trace norm of an operator
as $|\hat{A}|=\sqrt{\hat{A}\hat{A}^{\dagger}}$.

Note that when $\sigma\left[t,~{}\rho_{1}(0),~{}\rho_{2}(0)\right]>0$, we have
the backflow of information from the environment to the system and so, the
evolution is non-Markovian $\mathcal{N}\neq 0$. However, if
$\sigma\left[t,~{}\rho_{1}(0),~{}\rho_{2}(0)\right]<0$ then $\mathcal{N}=0$,
meaning that the evolution is Markovian, information irreversibly flows from
system to environment, and backflow does not occur. As a result, non-
Markovianity indicates the total backflow of information from the environment
and the system.

In ref. [51], it has been proved that the optimal pair of initial states that
maximize the degree of non-Markovianity $\mathcal{N}$ are
$\rho_{1}(0)=|e\rangle\langle e|$ and $\rho_{2}(0)=|g\rangle\langle g|$. Using
the findings of the previous section and taking into account the mentioned
optimal initial states, one can find
$\mathcal{D}\left[\rho_{1}(t),\rho_{2}(t)\right]=|\widetilde{h}(t)|^{2}$.
Therefore, the degree of non-Markovianity is obtained as

$\mathcal{N}=\frac{1}{2}\left[\int_{0}^{\tau}\partial_{t}|\widetilde{h}(t)|^{2}dt+|\widetilde{h}(\tau)|^{2}-1\right].$
(15)

##  IV Results and discussion

Now everything is ready to study the QSL time and non-Markovianity based on
the analytical expressions (12) and (15) obtained in the previous section.

In Fig. 2, the non-Markovianity (15) has been plotted as functions of atom-
pseudomode coupling $k/\Gamma_{1}$ and pseudomode-pseudomode coupling
$J/\Gamma_{1}$. From this figure, it can be seen that in the absence of
pseudomode coupling strength ($J=0$), the evolution is Markovian in a weak
coupling regime $k/\Gamma_{1}<0.25$ and non-Markovian in a strong coupling
regime $k/\Gamma_{1}>0.25$. In general, the non-Markovian feature of evolution
is revealed for strong atom-pseudomode coupling for any pseudomode-pseudomode
coupling strength. Therefore, Fig. 2 clearly shows the regions in which the
Markovian and non-Markovian dynamics take place.

Interestingly, by employing Eqs. (12) and (15), it is possible to connect QSL
time and degree of non-Markovianity as

$\tau_{\mathrm{QSL}}=\tau\times\left(\frac{2\mathcal{N}}{1-|\widetilde{h}(\tau)|^{2}}+1\right)^{-1}.$
(16)

From Eq. (16), if the evolution is Markovian ($\mathcal{N}=0$), then the
actual evolution time is equal to the QSL time, so we have a saturated bound
for the QSL time. On the other hand, in a situation where the evolution is
non-Markovian ($\mathcal{N}\neq 0$), the QSL time is shorter than the actual
evolution time.

Figure 2: The degree of non-Markovianity as functions of atom-pseudomode
coupling $k/\Gamma_{1}$ and pseudomode-pseudomode coupling $J/\Gamma_{1}$ with
$\Gamma_{2}/\Gamma_{1}=0.5$ and $\delta/\Gamma_{1}=0$.  Figure 3: (a) QSL time
and (b) non-Markovianity as a function of $\Gamma_{1}\tau$ in a weak coupling
regime ($k/\Gamma_{1}=0.1$) for different values of $J/\Gamma_{1}$ with
$\Gamma_{2}/\Gamma_{1}=0.5$ and $\delta/\Gamma_{1}=0$.

In Fig. 3, the QSL time and non-Markovianity have been plotted as a function
of actual evolution time $\tau$ in a weak coupling regime with
$k/\Gamma_{1}=0.1$ for different values of pseudomode-pseudomode coupling
$J/\Gamma_{1}$. From Fig. 3(a), one can see that the QSL time decreases with
increasing $J/\Gamma_{1}$. Besides, Fig. 3(b) shows that in this weak coupling
regime, the non-Markovianity feature of the evolution is revealed by
increasing $J/\Gamma_{1}$. Thus, we conclude that the QSL time decreases with
the emergence of the non-Markovian feature of quantum evolution. Notably, one
can also see that in Markovian evolution (red solid line), the QSL time is
equal to the actual evolution time ($\tau_{\mathrm{QSL}}=\tau$) because
$\mathcal{N}=0$, whose root is in Eq. (16). In other words, it can be said
that the reduction of QSL time originates from the non-Markovian feature of
the quantum evolution.

Figure 4: (a) QSL time and (b) non-Markovianity as a function of
$\Gamma_{1}\tau$ in a strong coupling regime ($k/\Gamma_{1}=0.75$) for
different values of $J/\Gamma_{1}$ with $\Gamma_{2}/\Gamma_{1}=0.5$ and
$\delta/\Gamma_{1}=0$.

Fig. 4 shows the QSL time and non-Markovianity as a function of actual
evolution time $\tau$ in a strong coupling regime with $k/\Gamma_{1}=0.75$ for
different values of $J/\Gamma_{1}$. It seems that the results are similar to
what we observed in the weak coupling regime, however, the non-Markovian
feature of evolution can be observed in this case even when $J/\Gamma_{1}$ is
zero, which is in good agreement with the results of Fig. 2.

Figure 5: (a) QSL time and (b) non-Markovianity as a function of
$k/\Gamma_{1}$ for zero and nonzero values of $\delta/\Gamma_{1}$ with
$\Gamma_{2}/\Gamma_{1}=0.5$, $J/\Gamma_{1}=1$, and $\tau=5$.

In Fig. 5(a), the QSL time has been plotted as a function of $k/\Gamma_{1}$
for both resonant ($\delta/\Gamma_{1}=0$) and off-resonant
($\delta/\Gamma_{1}=5$) cases. The QSL time decreases with increasing
$k/\Gamma_{1}$ for both resonant and off-resonant cases. Notice that the QSL
time in the off-resonant case is shorter than in the resonant case. In Fig.
5(b), we plot the non-Markovianity versus $k/\Gamma_{1}$. From this plot, we
observe that the degree of non-Markovianity increases by the strengthening of
the coupling between the atom and the first pseudomode for both resonant and
off-resonant cases. By comparing Figs. 5(a) and 5(b), it can be concluded that
there exists an inverse relation between QSL time and non-Markovianity.

Figure 6: (a) QSL time and (b) non-Markovianity as a function of
$J/\Gamma_{1}$ for two different values of $k/\Gamma_{1}$ with
$\Gamma_{2}/\Gamma_{1}=0.5$, $\delta/\Gamma_{1}=0$, and $\tau=5$.

The QSL time and non-Markovianity have been plotted as a function of
$J/\Gamma_{1}$ for two different values of $k/\Gamma_{1}$ in Fig. 6. We see
that in a strong coupling regime ($k/\Gamma_{1}=0.75$), the QSL time is
shorter than that of the weak coupling regime ($k/\Gamma_{1}=0.1$). On the
other hand, the non-Markovianity gets its maximum value in a strong coupling
regime. Again, we find that there is an inverse relation between non-
Markovianity and QSL time. These results are in good agreement with the
previous figures.

Figure 7: (a) QSL time and (b) non-Markovianity as a function of
$\delta/\Gamma_{1}$ for two different values of $k/\Gamma_{1}$ with
$\Gamma_{2}/\Gamma_{1}=0.5$, $J/\Gamma_{1}=1$, and $\tau=5$.

Finally, the QSL time and non-Markovianity have been illustrated as a function
of the detuning parameter $\delta/\Gamma_{1}$ for two different values of
$k/\Gamma_{1}$ in Fig. 7. It is worth noting that the results are similar to
the findings shown in the previous plots. Remarkably, the existence of an
inverse relation between non-Markovianity and QSL time is also confirmed in
this case.

Note that the ratio between the QSL time, $\tau_{\mathrm{QSL}}$, and the
actual time, $\tau$, quantifies the possibility of speeding up the quantum
evolution. QSL bounds are saturated with no potential for speedup when
$\tau_{\mathrm{QSL}}/\tau=1$. However, the QSL bound is unsaturated if
$\tau_{\mathrm{QSL}}/\tau<1$, and there is potential for speeding up the
evolution of the dynamical system. Our results show that
$\tau_{\mathrm{QSL}}/\tau\leq 1$, indicating that both of these scenarios
occurred in our study.

##  V Conclusion and outlook

Using a cavity-based engineered environment, we explored the possibility of
manipulating the QSL time. In particular, we have observed that by
constructing the Markovian environments consisting of two cavities, one can
switch between quantum speed-down and speedup regimes for the dynamics of an
artificial atom interacting with a pseudomode. Based on memory effects, the
dynamics of the open systems are divided into two categories: memory-less
(Markovian) and with-memory (non-Markovian) dynamics. Notably, we observed
that there exists an inverse relation between the QSL time and non-
Markovianity. In other words, the non-Markovianity feature of quantum
evolution can lead to quantum speedup in the dynamics. Therefore, this feature
leads to an intrinsic interest in the context of controlling the memory
effects that are associated with open systems. The parameters included atom-
pseudomode coupling, pseudomode-pseudomode coupling, and detuning between atom
and pseudomode were used to manipulate the non-Markovianity. We found that in
a strong coupling regime, for both atom-pseudomode and pseudomode-pseudomode
couplings, the evolution is non-Markovian and the quantum speedup is revealed.
As a result of this study, new insights are gained into the control of QSL
time for open systems, which opens the door to further experimental
development towards the development of quantum thermodynamic devices and
quantum computers.

## CRediT authorship contribution statement

Maryam Hadipour: Writing–original draft, Investigation. Soroush Haseli:
Writing–original draft, Writing–review & editing, Supervision, Software,
Investigation, Conceptualization. Saeed Haddadi: Writing–review & editing,
Investigation, Formal analysis.

## Acknowledgements

S. Haddadi was supported by Semnan University under Contract No. 21270.

## Disclosures

The authors declare that they have no known competing financial interests.

## Data availability

No datasets were generated or analyzed during the current study.

## References

  * [1] Deffner S and Campbell S 2017 J. Phys. A: Math. Theor. 50, 453001.
  * [2] Hilgevoord J 2002 Am. J. Phys. 70, 301.
  * [3] Deffner S 2017 New J. Phys. 19, 103018.
  * [4] Pfeifer P 1993 Phys. Rev. Lett. 70, 3365.
  * [5] Bukov M, Sels D and Polkovnikov A 2019 Phys. Rev. X 9, 011034.
  * [6] Funo K, Shiraishi N and Saito K 2019 New J. Phys. 21, 013006.
  * [7] Hegerfeldt G C 2013 Phys. Rev. Lett. 111, 260501.
  * [8] Garcia-Pintos L P and del Campo A 2019 New J. Phys. 21, 033012.
  * [9] Garcia-Pintos L P and del Campo A 2019 New J. Phys. 21, 033012.
  * [10] Kobe D H and Aguilera-Navarro V C 1994 Phys. Rev. A 50, 933.
  * [11] Jones P J and Kok P 2010 Phys. Rev. A 82, 022107.
  * [12] Xu Z-Y 2016 New J. Phys. 18, 073005.
  * [13] Deffner S and Lutz E 2013 J. Phys. A: Math. Theor. 46, 335302.
  * [14] Shao Y, Liu B, Zhang M, Yuan H and Liu J 2020 Phys. Rev. Res. 2, 023299.
  * [15] Russell B and Stepney S 2014 Phys. Rev. A 90, 012303.
  * [16] Hu X, Sun S and Zheng Y 2020 Phys. Rev. A 101, 042107.
  * [17] Cheneau M et al 2012 Nature 481, 484.
  * [18] del Campo A 2021 Phys. Rev. Lett. 126, 180603.
  * [19] Taddei M M, Escher B M, Davidovich L and de Matos Filho R L 2013 Phys. Rev. Lett. 110, 050402.
  * [20] Escher B M, de Matos Filho R L and Davidovich L 2011 Nat. Phys. 7, 406.
  * [21] Lam M R et al 2021 Phys. Rev. X 11, 011035.
  * [22] Poggi P M 2019 Phys. Rev. A 99, 042116.
  * [23] Caneva T, Murphy M, Calarco T, Fazio R, Montangero S, Giovannetti V and Santoro G E 2009 Phys. Rev. Lett. 103, 240501.
  * [24] Zhang Y-J, Wei H, Yan W-B, Man Z-X, Xia Y-J and Fan H 2021 New J. Phys. 23, 113004.
  * [25] Marvian I and Lidar D A 2015 Phys. Rev. Lett. 115, 210402.
  * [26] Bai S Y and An J H 2020 Phys. Rev. A 102, 060201.
  * [27] Campaioli F, Pollock F A, Binder F C, Celeri L, Goold J, Vinjanampathy S and Modi K 2017 Phys. Rev. Lett. 118, 150601.
  * [28] Hovhannisyan K V, Perarnau-Llobet M, Huber M and Acın A 2013 Phys. Rev. Lett. 111, 240401.
  * [29] Binder F C, Vinjanampathy S, Modi K and Goold J 2015 New J. Phys. 17, 075015.
  * [30] Okuyama M and Ohzeki M 2018 Phys. Rev. Lett. 120, 070402.
  * [31] Shanahan B, Chenu A, Margolus N and del Campo A 2018 Phys. Rev. Lett. 120, 070401.
  * [32] Bolonek-Lason K, Gonera J and Kosinski P 2021 Quantum 5, 482.
  * [33] Poggi P M, Campbell S and Deffner S 2021 PRX Quantum 2, 040349.
  * [34] Campaioli F, Yu C-s, Pollock F A and Modi K 2022 New J. Phys. 24, 065001.
  * [35] Mandelstam L and Tamm I 1991 The uncertainty relation between energy and time in non-relativistic quantum mechanics (Selected Papers, Berlin: Springer pp. 115).
  * [36] Margolus N and Levitin L B 1998 Physica D 120, 188.
  * [37] Levitin L B and Toffoli T 2009 Phys. Rev. Lett. 103, 160502.
  * [38] Jozsa R 1994 J. Mod. Opt. 41, 2315.
  * [39] Uhlmann A 1976 Rep. Math. Phys. 9, 273.
  * [40] Luo S and Zhang Q 2004 Phys. Rev. A 69, 032106.
  * [41] Wu S-x and Yu C-s 2020 Sci. Rep. 10, 5500.
  * [42] Deffner S and Lutz E 2013 Phys. Rev. Lett. 111, 010402.
  * [43] Wu S-x and Yu C-s 2018 Phys. Rev. A 98, 042132.
  * [44] Cai X and Zheng Y 2017 Phys. Rev. A 95, 052104.
  * [45] del Campo A, Egusquiza I L, Plenio M B and Huelga S F 2013 Phys. Rev. Lett. 110, 050403.
  * [46] Ektesabi A, Behzadi N and Faizi E 2017 Phys. Rev. A 95, 022115.
  * [47] Liu H B, Yang W L, An J H and Xu Z Y 2016 Phys. Rev. A 93, 020105.
  * [48] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, Oxford, UK: Cambridge University Press, 2007).
  * [49] Davies, E.B.: Quantum theory of open systems (Academic Press, London 1976).
  * [50] Wolf, M M, Eisert J, Cubitt T S and Cirac J I 2008 Phys. Rev. Lett. 101, 150402.
  * [51] Breuer H P, Laine E M and Piilo J 2009 Phys. Rev. Lett. 103, 210401.
  * [52] Rivas A, Huelga S F and Plenio M B 2010 Phys. Rev. Lett. 105, 050403.
  * [53] Luo S, Fu S and Song H 2012 Phys. Rev. A 86, 044101.
  * [54] Zeng H S, Tang N, Zheng Y P and Wang G Y 2011 Phys. Rev. A 84, 032118.
  * [55] He Z, Zeng H S, Li Y, Wang Q and Yao C 2017 Phys. Rev. A 96, 022106.
  * [56] Zhang Y J, Han W, Xia Y J, Cao J P and Fan H 2014 Sci. Rep. 4, 4890.
  * [57] Cimmarusti A D, Yan Z, Patterson B D, Corcos L P, Orozco L A and Deffner S 2015 Phys. Rev. Lett. 114, 233602.
  * [58] Zhang Y J, Han W, Xia Y J, Cao J P and Fan H 2015 Phys. Rev. A 91, 032112.
  * [59] Cianciaruso M, Maniscalco S and Adesso G 2017 Phys. Rev. A 96, 012105.
  * [60] Sun Z, Liu J, Ma J and Wang X 2015 Sci. Rep. 5, 8444.
  * [61] Xu K, Zhang G F and Liu W M 2019 Phys. Rev. A 100, 052305.
  * [62] Xu Z Y, Luo S, YangWL, Liu C and Zhu S 2014 Phys. Rev. A 89, 012307.
  * [63] Teittinen J, Lyyra H and Maniscalco S 2019 New J. Phys. 21, 123041.
  * [64] Berrada K 2018 Physica E 95, 6.
  * [65] Shahri Y, Hadipour M, Haddadi S, Dolatkhah H, and Haseli S 2023 Phys. Lett. A 470, 128783.
  * [66] Gholizadeh A, Hadipour M, Haseli S, Haddadi S, Dolatkhah H 2023 Commun. Theor. Phys. 75, 075101.
  * [67] Hadipour M, Haseli S, Dolatkhah H, Haddadi S, Czerwinski A 2022 Photonics 9, 875.
  * [68] Mazzola L, Maniscalco S, Piilo J, Suominen K A, and Garraway B M 2009 Phys. Rev. A 80, 012104.
  * [69] Xu K, Zhu H J, Zhang G F, Liu W M 2021 Phys. Rev. E 104, 064143.
  * [70] Garraway B M 1997 Phys. Rev. A 55, 2290.
  * [71] Garraway B M 1997 Phys. Rev. A 55, 4636.
  * [72] Mazzola L, Maniscalco S, Piilo J, Suominen K A, and Garraway B M 2009 Phys. Rev. A 80, 012104.
  * [73] Pleasance G, Garraway B M, and Petruccione F 2020 Phys. Rev. Research 2, 043058.


# Interface Retaining Coarsening of Multiphase Flows

Xianyang Chen  Johns Hopkins University, MD, USA  Jiacai Lu and Gretar
Tryggvason  Johns Hopkins University, MD, USA

###### Abstract

Multiphase flows are characterized by sharp moving interfaces, separating
different fluids or phases. In many cases the dynamics of the interface
determines the behavior of the flow. In a coarse, or reduced order model, it
may therefore be important to retain a sharp interface for the resolved
scales. Here, a process to coarsen or filter fully resolved numerical
solutions for incompressible multiphase flows while retaining a sharp
interface is examined. The different phases are identified by an index
function that takes different values in each phase and is coarsened by solving
a constant coefficient diffusion equation, while tracking the interface
contour. Small flows scales of one phase, left behind when the interface is
moved, are embedded in the other phase by solving another diffusion equation
with a modified diffusion coefficient that is zero at the interface location
to prevent diffusion across the interface, plus a pressure like equation to
enforce incompressibility of the coarse velocity field. Examples of different
levels of coarsening are shown. A simulation of a coarse model, where small
scales are treated as a homogeneous mixture, results in a solution that is
similar to the filtered fully resolved field for the early time Rayleigh-
Taylor instability.

##  1 Introduction

Multiphase flows, just like single phase flows, are usually unsteady, with a
large range of spatial and temporal scales. Although the governing equations
are known and it is in principle possible to numerically simulate the time
dependent evolution, for most problems of practical interest the large range
of scales makes doing so impractical. Experience also suggest it is unlikely
to be necessary since in most cases the flows exhibit a fair degree of
universality at the smallest scales. For single phase flows it is by now
fairly well recognized that while it is difficult to develop general closure
laws for mathematical models of the fully stationary average flow, models
where the unsteady motion of the large scales is simulated and models are only
used to describe the average motion of the small scales, can lead to much
improved predictions. Indeed, in many cases it is found that models for the
small scales can be much less elaborate than closure models for the average
motion, since the range of scales that needs to be modeled is greatly reduced.
This approach is often referred to as Large Eddy Simulations (LES). LES
originally referred to simulations where only very small scales were modeled,
but it is now often used indiscriminately for any unsteady simulation
incorporating models for unresolved scales. For multiphase flows, where a
sharp interface separating large regions of different phases, it is likely
that a similar approach, where the unsteady motion of the large scales is
tracked and small scales are modeled, is the most realistic way to develop
accurate predictive strategies.

Considerable effort has been put into developing LES like models for
multiphase flows. Sagaut and Germano [2005] discussed filtering flow fields
near interfaces, noting that the filtering should smooth the turbulence but
preserve the discontinuity, and derived consistency conditions for possible
subgrid models. Filtered LES equations were derived in Labourasse et al.
[2007], who identified the unresolved terms, and showed results of an a priori
test, where proposed closure terms are compared with filtered DNS results for
two two-dimensional flows, one with small vortices interacting with a bubble
and another for “phase inversion” where a light fluid at the bottom of the
computational domain moves to the top. This study was extended to three-
dimensional flows by Toutant et al. [2008] who proposed closure relations for
some of the subgrid terms and examined the interaction of small scale vortices
with a single deformable bubble. Further studies of the fully three-
dimensional phase inversion problem was done by Vincent et al. [2008] who also
derived the averaged LES equations and examined the magnitude of the various
subgrid terms. The jump conditions for interfaces in filtered fields were
revisited by Toutant et al. [2009a] and Toutant et al. [2009b] who in addition
to filtering the turbulent flow simplified the interface, and proposed subgrid
models for some of the unresolved terms. Model equations were also derived by
Liovic and Lakehal [2007] who used “component-weighted volume averaging” to
filter the equations and develop models for unresolved subgrid terms near
deforming interfaces, based on the distance to the interface. Aniszewski et
al. [2012] used approximate deconvolution to model the surface tension tensor
for the filtered equations and Wacławczyk and Oberlack [2011] discussed
modeling in the context of ensemble-averaged fields. A very different approach
to subgrid modeling was introduced by Hermann and Gorokhovski [2008, 2009] and
Herrmann [2013] who used a dual scale strategy where the flow field is
coarsened but the full interface is retained on a finer grid. The filtered
equations are similar to those in Toutant et al. [2009a] but the subgrid terms
for surface tension can be computed directly for the well resolved interface.
A subgrid model is, hoowever, need for the advection of the interface on the
finer grid. Recent development of closure terms for the model equations
derived in Labourasse et al. [2007], focusing on atomization, can be found in
Ketterl and Klein [2018], Klein et al. [2019] and Ketterl et al. [2019].
Studies of the use of deconvolution models for a two-dimensional phase
inversion are described in Gouénard et al. [2019] and for attempts to deal
with non-isothermal flow see Saeedipour et al. [2019]. A large number of
authors have developed numerical models where small bubbles and drops are
represented as Lagrangian point particles while interfaces are tracked using
any of the many numerical method for flows with sharp interfaces (Tomar et al.
[2010], Ma et al. [2017], Hsiao et al. [2017], Zuzio et al. [2018], for
example). Often the point particles are formed by breaking off a small fluid
blob from an interface, usually in an ad-hoc manner. The interface separating
two fluids consists, however, not just of relatively “resolvable” part and
fluid particles, but usually has waves, splashes, crowns, finger and
folds—some of which may lead to bubbles or drops and others that do not—but
that exists at scales too small to be resolved in simulations of real systems.
Another effort to retain large scale structure but model small scales can be
found in Nykteri et al. [2020] who use a hybrid method where a sharp interface
method is used for some part of the computational domain, and a $\Sigma-Y$
two-fluid model, introduced by Navarro-Martinez [2014], is used for the
unresolved parts, with a dynamic switching between the different approaches.
Recent discussions of the status of LES modeling for multiphase flows can be
found in Vincent et al. [2018], Lakehal [2018], Mukundan et al. [2020] and
Nykteri et al. [2020].

As Toutant et al. [2009a] point out, coarse models are significantly different
from single phase LES models because of the presence of the index function and
a few author have sought to make that clear by suggesting different names.
Thus, Toutant et al. [2009a] talk about the Interfaces and Sub-grid Scales
(ISS) approach and Lakehal and Labois [2011] about Large-Eddy & Interface
Simulation (LEIS) models. Later papers seem, however, to have settled on
referring to methods based on filtering as “LES for multiphase flows.” Here we
will refer to the coarse field and talk about coarse flow models.

For single phase flows, coarsening the flow field by filtering to eliminate
high wave number content is a well-established procedure, although accounting
for the effect of the small scales by a closure model is still ongoing
research. The main challenge for multiphase flow is that filters for single
phase flows smooth everything, including the phase boundaries. Since the
interface separating the fluids is often a dominant feature of the flow, it
seems that a better approach would be to coarsen the flow in a way that
retains the interface, at least in many parts of the flow. Interface retaining
coarsening is a common issue in many areas such as cartography where a sharp
shoreline, although simplified, is retained on large-scale maps, and in image
processing where it is found that retaining interfaces and steep gradients
gives “better” coarsened images (Pal et al. [2015]). Apart from simulations
where small droplets are broken off the liquid as Lagrangian drops in
simulations of atomization (Tomar et al. [2010] and other references cited
above) and the work of Nykteri et al. [2020] where small scales are
represented by a two-fluid model, we are only aware of two efforts to retain
sharp interfaces in modeling of multipahse flows. In Toutant et al. [2009a]
the index field identifying the different phases is first smoothed by applying
a Gaussian filter and then the interface is reconstructed by identifying the
contour that originally coincided with the interface. Recognizing that the
interface boundary condition can change, Toutant et al. [2009b] used
asymptotic expansion to derive new boundary conditions for flows with bubbles,
for situations where the filter size is much smaller than the bubbles,
justified by assuming that although the bubble surface is disturbed by the
turbulence, the Kolmogorov length scale is much smaller than the bubble
diameter. The other discussion of interface retaining coarsening is a very
short section in Lakehal [2018] where an “All-Regime Multi-fluid model” is
introduced. Although preliminary results are shown, little details and no
reference are provided and no further development seems to have taken place.
In all other studies of coarse models for multiphase flows, referenced above,
the index function seems to have been filtered along with the rest of the flow
field.

Development of closure models in single phase flows have traditionally relied
on analytical models that are calibrated with results from direct numerical
simulations where the unresolved terms can be computed exactly, since their
form, in terms of the resolved and the filtered variables is known. Most
authors developing models for multiphase flows have taken a similar approach.
Recent developments in a number of areas, which we will collectively refer to
as data-driven modeling, however, suggests a different approach. Several
authors have, in particular, suggested building partial differential equations
directly from data (Raissi and Karniadakis [2018], Long et al. [2018], Lee et
al. [2020]). In this approach the “learned” equations ensure that the coarse
field evolves correctly, and the fully resolved field plays no role except as
the starting point for the coarsening. For fluid flows we have, of course, a
fairly good idea what the overall structure of the equations are, so instead
of having to learn the full equations, we should only have to learn how to
modify them or add extra terms. Machine learning can, of course, also be used
to find closure terms from the fully resolved solution (see Ma et al. [2015,
2016]), but the possibility of obtaining the modifications directly from the
coarse data offers possible new strategies. Indeed, one of the main
consequences of approaching closure modeling from this perspective is that the
coarsening can be done in many ways and we do not need an explicit connections
between the extra terms and the fully resolved field. For a simple
demonstration of this approach for 2D single phase turbulence, see Chen et al.
[2021a].

Here we explore coarsening the phase distribution and the flow field where a
sharp but simplified interface is retained. The index function is coarsened by
solving a constant coefficient diffusion equation, while tracking the
interface contour. Patches of one phase, left behind when the interface is
moved, are embedded in the other phase by solving another diffusion equation
with a modified diffusion coefficient that is zero at the interface location
to prevent diffusion across the interface. When smoothing the the momentum, we
enforce incompressibility of the coarsened flow field by also including a
pressure field. The small scales can be treated in many ways, including by the
models referred to above, but here we only give one example, using a simple
homogeneous mixture model with closure models determined by inspection.

##  2 Scope

We consider unsteady incompressible flow consisting of different fluids or
phases governed by the Navier Stokes equations

$\frac{\partial\rho{\bm{u}}}{\partial
t}+\nabla\cdot(\rho{\bm{u}}{\bm{u}})=-\nabla
p+\rho{\bm{g}}+\nabla\cdot\mu(\nabla{\bm{u}}+\nabla{\bm{u}}^{T})+{\bm{f}}_{\sigma}\quad\hbox{and}\quad\nabla\cdot{\bm{u}}=0.$
(1)

Here, $\rho$ is the density, $\mu$ is the viscosity, ${\bm{u}}$ is the
velocity, $p$ is the pressure, ${\bm{g}}$ is the gravity acceleration and
${\bm{f}}_{\sigma}$ is the surface tension term. Both the density and the
viscosity are different in the different fluids. Solving these equations
accurately gives the fully resolved flow field at any given time and spatial
location. The different phases are identified by an index or marker function
$\chi$ and assuming, for simplicity, that only two phases are involved, we
have:

$\chi({\bm{x}})=\left\\{\begin{array}[]{l l}\text{0 in fluid 0}\\\ \text{1 in
fluid 1}.\end{array}\right.$ (2)

The various flow quantities, such as density, are then given by
$\rho=\chi\rho_{1}+(1-\chi)\rho_{0}$.

The full solution generally contains a large range of scales and the purpose
of a coarser model is to remove the smallest ones, yet account for their
effects. For a single phase flow where filtering to remove the high wavenumber
components is generally used, the coarse flow is simply a smoother version of
the full solution. For multiphase flows, where a sharp interface separates the
phases, one possibility is to smooth the index function along with the
velocity field, thus removing the sharp interface. Here, however, we pursue a
different strategy and retain a sharp interface, although its shape will
generally be simplified. Thus, we seek a coarsening strategy that has the
following characteristics:

  * •

The large and the small scale phase distribution are separated by coarsening
the index function. The interface is simplified but stays sharp, with some
parts “collapsing” into points, lines or sheets of zero volume.

  * •

The large and the small scale flow on either side of the interface are
separated by coarsening. In regions where the smoothing of the index function
causes fluid to “switch sides,” the fluid is mixed with or embedded in the
other fluid.

  * •

As the coarsening is reduced, the flow field approaches the results given by
direct numerical simulations (DNS).

Ideally, the coarsening strategy should be general in that we can coarsen just
a little bit as well as very aggressive. Figure 1 shows the process
schematically. Here we assume a very aggressive filtering so that the coarse
index function is much simpler than for the fully resolved flow. As we filter
the index function, we also filter the flow variables, creating mixed zoned on
either side of the interface, as shown on the right. Notice that while the
densities and other material properties in the original fully resolved flow
may be constant, that is not generally true for the coarsened flow.

The coarsening is a precursor to the development or implementation of
multiscale models to evolve the coarse flow field. The flow in the mixed zones
can be modeled in a variety of ways, ranging from simple mixture models, drift
flux model, full two-fluid models or by using Lagrangian point particles to
account for bubbles in liquids and drops in gas. Away from interfaces, models
for the flow field should become standard models for multiphase and/or
turbulent flows but close to the interface models need to developed for the
additional terms, using analysis and/or machine learning.

Figure 1: Aggressive filtering of an idealized starting jet. The original
index function is shown in the left frame and the filtered index function is
shown on the right. The large scale flow is represented by the gray area and
the mixed zones by the blue regions. The original interface is shown by a thin
solid line.

##  3 Coarsening the Fully Resolved Field

We start by the observation that the solution of the linear diffusion equation

$\frac{\partial g(t,{\bf x})}{\partial t}=D\nabla^{2}g(t,{\bf x}),$ (3)

in an unbounded two-dimensional domain, is given by

$g({\bf x},T)=\frac{1}{4\pi DT}\int_{Area}e^{-\frac{||{\bf x}-{\bf
x}^{\prime}||^{2}}{4DT}}g_{o}({\bf x}^{\prime})da^{\prime},$ (4)

at time $T$, where $g_{o}$ is the initial condition. Thus, filtering a field
by a Gaussian kernel is equivalent to evolving it by a diffusion equation for
a time that is related to the length scale of the filter $\Delta$ by
$4DT=\Delta^{2}/6$. The reason we find it more convenient to coarsen the field
by diffusion rather than filtering, should become clear in the rest of this
section. The connection between filtering and diffusion for single phase
turbulence has been pointed out by Johnson [2020, 2021], and Capecelatro and
Desjardins [2013] use diffusion to distribute the effect of small particles
onto a fixed grid, as examples of other uses of diffusion instead of explicit
filtering. We have taken advantage of the flexibility of doing smoothing by
diffusion rather than filtering in Chen et al. [2021b], where we used it to
sharpen a smooth distribution.

Figure 2: The motion of an interface following the interface contour of the
index function. Left: the interface, where the index function has the value
$\chi^{*}$, needs to catch up to $\chi_{I}$ and must move a distance $\Delta
x=u_{n}\Delta\tau$, computed from $\chi_{I}$ and the gradient. Right: Once
interfaces are close enough, there is no $\chi_{I}$ contour close by.

###  3.1 Smoothing the Index Function

To simplify the interface we diffuse the index function by solving

$\frac{\partial\chi}{\partial\tau}=\nabla^{2}\chi,$ (5)

in “pseudo-time” $\tau$, until $\tau=T$, at every instance in “real time” when
we want a smoothed solution. Notice that we set the diffusion coefficient to
unity ($D=1$) since a different value simply rescales $\tau$. To move the
interface with the contour identifying the interface, $\chi=\chi_{I}$, we find
the “diffusion velocity,” using that the value of the function at the
interface does not change, so the material derivative of $\chi$ is zero:

$\frac{D\chi}{D\tau}=\frac{\partial\chi}{\partial\tau}+{\bm{u}}\cdot\nabla\chi=0.$
(6)

We only need the velocity in the normal direction to the interface, and a
normal is defined by ${\bm{n}}=\nabla\chi/|\nabla\chi|$, so we write
${\bm{u}}=u_{n}{\bm{n}}=u_{n}(\nabla\chi/|\nabla\chi|)$. Rearranging and
dividing by $|\nabla\chi|$ gives

$\frac{1}{|\nabla\chi|}\frac{\partial\chi}{\partial\tau}=-u_{n}\frac{\nabla\chi\cdot\nabla\chi}{|\nabla\chi|^{2}}=-u_{n},$
(7)

since ${\bm{n}}\cdot{\bm{n}}=1$. After the index field has been updated by
taking one step in pseudo-time, the interface usually no longer coincides with
the interface contour $\chi_{I}$ and for an interface point that is not
exactly on the interface contour, the value differs by
$\Delta\chi=\chi_{I}-\chi^{*}$ from the interface value, where $\chi^{*}$ is
the old value of the index function at the interface. See figure 2. In a time
step $\Delta\tau$ the interface thus needs to “catch up,” so we approximate
$\partial\chi/\partial t\approx\Delta\chi/\Delta\tau$ and write

${\bm{u}}_{I}=u_{n}{\bm{n}}=-\frac{(\chi_{I}-\chi^{*})}{|\nabla\chi|^{2}\Delta\tau}\nabla\chi.$
(8)

We interpolate $\chi$ from the grid to find $\Delta\chi$ and $\nabla\chi$,
from which we can find $|\nabla\chi|^{2}=\nabla\chi\cdot\nabla\chi$. Each
interface point is then moved by solving

$\frac{d{\bm{x}}_{I}}{d\tau}={\bm{u}}_{I}.$ (9)

Figure 3: The evolution of an interface simplified by diffusing the index
function at times $t=0.175$ (top row) and $t=0.25$ (bottom row). In each row
the results are shown for $\tau=0$ (original), $\tau=0.0004$ ($\Delta=0.1$),
$\tau=0.0026$ ($\Delta=0.25$), and $\tau=0.0104$ ($\Delta=0.5$). The interface
is the red line and contours of the index function are also shown. “Collapsed”
interfaces are shown as gray lines.

An interface moves when the index function is smoothed by diffusion for two
reasons: Curvature of the interface and the presence of a close, parallel,
interface. Solving a diffusion equation accounts for motions due to both
reasons, but whereas an interface will simply stop moving when it becomes
straight, an interface moving toward another interface is likely to eventually
cross since the extreme value of the index function will be smaller or larger
than the interface value and there is therefore no nearby point where the
value of $\chi$ is equal to the interface value. In those cases it is
incorrect to approximate the time derivative as we did above and we simply
want the interface to stop moving. To check for this, we compute a normal to
the interface and examine the value of $\chi$ at a point slightly ahead of the
interface. In most cases this value will be larger than the current value of
$\chi$, if $\chi$ is smaller than the interface value, and smaller if $\chi$
is larger than the interface value, but if both values are on the same side of
the interface value, then it is likely that another interface is nearby and we
put the interface velocity to zero. We usually also flag such “collapsed”
interfaces since often we may want to ignore them when modeling the mixed
zone. We note that our procedure for smoothing the index function is
essentially identical to the one used by Toutant et al. [2009a], except that
we follow the interface as the field is smoothed, instead of first smoothing
and then restoring the interface afterwards from the contour line identified
with the interface.

Figure 4: A plot of interface length versus pseudo-diffusion-time. The solid
line is the large scale interface length while the dashed line is the
collapsed interface length.

There are, of course, other possibilities to move the interface as the index
function is smoothed. Instead of approximating the time derivative as done
above, we can use equation (5) and replace the time derivative by the
Laplacian. Computing the Laplacian on the fixed grid and interpolating it to
the interface works and should, in principle, handle interacting interfaces
correctly. We have, however, not found it as robust as the process described
above. Since the velocity of an interface in a diffusion field is given by its
curvature, we could also move it that way. It is, however, first of all,
generally difficult to accurately evaluate the curvature and secondly, an
interface evolved by its curvature is not guarantied to follow the appropriate
contour as is the case with our approach, and numerical errors can therefore
accumulate.

To show how the coarsening works, we apply it to the large amplitude stage of
a two-dimensional Rayleigh-Taylor instability. The computational domain is a
rectangle of dimensions $1\times 2$ with periodic side boundaries and rigid
top and bottom, and resolved by a regular $128\times 256$ grid. The density
and viscosity of the heavy and light fluids are $\rho_{1}=10$, $\mu_{1}=0.02$,
and $\rho_{0}=5$ and $\mu_{0}=0.01$, respectively. Gravity acceleration is
$g_{y}=-100$ and surface tension is $\sigma=0.25$. At the initial time the
velocities are zero but the elevation of the interface is perturbed by
$y(x)=1+0.1\sin(2\pi x)+0.1\sin(4\pi x)$. Figure 3 shows the simplification of
the interface at times $t=0.175$ and $t=0.25$. The filtering is shown as an
evolution in pseudo time $\tau$, starting with the unfiltered interface on the
left and progressing to the right. Stopping at a given pseudo-time correspond
to a given filter size. The red line is the tracked interface and we also show
contours of the diffused index function. The gray traces are inactive (or
“collapsed”) interfaces left over as the simplification progresses. In the
second frame, at $\tau=0.0004$, only the very smallest scales have been
eliminated and the large scale structure of the interface is intact, except
that at the later time the upward moving protrusion seen at the earlier time
has just detached from the lower fluid. In the third frame significantly more
simplification has taken place and the downward protrusion moving downward at
the earlier time has separated from the rest of the heavy fluid at the later
time, so at the later time the flows consists of a relatively flat interface
and light fluid blob moving up and heavy blob moving down. In the forth frame,
at $\tau=0.0104$, the interface has been smoothed significantly at the earlier
time and at the later time it consists of a nearly flat interface and one blob
of the light fluid moving upward and a downward moving mixed region. If we
continued to smooth the interface, eventually we end up with a flat interface
and mixed regions moving up and down.

The interface length is shown versus pseudo-time in figure 4, for both times
shown in figure 3. The interface length decreases as the interface is
simplified, as expected. In addition to the length of the interface separating
regions of different phases (solid lines) we also plot the length of
interfaces that have “collapsed” and no longer serve as a phase boundaries
(dashed lines). The length of those increases with time and the figure shows
that the collapse takes place rapidly, leading to kinks in the curves. The
interface length is reduced in two ways: by straightening the interface since
diffusion corresponds to interface motion by mean curvature and by interfaces
disappearing as close by interfaces collapse onto each other. For this
problem, the plot shows that most of the shortening of the interface is due to
interfaces collapsing.

Figure 5: The evolution of the interface in pseudo-time as the interface is
smoothed by diffusion for a three-dimensional Rayleigh-Taylor instability. The
original interface is shown in the frame on the left and the contours for the
index function in the gray plane are shown in the four frames on the right for
$\tau=0$ (original), $\tau=0.0026$ ($\Delta=0.25$), $\tau=0.0204$
($\Delta=0.7$), and $\tau=0.0267$ ($\Delta=0.8$). The red line is the
intersection of the interface with the plane.

Although we have chosen to show the simplification for two-dimensional flows,
the process described above applies to fully three-dimensional flows as well
and in figure 5 we show the simplification for a Rayleigh-Taylor instability
at one “real” time, versus pseudo-time. In the frame on the left the three-
dimensional interface is shown along with a plane cutting diagonally through
the computational domain. The interface and contours of the index functions in
the gray plane are shown for three pseudo-times in the frames on the right.
The computational domain is a periodic hexahedron of dimensions $1\times
1\times 3$ with periodic side boundaries and rigid top and bottom, resolved by
a regular $64\times 64\times 192$ grid. The density and viscosity of the heavy
and light fluids are $\rho_{1}=10$, $\mu_{1}=0.08$, and $\rho_{0}=5$ and
$\mu_{0}=0.04$, respectively. Gravity acceleration is $g_{y}=-100$ and surface
tension is zero. At the initial time the velocities are zero but the elevation
of the interface is perturbed by $z(x,y)=1.75+0.12\times[\cos(2\pi
x)+\cos(2\pi y)]$. The small scales have been eliminated at $\tau=0.0026$, at
$\tau=0.0204$ a downward moving drop has separated from the heavy fluid, and
in the last frame ( $\tau=0.0267$) this drop has been mixed with the lighter
fluid, as seen for the two-dimensional case.

###  3.2 Smoothing the Flow Variables

As the interface moves, field variables such as density and momentum are left
behind on the opposite side of the interface. Variables that switch side are
taken to belong to the small scales and we can account for their evolution in
many ways. Here we smooth those and “blend” them with the original field.
Thus, we need to diffuse those variables, but only on one side of the
interface. To do so we again solve a diffusion equation, but to prevent
diffusion across the interface, we put the diffusion coefficients around the
interface equal to zero in a thin band containing the interface. We denote the
new diffusion coefficient by $\tilde{D}$ and since the diffusion coefficient
is now not constant, we need to solve

$\frac{\partial\phi}{\partial\tau}=\nabla\cdot\tilde{D}\nabla\phi,$ (10)

where $\phi$ is the flow variable that we are smoothing. Smoothing the index
function again, but using $\tilde{D}$ instead of a constant $D$, gives the
volume fraction $\alpha$ of the phase where $\chi=1$ for the unfiltered field.

Unlike for single phase flow, where we generally filter the velocity field,
here we filter the momentum, since it is the conserved quantity. This can lead
to nonzero divergence at the interface and to correct for that we recompute a
pressure field needed to make the velocity divergence free. For the mixture
model used here, where both phases have the same velocity, the mixture
velocity should be divergence free. For a more advanced model, such as the
drift flux model, we would have slip and the velocity used here becomes the
mass weighted velocity which is not necessarily divergence free. However, once
the slip is known, the volume source can be computed.

The complete smoothing process involves evolving the index function and the
flow variables in pseudo-time $\tau$ for long enough to achieve the desired
smoothing. In addition to diffusing the index function to move the interface
(equations 5 and 8), we also evolve a copy of the index function to find the
void fraction by using $\tilde{D}$, which consists of the constant diffusion
coefficient $D$ modified by putting it to zero at and around the interface.
Similarly, we evolve the momentum using a pseudo-pressure field to enforce
incompressibility.

Denoting the filtered variables by a tilde, the equations to be evolved in
pseudo time ($\tau\leq T$) to generate a coarse flow field are:

$\displaystyle\tilde{\chi}(\tau=0)=\chi;\quad\alpha(\tau=0)=\chi;\quad D=1.0;$
(11)
$\displaystyle\frac{\partial\tilde{\chi}}{\partial\tau}=\nabla^{2}\tilde{\chi};\qquad{\bm{u}}_{f}=-\frac{1}{\Delta\tau}\Biggl{(}\frac{(\tilde{\chi}_{I}-\tilde{\chi})}{|\nabla\tilde{\chi}|^{2}}\Biggr{)}\nabla\tilde{\chi};\qquad\frac{d{\bm{x}}_{f}}{d\tau}={\bm{u}}_{f};$
(12) $\displaystyle\tilde{D}({\bm{x}},\tau)=D\hbox{ modified by setting
}D({\bm{x}}_{f})=0;$ (13)
$\displaystyle\frac{\partial\alpha}{\partial\tau}=\nabla\cdot\tilde{D}\nabla\alpha;$
(14)
$\displaystyle\frac{\partial}{\partial\tau}(\widetilde{\rho{\bm{u}}})=-\nabla\tilde{p}+\nabla\cdot\tilde{D}\nabla(\widetilde{\rho{\bm{u}}});\quad\hbox{with}\quad\nabla\cdot\tilde{\bm{u}}=0.$
(15)

We assume that $\widetilde{\rho{\bm{u}}}=\tilde{\rho}\tilde{\bm{u}}$, that the
density is given by $\tilde{\rho}=\alpha\rho_{1}+(1-\alpha)\rho_{o}$, and that
the filtered velocity is incompressible. The second equation (12) is the
smoothing of the index function described above. The third equation (13) is
the modification of the diffusion coefficient by setting it to zero where the
interface is. In equation (14) we diffuse a copy of the index function to find
the volume fraction, and the last equation (15) is the evolution of the
momentum using a projection method to enforce incompressibility by finding the
appropriate pressure. Enforcing incompressibility of the coarse field would be
difficult if the flow field is smoothed directly using a filter (Toutant et
al. [2009b]) but seems both important and reasonable to do. If we take the
slip velocity to be zero, as in the simple mixture model used in the next
section, there is only one velocity and the divergence is zero. Once the
interface has been moved by the smoothing, the coarse filtered index function
$\tilde{\chi}$ identifying the large scale phase distribution is reconstructed
by putting it equal to $1$ in one fluid and $0$ in the other. When the
coarsening leaves most of the phases separated by a sharp interface, we assume
that relatively small amount of phase 0 is mixed with phase 1 and vice versa.
For the coarse field the index function $\tilde{\chi}$ identifies the
different fluids but $\alpha$ is equal to the diffused index function $\chi$.
If there is no mixing, $\alpha=\tilde{\chi}=0$ in fluid $0$ and
$\alpha=\tilde{\chi}=1$ in fluid $1$. We define a perturbation void fraction
$\alpha^{\prime}=\alpha-\tilde{\chi}$ which takes negative values in fluid 1
and a positive values in fluid 0 and quantifies how much of one fluid is mixed
in the other fluid. In general we expect $\alpha^{\prime}$ to consists of
isolated regions of positive and negative values and to be zero for most of
the flow field.

Figure 6 shows the original sharp interface (left frame) and the perturbation
volume fraction at three stages of smoothing, for the second time shown in
figure 3 and at the same pseudo times. Where the interface is moved by the
smoothing, one phase is “left behind” and is mixed with the phase originally
on that side of the interface. Initially most of the mixing takes place near
high curvature regions of the interface, where it is retreating rapidly, but
as small scale features are eliminated, we see regions of high mixture
fractions inside each domain. We note that for interfaces that have collapsed
into thin sheets, we do not put the diffusion to zero, although we have the
option of doing so.

In figure 7 the streamfunction computed from the velocity field is shown,
first for the unfiltered velocity (left frame) and then for the smoothed
velocity (frames 2-4) at the same pseudo times as in figure 3 and 6. Since we
enforce incompressibility when we smooth the flow, the normal velocity remains
continuous across the simplified interface. For modest smoothing, where only
the smallest scales have been eliminated the flow field remains close to the
original one (second and third frame) but for aggressive smoothing as in the
last frame the flow field above and below the interface has been simplifies to
a pair of counter rotating vortices, driven by the upward moving blob and the
downward moving mixed region.

Figure 6: The evolution of the perturbation volume fraction field
($\alpha^{\prime}$) in pseudo-time by nonlinear diffusion, for the later time
shown in figure 3 ($t=0.25$) and the same values of $\tau$.

Figure 7: The evolution of the stream function in pseudo-time as the interface
is smoothed by diffusion, for the later time shown in figure 3 ($t=0.25$) and
the same values of $\tau$.

Figure 8: The evolution of the vorticity in pseudo-time as the interface is
smoothed by diffusion, for the later time shown in figure 3 ($t=0.25$) and the
same values of $\tau$.

Figure 8 shows the vorticity field computed from the velocity field, first for
the original velocity (left frame) and then for the smoothed velocity (frames
2-4) at the same pseudo times as in figures 3 – 7. Since the flow is driven by
baroclinic vorticity generated at the interface, the concentration is highest
there and indeed, nearly all the original un-smoothed vorticity is at the
interface. As small scale features are eliminated, the vorticity is mixed with
the bulk, although the concentration remains highest at the interface. For the
most aggressive smoothing in the right frame, we see pairs of faint counter-
rotating vortices in the bulk on either side of the interface, corresponding
to the vortices seen in the plot of the streamfunction in figure 7.

We can generate field values for other quantities than those show here. The
interfacial area, for example, is one important measure of the simplification
of the phase distribution and we can construct an area field by
differentiating the unfiltered index function. By filtering it we get the area
concentration. Similarly, we can filter components of the area tensor. The
volume fraction and the area concentration would allow us to obtain the Sauter
Mean diameter and equivalent number density, if we assume that the smoothed
phase consists of spherical drops or bubbles.

##  4 Evolving the Coarse Field

The main purpose of coarsening the flow field is to provide data for modeling
the evolution of the coarse flow. Such models come in various forms and the
purpose of the present paper is not to explore the many possibilities in any
detail. We do, however, include one simple example in this section.

The homogeneous mixture model, where the phases are assumed to be completely
mixed and move with the large scale fluid velocity is probably the simplest
model imaginable. In our implementation a sharp interface is assumed to
separate the large scales, and the small scales are represented by mixtures
embedded in the large scale regions, so the conversion between large resolved
scales and modeled small scales must be captured. The coarse interface moves
with the coarse velocity when no conversion happens, but when transfer of the
index function between large and small scales takes place, the interface and
fluid velocity are different. We generally expect this to happen in high
curvature regions. There is only one large scale velocity which is found by
solving a momentum equation, supplemented by the incompressibility conditions:

$\frac{\partial({\tilde{\rho}\tilde{\bm{u}}})}{\partial
t}+\nabla\cdot(\tilde{\rho}\tilde{\bm{u}}\tilde{\bm{u}})=-\nabla\tilde{p}+\tilde{\rho}{\bm{g}}+\nabla\cdot{\bm{\tau}}_{e}\quad\hbox{
and }\quad\nabla\cdot\tilde{\bm{u}}=0.$ (16)

Here we have lump into ${\bm{\tau}}_{e}$ all transport terms due to viscous
stresses, surface tension and small scale velocity fluctuations. The density
field everywhere is given by

$\tilde{\rho}=(\alpha^{\prime}+\tilde{\chi})\rho_{1}+(1-\alpha^{\prime}-\tilde{\chi})\rho_{o},$
(17)

where $\alpha^{\prime}$ is the perturbation volume fraction introduced earlier
and is evolved by an advection-diffusion equation:

$\frac{\partial\alpha^{\prime}}{\partial
t}+{\bm{u}}\cdot\nabla\alpha^{\prime}=\nabla\cdot\tilde{D}_{e}\nabla\alpha^{\prime}+S_{I}.$
(18)

$\tilde{D}_{e}$ is an effective diffusion term, adjusted to match the
spreading of the mixed region in the coarsened field and modified to prevent
diffusion across the interface by setting $\tilde{D}_{e}({\bm{x}}_{I})=0$.
$\tilde{D}_{e}$ can depend on $\alpha$, $\tilde{\chi}$, and other flow
variables, but here we will take it to be a constant. The interface source
term is related to the relative speed of the interface by

$S_{I}=[\chi]({\bm{u}}_{f}-{\bm{u}}_{I})\cdot{\bm{n}}=[\chi]\Delta u_{I}$ (19)

where ${\bm{u}}_{f}$ is the fluid velocity and ${\bm{u}}_{I}$ is the interface
velocity, and $\Delta u_{I}=({\bm{u}}_{f}-{\bm{u}}_{I})\cdot{\bm{n}}$. The
interface moves by ${d{\bm{x}}_{I}/dt}={\bm{u}}_{f}-\Delta u_{I}{\bm{n}}$ so
$\Delta u_{I}$ determines the “production” of small scales where the interface
does not exactly follow the fluid velocity.

To close the equations for the coarse field, we must develop relationships
that describe how $\Delta u_{I}$, $D_{e}$ and ${\bm{\tau}}_{e}$ depend on the
coarse variables. In general, we expect that those would be found by comparing
the evolution of the solution for the coarser field with the evolution of the
coarsened fully resolved solution and correlated by machine learning, for
example. Here, however, we have simply guess a value for $\tilde{D}_{e}$, and
taken the stress term in the momentum equation to be the same as for the
unfiltered flow, including the surface tension. Finding $\Delta u_{I}$
requires a more complex approach. We expect most of the small scale production
to take place in high curvature regions and since motion by mean curvature is
equivalent to diffusion, as noted earlier, the simplest approach is to use the
same approach as we used for the original coarsening and solve a diffusion
equation in pseudo time, after the solution has been advanced assuming the
interface moves with the fluid velocity. Once the interface has been moved,
the phase “left behind” becomes inside the other phase becomes small scale
flow.

Figure 9: The Rayleigh-Taylor instability at time $t=0.2$. The left frame
shows the unfiltered interface, the middle frame shows the interface and the
perturbation volume fraction obtained by smoothing the flow by diffusion. The
solution in the right frame is obtained by the simple mixture model described
in the text.

Figure 9 shows a large amplitude stage of two fluids undergoing a Rayleigh
Taylor instability at time $t=0.2$. The first frame shows the original
unfiltered interface. the middle frame shows the interface and the
perturbation volume fraction obtained by smoothing the flow by diffusion using
$\Delta=0.18$, and the right frame is the solution evolved from the same
initial conditions using the simple homogeneous mixture model with
$\tilde{D}_{e}=0.05$ and $\Delta u_{I}$ found by diffusion in pseudo time,
every twenty time steps, with the modification that we leave the interface
untouched if it moves less than $0.15\times h$, where $h$ is the grid spacing.
Here, $\Delta t=6.1\times 10^{-5}$ and $\Delta\tau=6.1\times 10^{-6}$ and we
take ten steps in pseudo time. The model equations are solved by the same
method as used for the fully original flow and on the same grid. The frequency
of the interface modification was selected such that the interface matches
roughly with the filtered interface, but no effort was made to match the
effective diffusivity closely. At earlier times the agreement is also good,
but after the time shown here the filtered interface starts to undergo
topology changes and a more sophisticated selection of the model parameters
would be necessary.

##  5 Conclusion

We discuss a strategy to coarsen multiphase flows in a consistent way while
retaining a sharp, but simplified interface. The motivation is the possibility
of developing closure modeling through machine learning by working directly
with the coarsened field where an explicit form relating the closure terms to
the fully resolved flow is not needed. Regions without any mixing can be
treated using standard turbulence models and mixed zones, where one phase is
embedded in the other, can be treated in a variety of ways, such as assuming a
homogeneous mixture (as done here), by a drift flux model, an Eularian two
fluid model, or representing the embedded phase as Lagrangian particles.

We note that here we represent the interface by connected marker points that
are moved as the smoothing progresses. To simplify the interface we could also
simply smooth the index function and then reconstruct the interface from the
contour of the average value. For the other steps, where we separate the large
and small flow scales by diffusion and we need to evolve the interface along
with the field values, having the interface makes this relatively
straightforward. This could, however, presumably also be done in the absence
of explicit tracking of the interface by reconstructing the interface from the
average contour at every step in pseudo time.

Although the motivation for the coarsening strategy presented here is reduced
order modeling of the flow, it also provides us with a tool that should be
useful in analyzing the distribution of scales in multiphase flows. As we
coarsen the flow, we separate it into small and large scales, which can be
examined separately as the coarsening is varied. Other possibilities include
following a different contour than the average one to “skeletonize” the flow.
We hope to explore some of those possibilities in later studies.

Acknowledgement: This research was supported in part by NSF grant
CBET-1953082.

## References

  * Aniszewski et al. [2012] W. Aniszewski, A. Bogusławski, M. Marek, and A. Tyliszczak.  A new approach to sub-grid surface tension for LES of two-phase flows.  _Journal of Computational Physics_ , 231:7368–7397, 2012\.
  * Capecelatro and Desjardins [2013] J. Capecelatro and O. Desjardins.  An Euler-Lagrange strategy for simulating particle-laden flows.  _Journal of Computational Physics_ , 238:1–31, 2013.
  * Chen et al. [2021a] X. Chen, J. Lu, and G. Tryggvason.  Finding closure terms directly from coarse data for 2D turbulent flow.  _Submitted for publication_ , 2021a.
  * Chen et al. [2021b] X. Chen, J. Lu, and G. Tryggvason.  Condensing smooth layers into singular sheets by weighted coordinate smoothing.  _Journal of Computational Physics_ , 431:110140, 2021b.
  * Gouénard et al. [2019] S. Gouénard, S. Vincent, and S. Mimouni.  LES of turbulent non-isothermal two-phase flows within a multifield approach.  _Thermodynamics of Interfaces and Fluid Mechanics_ , 3:29 pages, 2019.
  * Hermann and Gorokhovski [2008] M. Hermann and M. Gorokhovski.  An outline of a LES subgrid model for liquid/gas phase interface dynamics.  _Proceedings of the 2008 CTR summer program, center for turbulence research_ , pages 171–181, 2008.
  * Hermann and Gorokhovski [2009] M. Hermann and M. Gorokhovski.  A large eddy simulation subgrid model for turbulent phase interface dynamics.  _11th triennial international annual conference on liquid atomization and spray systems_ , 2009.
  * Herrmann [2013] M. Herrmann.  A sub-grid surface dynamics model for sub-filter surface tension induced interface dynamics.  _Computers & Fluids_, 87:92–101, 2013.
  * Hsiao et al. [2017] C.-T. Hsiao, J. Ma, and G. L. Chahine.  Multiscale two-phase flow modeling of sheet and cloud cavitation.  _International Journal of Multiphase Flow_ , 90:102–117, 2017.
  * Johnson [2020] P. L. Johnson.  Energy transfer from large to small scales in turbulence by multiscale nonlinear strain and vorticity interactions.  _Physical Review Letters_ , 124:104501, 2020.
  * Johnson [2021] P. L. Johnson.  On the role of vorticity stretching and strain self-amplification in the turbulence energy cascade.  _Journal of Fluid Mechanics_ , In Press, 2021.
  * Ketterl and Klein [2018] S. Ketterl and M. Klein.  A-priori assessment of subgrid scale models for large-eddy simulation of multiphase primary breakup.  _Computers and Fluids_ , 165:64–77, 2018.
  * Ketterl et al. [2019] S. Ketterl, M. Reißmann, and M. Klein.  Large eddy simulation of multiphase flows using the volume of fluid method: Part 2-a-posteriori analysis of liquid jet atomization.  _Experimental and Computational Multiphase Flow_ , 1:201–211, 2019.
  * Klein et al. [2019] M. Klein, S. Ketterl, and J. Hasslberger.  Large eddy simulation of multiphase flows using the volume of fluid method: Part 1-governing equations and a priori analysis.  _Experimental and Computational Multiphase Flow_ , 1:130–144, 2019.
  * Labourasse et al. [2007] E. Labourasse, D. Lacanette, A. Toutant, P. Lubin, S. Vincent, O. Lebaigue, J.-P. Caltagirone, and P. Sagaut.  Towards large eddy simulation of isothermal two-phase flows: Governing equations and a priori tests.  _International Journal of Multiphase Flow_ , 33:1–39, 2007\.
  * Lakehal [2018] D. Lakehal.  Status and future developments of large-eddy simulation of turbulent multi-fluid flows (LEIS and LESS).  _International Journal of Multiphase Flow_ , 104:322–337, 2018.
  * Lakehal and Labois [2011] D. Lakehal and M. Labois.  A new modelling strategy for phase-change heat transfer in turbulent interfacial two-phase flow.  _International Journal of Multiphase Flow_ , 37:627–639, 2011.
  * Lee et al. [2020] S. Lee, M. Kooshkbaghi, K. Spiliotis, C. I. Siettos, and I. G. Kevrekidis.  Coarse-scale PDEs from fine-scale observations via machine learning.  _Chaos_ , 30:013141, 2020.
  * Liovic and Lakehal [2007] P. Liovic and D. Lakehal.  Multi-physics treatment in the vicinity of arbitrarily deformable gas–liquid interfaces.  _Journal of Computational Physics_ , 222:504–535, 2007\.
  * Long et al. [2018] Z. Long, Y. Lu, X. Ma, and B. Dong.  Pde-net: Learning pdes from data.  In _Proceedings of the 35th International Conference on Machine Learning, Stockholm, Sweden_ , PMLR 80, 2018.
  * Ma et al. [2017] J. Ma, C.-T. Hsiao, and G. L. Chahine.  A physics based multiscale modeling of cavitating flows.  _Computer & Fluid_, 145:68–84, 2017.
  * Ma et al. [2015] M. Ma, J. Lu, and G. Tryggvason.  Using statistical learning to close two-fluid multiphase flow equations for a simple bubbly system.  _Physics of Fluids_ , 27:092101, 2015.  doi: https://doi.org/10.1063/1.4930004.
  * Ma et al. [2016] M. Ma, J. Lu, and G. Tryggvason.  Using statistical learning to close two-fluid multiphase flow equations for bubbly flows in vertical channels.  _International Journal of Multiphase Flow_ , 85:336–347, 2016.  doi: https://doi.org/10.1016/j.ijmultiphaseflow.2016.06.021.
  * Mukundan et al. [2020] A. A. Mukundan, G. Tretola, T. Ménard, M. Herrmann, S. Navarro-Martinez, K. Vogiatzaki, J. C. B. de Motta, and A. Berlemont.  DNS and LES of primary atomization of turbulent liquid jet injection into a gaseous crossflow environment.  _Proceedings of the Combustion Institute_ , page In Press, 2020.
  * Navarro-Martinez [2014] S. Navarro-Martinez.  Large eddy simulation of spray atomization with a probability density function method.  _International Journal of Multiphase Flow_ , 63:11–22, 2014\.
  * Nykteri et al. [2020] G. Nykteri, P. Koukouvinis, S. R. G. Avila, C.-D. Ohl, and M. Gavaises.  A $\Sigma-\Upsilon$ two-fluid model with dynamic local topology detection: Application to high-speed droplet impact.  _Journal of Computational Physics_ , 408:109225, 2020.
  * Pal et al. [2015] C. Pal, A. Chakrabarti, and R. Ghosh.  A brief survey of recent edge-preserving smoothing algorithms on digital images.  _Procedia Computer Science_ , 00:1–40, 2015.
  * Raissi and Karniadakis [2018] M. Raissi and G. E. Karniadakis.  Hidden physics models: Machine learning of nonlinear partial differential equations.  _Journal of Computational Physics_ , 357:125–141, 2018\.
  * Saeedipour et al. [2019] M. Saeedipour, S. Vincent, and S. Pirker.  Large eddy simulation of turbulent interfacial flows using approximate deconvolution model.  _International Journal of Multiphase Flow_ , 112:286–299, 2019.
  * Sagaut and Germano [2005] P. Sagaut and M. Germano.  On the filtering paradigm for LES of flows with discontinuities.  _Journal of Turbulence_ , 6:N23, 2005.
  * Tomar et al. [2010] G. Tomar, D. Fuster, S. Zaleski, and S. Popinet.  Multiscale simulations of primary atomization.  _Computers & Fluids_, 39:1864–1874, 2010.
  * Toutant et al. [2008] A. Toutant, E. Labourasse, O. Lebaigue, and O. Simonin.  DNS of the interaction between a deformable buoyant bubble and spatially decaying turbulence: a priori tests for LES two-phase flow modelling.  _Computers and Fluids_ , 37:877–886, 2008.
  * Toutant et al. [2009a] A. Toutant, M. Chandesris, D. Jamet, and O. Lebaigue.  Jump conditions for filtered quantities at an under-resolved interface. part 1: Theoretical developments.  _International Journal of Multiphase Flow_ , 35(1100-1118), 2009a.
  * Toutant et al. [2009b] A. Toutant, M. Chandesris, D. Jamet, and O. Lebaigue.  Jump conditions for filtered quantities at an under-resolved interface. part 2: A priori tests.  _International Journal of Multiphase Flow_ , 35:1119–1129, 2009b.
  * Vincent et al. [2008] S. Vincent, J. Larocque, D. Lacanette, A. Toutant, P. Lubin, and P. Sagaut.  Numerical simulations of phase separation and a priori two-phase LES filtering.  _Computers and Fluids_ , 37:898–906, 2008.
  * Vincent et al. [2018] S. Vincent, M. Tavares, S. Fleau, S. Mimouni, M. Ould-Rouiss, and J.-L. Estivalezes.  A priori filtering and LES modeling of turbulent two-phase flows application to phase separation.  _Computers and Fluids_ , 176:245–259, 2018.
  * Wacławczyk and Oberlack [2011] M. Wacławczyk and M. Oberlack.  Closure proposals for the tracking of turbulence-agitated gas–liquid interfaces in stratified flows international.  _Journal of Multiphase Flow_ , 37:967–976, 2011.
  * Zuzio et al. [2018] D. Zuzio, J.-L. Estivalèzes, and B. DiPierro.  An improved multiscale Eulerian-Lagrangian method for simulation of atomization process.  _Computers and Fluids_ , 176:285–301, 2018.