[3]Santanu Sinha

1] PoreLab, Department of Chemistry, Norwegian University of Science and
Technology (NTNU), Trondheim, N-7491, Norway

2] PoreLab, Department of Physics, Norwegian University of Science and
Technology (NTNU), Trondheim, N-7491, Norway

[3] PoreLab, Department of Physics, University of Oslo (UiO), Oslo, N-0371,
Norway

# Steady-state two-phase flow of compressible and incompressible fluids in a
capillary tube of varying radius

Hyejeong L. Cheon  Hursanay Fyhn  Alex Hansen  Øivind Wilhelmsen
<EMAIL_ADDRESS>[  [  *

###### Abstract

We study immiscible two-phase flow of a compressible and an incompressible
fluid inside a capillary tube of varying radius under steady-state conditions.
The incompressible fluid is Newtonian and the compressible fluid is an
inviscid ideal gas. The surface tension associated with the interfaces between
the two fluids introduces capillary forces that vary along the tube due to the
variation in the tube radius. The interplay between effects due to the
capillary forces and the compressibility results in a set of properties that
are different from incompressible two-phase flow. As the fluids move towards
the outlet, the bubbles of the compressible fluid grows in volume due to the
decrease in pressure. The volumetric growth of the compressible bubbles makes
the volumetric flow rate at the outlet higher than at the inlet. The growth is
not only a function of the pressure drop across the tube, but also of the
ambient pressure. Furthermore, the capillary forces create an effective
threshold below which there is no flow. Above the threshold, the system shows
a weak non-linearity between the flow rates and the effective pressure drop,
where the non-linearity also depends on the absolute pressures across the
tube.

###### keywords:

two-phase flow, compressibility, bubble-growth, rheology

##  1 Introduction

Hydrodynamic properties of the flow of multiple immiscible and incompressible
fluids, otherwise known as two-phase flow b88 ; d92 ; b17 ; ffh22 , are
controlled by a number of different factors: fluid properties such as the
viscosity contrast and surface tension between the fluids, driving parameters
such as the applied pressure drop or the flow rate, and geometrical properties
of the medium such as the size and shape of the pore space. The combined
effects of these factors make two-phase flow different and more complex than
single phase flow. The dimensionless parameters that play a key role to define
the flow properties are the ratio between the viscous and capillary forces,
referred to as the capillary number, and the ratio between the viscosities of
the two fluids. Depending on the values of these parameters, the flow
generates different types of fingering patterns cw85 ; lz85 ; mfj85 ; lmt04 ;
zmp19  or stable displacement fronts ltz88  during invasion processes where
one fluid displaces another in the porous medium.

Displacement processes are transient. If one continues to inject after
breakthrough, the flow enters a steady state characterized by a situation
where the macroscopic flow properties fluctuate or remain constant around
well-defined averages. A more general form of steady-state flow can be
achieved by continuously injecting both fluids simultaneously. In this case,
the dynamics at the pore scale might have fluid clusters breaking up and
forming, while the macroscopic flow parameters still have well-defined
averages.

Over the last decade, it has become clear that steady-state flow deviates from
the linear Darcy relationship d56  between the total flow rate and pressure
drop over a range of parameters. Rather, one finds a power law relationship
between pressure drop and the volumetric flow rate tkr09 ; rcs11  in that
range. In terms of the capillary number, this range is intermediary, with
linearity appearing both for lower and higher values sbd17 ; glb20 ; zbg21 .
Theoretical work to understand the physics behind the non-linearity has
appeared in e.g. tlk09 ; sh12 ; zbg21 , and computational studies have been
performed using Lattice Boltzmann simulations yts13  and dynamic pore network
modeling sgv21 ; sbd17 . It is now believed that a fundamental mechanism
behind this non-linearity is the capillary barriers at the pore throats, which
create an effective yield threshold. When the viscous forces increase, they
overcome the capillary barriers creating new flow paths. This increases the
effective mobility and thus the non-linear behavior appears rh87 . The
disorder in the pore-space properties, such as the pore-size distribution
rsh21  and the wetting angle distribution fsr21 , therefore play key roles in
determining the value of the exponent relating the volumetric flow rate and
the pressure drop in the non-linear regime.

The majority of the analytical and numerical approaches mentioned above
consider the two fluids to be incompressible, whereas many of the experiments
and applications use air as one of the fluids. Air is strongly compressible,
which introduces complex pore-scale mechanisms such as trapping and
coalescence l41 ; ly94 . Compressibility is relevant to a wide range of
applications with liquid and gas transport in porous media, for example, ${\rm
CO}_{2}$ transport and storage rk15 ; lkd15 ; ipr19  and the transport in fuel
cells nmn20 . Another class of applications where the compressibility plays a
key role are those involving phase transitions of the fluids such as boiling
and condensation. There are industrial applications where such processes are
of high importance, for example aerospace vehicle thermal protection hzl17 ,
high power electronics cooling systems gzk11 ; lww12 ; lwy20  and chemical
reactors bs19 . These applications utilize the high specific surface area of a
porous medium with fluid flowing inside, which enhances the heat and mass
transfer rates sgd16 ; szx11 . There are also natural processes such as drying
of soil rn94  where a liquid to gas transition takes place.

In this paper we present a study of two-phase flow of a mixture of
compressible and incompressible fluids in a capillary tube with varying
radius. We consider two fluids, one is an incompressible Newtonian fluid
obeying Poiseuille flow whereas the other is a compressible ideal gas, where
the viscosity is assumed to be negligible. The fluids flow as a series of
bubbles and droplets under a constant pressure drop along the tube.

In case of two-phase flow of two incompressible fluids in a corresponding
capillary tube, it has been found that the volumetric flow rate depends on the
square root of the pressure drop along the tube minus a threshold pressure
shb13 . The primary goal of the present work is to determine how this
constitutive equation changes when one of the two fluids is compressible.

A secondary goal of this work is to provide a basis for dynamic pore network
modeling b01 ; mt09 ; jh12 ; sgv21  of compressible-incompressible fluid
mixtures. This opens the possibility for incorporating thermodynamic effects
in such models such as boiling. We note that the other dominating
computational model in this context, the Lattice Boltzmann model grz91 ; rin12
, can only incorporate fluids that are weakly compressible qyz17 ; gfj20 .

We describe in Section 2.1 the equations that govern the flow through the
capillary tube. In Section 2.2 we introduce the boundary conditions used,
i.e., how we inject alternate compressible and incompressible fluid into the
tube. Section 2.3 describes how the governing equations are integrated in
time.

Section 3 presents the results of our investigation. Section 3.1 defines what
we mean by steady-state flow in the context of expanding bubbles. In section
3.2 we investigate how the compressible bubbles grow as they advance along the
tube, thus increasing the overall flow rate of the fluids. Section 3.3
presents the relation between volumetric flow rate and pressure drop at both
the inlet and outlet.

We summarize our results in Section 4. Section 5 contains the description of
the videos provided in the electronic supplementary material.

##  2 Methodology

The capillary tube considered in this work is filled with an incompressible
and a compressible fluid, immiscible to each other, which flow through it. The
radius of the capillary tube varies along the flow direction, $x$. The fluids
are separated by menisci, generating a surface tension. The incompressible
fluid is a viscous Newtonian liquid obeying Hagen-Poiseuille flow whereas the
compressible fluid is an inviscid ideal gas. The flow occurs as a plug flow
with a series of alternate bubbles and droplets of the two fluids as
illustrated in Figure 1. There is no fluid film along the tube walls and
therefore no coalescence or snap off taking place inside the tube during the
flow. We will refer the compressible and the incompressible fluid segments as
bubbles and droplets respectively.

Figure 1: Illustration of the tube geometry and the indexed variables. The
shaded fluid represents the non-wetting compressible fluid and the white fluid
represents the more wetting incompressible liquid. There are $N=6$ bubbles
here indicated by the numbers $i=1,\ldots,6$. The indexed variables $P_{i}$,
$V_{i}$ and $n_{i}$ respectively correspond to the pressure, volume and moles
of the $i$th bubble whereas $Q_{i}$ corresponds to the flow rate of the
droplet between $i$th and $(i+1)$th bubbles.

###  2.1 Governing equations

We assume that at a given time the system contains $N$ compressible bubbles
denoted by $i=1,2,\ldots,N$ from left to right as shown in Figure 1. The
volume $V_{i}$ and the pressure $P_{i}$ of the $i$th bubble are connected
through the ideal gas law,

$\displaystyle P_{i}V_{i}=n_{i}RT\;,$ (1)

where $n_{i}$ is the number of moles of gas present inside the bubble, $R$ is
the ideal gas constant and $T$ is the temperature. The tube is assumed to have
a constant average cross-sectional area ($A$) in terms of the fluid volume and
therefore $V_{i}=A(x_{i}^{r}-x_{i}^{l})$ where $x_{i}^{l}$ and $x_{i}^{r}$ are
the positions of the left and right menisci of the $i$th bubble respectively.

The volume of an incompressible droplet on the other hand will remain constant
throughout the flow and the flow rate will depend on the pressures of the two
compressible bubbles bordering it. The volumetric flow rate of the
incompressible droplet between $i$ and $i+1$ is denoted by $Q_{i}$, and
follows the constitutive equation w21 ,

$\displaystyle
Q_{i}=\frac{A^{2}}{8\pi\mu(x_{i+1}^{l}-x_{i}^{r})}\left[P_{i}-P_{\rm
c}(x_{i}^{r})-P_{i+1}+P_{\rm c}(x_{i+1}^{l})\right]\;,$ (2)

where $\mu$ is the viscosity of the incompressible fluid and $P_{c}(x)$ is the
capillary pressure at $x$. Here we consider the incompressible fluid to be
more wetting with respect to the pore walls than the compressible fluid, thus
determining the sign of $P_{c}$ in Equation 2. We model $P_{c}$ by using the
Young-Laplace equation d92 ,

$\displaystyle P_{\rm{c}}(x)=\frac{2\gamma}{r(x)}\;,$ (3)

where $r(x)$ is the radius of the tube at $x$. Here
$\gamma=\sigma\cos(\theta)$ where $\sigma$ is the surface tension between the
fluids and $\theta$ is the wetting angle of the fluid with respect to the tube
wall. The variation in the radius of the tube shown in Figure 1 is modeled by

$\displaystyle r(x)=\frac{1}{2}\left[w+2a\cos\left(\frac{2h\pi
x}{L}\right)\right]$ (4)

where $L$ is the tube length, $w$ is the average radius, $a$ is the amplitude
of oscillation and $h$ is the number of periods.

###  2.2 Boundary conditions

The system is driven by a constant pressure drop $\Delta P=P_{0}-P_{L}$ where
$P_{0}$ and $P_{L}$ are the pressures at the inlet ($x=0$) and outlet ($x=L$)
respectively. The two fluids are injected alternatively at the inlet.
Depending on which fluid that is being injected and which fluid that is
leaving the tube, there will be different configurations as illustrated in
Figure 2. When a bubble is entering at the inlet [Figure 2(a)] or leaving at
the outlet [Figure 2(c)], the pressure in that bubble is given by $P_{0}$ or
$P_{L}$ respectively. This is because the compressible fluid has no viscosity
and thus the pressure inside a bubble is uniform. The pressures inside all
other bubbles are calculated using Equation 1. When a droplet is entering at
the inlet [Figure 2(b) and (c)] or leaving at the outlet [Figure 2(a) and
(b)], the respective flow rates $Q_{0}$ and $Q_{N}$ are given by,

$\begin{split}\displaystyle Q_{0}&=\frac{A^{2}}{8\pi\mu
x_{1}^{l}}\left[P_{0}-P_{1}+P_{\rm c}(x_{1}^{l})\right]\;{\rm and}\\\
Q_{N}&=\frac{A^{2}}{8\pi\mu(L-x_{N}^{r})}\left[P_{N}-P_{\rm
c}(x_{N}^{r})-P_{L}\right]\;,\end{split}$ (5)

whereas the flow rates of the remaining droplets are calculated using Equation
2.

(a) (b) (c)

Figure 2: Illustration of different configurations where bubbles and droplets
are colored as gray and white respectively. In (a), a bubble is entering at
the inlet and therefore $P_{1}=P_{0}$ there. In (c), a bubble is leaving at
the outlet, therefore $P_{N}=P_{L}$ there. A droplet is entering at the inlet
in (b) and (c), and leaving at the outlet in (a) and (b). The flow rates of
such droplets are calculated using Equation 5.

The simulation is started with the tube completely filled with the
incompressible fluid. The two fluids are then injected alternately through the
inlet using small time steps. Whenever the injection is switched to a
different fluid, a new menisci is created and the injection is continued for
that fluid until the bubble or the droplet being injected has reached a given
length, $b_{\rm C}$ or $b_{\rm I}$ respectively. For each new bubble or
droplet, a new value for $b_{\rm C}$ or $b_{\rm I}$ is determined using the
following scheme:

$\displaystyle b_{\rm C}=b_{\rm min}+kF_{\rm C}b_{\rm max}\quad\;{\rm
and}\quad b_{\rm I}=b_{\rm min}+kF_{\rm I}b_{\rm max}\;,$ (6)

where $k$ is chosen from a uniform distribution of random numbers between $0$
and $1$. $F_{\rm C}$ and $F_{\rm I}$ are the tentative values of the
fractional flows for the bubbles and droplets respectively. The two parameters
$b_{\rm max}$ and $b_{\rm min}$ set the smallest and largest allowed sizes of
any bubble or droplet. We consider here $b_{\rm min}=L/10^{4}$ and $b_{\rm
max}=L/50$. The parameters $b_{\rm C}$ and $b_{\rm I}$ decide the initial
sizes of the bubbles and droplets just after they detach from the inlet. For
the compressible fluid, this determines the number of moles $n_{i}$ inside a
bubble,

$\displaystyle n_{i}=\frac{Ab_{\rm C}P_{0}}{RT}\;,$ (7)

which remains constant for that bubble throughout the flow after it gets
detached from the inlet.

###  2.3 Updating the menisci

At any time, the two menisci bordering a droplet inside the tube move with the
same velocities. The velocities of the menisci are calculated from the
velocities $v_{i}$ of the droplets using Equations 2 and 5,

$\displaystyle\dfrac{\operatorname{d\\!}{}{x_{i}^{r}}}{\operatorname{d\\!}{t}}=\dfrac{\operatorname{d\\!}{}{x_{i+1}^{l}}}{\operatorname{d\\!}{t}}=v_{i}=\frac{Q_{i}}{A}\;.$
(8)

We solve these ordinary differential equations using an explicit Euler scheme,
thus updating positions of all menisci by choosing a small time step $\Delta
t$.

Depending on the position of the menisci and the corresponding capillary
pressures, the bubbles may compress or expand. If a bubble compresses at any
time step, it means the left and right interfaces of that bubble approach each
other. This necessitates the choice of time step $\Delta t$ to be sufficiently
small, as otherwise, the two menisci around that bubble will collapse after
the time step. We deal with this situation in the following way. First we
calculate a time $\Delta t_{1}$ that is needed to pass one pore-volume of
incompressible fluid through the tube,

$\displaystyle\Delta t_{1}=\frac{8\pi\mu L^{2}}{A(P_{0}-P_{L})}\;.$ (9)

Next, we check for every bubble $i$ if $(v_{i-1}-v_{i})>0$, that is, whether
the two menisci bordering the bubble are approaching each other in that time
step. If this criterion is found to be true for any of the bubbles $j$, we
measure the time it will take for the two menisci to collapse,

$\displaystyle\Delta t_{2}^{j}=\frac{x_{j}^{r}-x_{j}^{l}}{v_{j-1}-v_{j}}\;.$
(10)

After calculating $\Delta t_{1}$ and $\Delta t_{2}^{j}$, we determine a time
$\Delta t$ for that step from,

$\displaystyle\Delta t={\rm min}(a^{*}\Delta t_{1},b^{*}\Delta t_{2}^{j})\;,$
(11)

which means that if there is a possibility for a bubble to collapse during the
time step, we chose $\Delta t$ from the minimum of $a^{*}\Delta t_{1}$ and all
of $b^{*}\Delta t_{2}^{j}$. If there is no possibility of collapse, we use
$\Delta t$ equal to $a^{*}\Delta t_{1}$. For the simulations presented in this
paper, we set $a^{*}=10^{-8}$ and $b^{*}=10^{-6}$.

##  3 Results and discussions

We perform steady-state simulations considering a tube of length $L=100\,{\rm
cm}$ with $w=1\,{\rm cm}$, $a=0.25\,{\rm cm}$ and $h=30$ (Equation 4). The
viscosity of the incompressible fluid is $\mu=0.001\,{\rm Pa.s}$, the ideal
gas constant is $R=8.31\,{\rm J/(mol.K)}$ and the temperature is kept fixed
throughout the simulation at $T=293\,{\rm K}$. We fix $F_{c}=0.4$ (Equation 6)
which sets the volumetric fractional flow of the compressible fluid at the
inlet around that value. We perform simulations varying the pressure drops
($\Delta P=P_{0}-P_{L}$) as well as the absolute outlet pressure with
different values of the surface tension, $\gamma$.

Figure 3: Total volumetric flow rate $q_{\rm T}^{\rm i}$ at the inlet as a
function of the injected pore volume $V_{p}$ for the outlet pressure
$P_{L}=1\,{\rm kPa}$ and the surface tension $\gamma=0.09\,{\rm N/m}$ and for
the pressure drops $\Delta P=1,2,4$ and $8\,{\rm kPa}$ respectively. The
steady-state values of the flow rates are measured by taking averages in the
range of $20$ to $40$ pore volumes as indicated by the dashed lines.

###  3.1 Steady-state flow

The steady state is defined by the volumetric flow rates of the fluids
fluctuating around a stable average. Due to the expansion of the compressible
fluid, which we will discuss in a moment, the volumetric flow rate of the
fluids changes as the fluids flow towards the outlet. We define the quantities
$Q_{\rm T}^{\rm i}$, $Q_{\rm C}^{\rm i}$, $Q_{\rm I}^{\rm i}$ as the average
steady-state flow rates for the total, compressible and incompressible fluids
at the inlet and $Q_{\rm T}^{\rm o}$, $Q_{\rm C}^{\rm o}$, $Q_{\rm I}^{\rm o}$
as those at the outlet. The inlet and outlet flow rates are measured by
tracking the displacement of the first meniscus nearest to the inlet and the
last meniscus near the outlet, which are either the left or the right meniscus
of the first ($i=1$) and the last ($i=N$) bubbles. The instantaneous flow
rates of the bubbles and droplets are measured as $q_{\rm C}^{\rm
i}=A\sum\Delta x_{1}^{r}/\sum\Delta t$ for $x_{1}^{l}=0$, $q_{\rm I}^{\rm
i}=A\sum\Delta x_{1}^{l}/\sum\Delta t$ for $x_{1}^{l}>0$ and $q_{\rm C}^{\rm
o}=A\sum\Delta x_{N}^{l}/\sum\Delta t$ for $x_{N}^{l}=L$, $q_{\rm I}^{\rm
o}=A\sum\Delta x_{N}^{r}/\sum\Delta t$ for $x_{N}^{R}<L$. This measurement is
performed after every $0.05$ pore-volumes of fluid are injected and the sum is
therefore over the time steps in between. The total flow rates are therefore
given by, $q_{\rm T}^{\rm i,o}=q_{\rm C}^{\rm i,o}+q_{\rm I}^{\rm i,o}$. This
provides the measurement of the injected and outlet flow rates as a function
of the injected pore volumes or of the time. In Figure 3, we plot $q_{\rm
T}^{\rm i}$ as a function of the pore-volumes ($V_{p}$) injected for
$P_{L}=1\,{\rm kPa}$. The pore-volume $V_{p}$ is defined as the ratio between
the volume of the inject fluids and the volume of the total pore space of the
tube, which provides an estimate of how many times the pore space was flushed
with the fluids. The plots show that $q_{\rm T}^{\rm i}$ increases with time
at the beginning of the flow. This increase in $q_{\rm T}^{\rm i}$ is due to
the decrease in the effective viscosity of the system caused by the injection
of inviscid compressible gas into the tube filled with viscous incompressible
fluid. After the injection of a few pore volumes, $q_{\rm T}^{\rm i}$
fluctuates around a constant average ($Q_{\rm T}^{\rm i}$) shown by the
horizontal dashed lines which defines the steady state. We run our simulations
for $40$ pore volumes of fluid where the steady-state averages are taken after
$20$ pore volumes injected to ensure that a steady state has been reached.

###  3.2 Bubble growth

As a compressible bubble moves along the tube, the volume of the bubble
increases due to the decrease in the pressure towards the outlet vls10 . The
bubble can also grow due to other mechanisms, such as the increase in
temperature or a phase transition between liquid and gas phases w98 ; kwd06 ,
but these phenomena are not studied here. A simulation with a single bubble
inside a short tube is shown in the supplementary material which illustrates
that the bubble increases in size as it flows towards the outlet. To
understand how this growth depends on different flow parameters in the steady
state, we define the growth function $G_{\rm C}(x)$ by,

$\displaystyle G_{\rm C}(x)=\frac{V(x)-V_{\rm 0}}{V_{0}}\;,$ (12)

where $V_{0}$ and $V(x)$ are the volume of a given bubble initially after
detaching from the inlet and when its center is at $x$. We measure $G_{\rm C}$
by including all the bubbles that are not attached to the inlet or outlet and
calculate the time average value of $(V(x)-V_{0})/V_{0}$ in the investigated
time interval, where $x$ is the center of the bubble.

(a) (b)

Figure 4: Plot of the bubble growth $G_{\rm C}(x)$ in the steady state as a
function of the scaled position $x/L$ inside the tube for zero surface
tension, $\gamma=0$. The two plots show the results for the same set of
pressure drops $\Delta P$ with different outlet pressures $P_{L}$.

(a) (b)

Figure 5: Variation of the pressure $P(x)$ [kPa] inside a compressible bubble
along the tube during steady state flow. $P(x)$ shows a linear behavior for
different values of $\Delta P$ and $P_{L}$.

Figure 4 shows the variation of $G_{\rm C}(x)$ along the tube for two
different outlet pressures, $P_{L}=1\,{\rm kPa}$ and $100\,{\rm kPa}$ where we
plot the results for the same set of pressure drops $\Delta P$. These results
are with zero surface tension, $\gamma=0$. There are a few details to note
here. First, the plots show that $G_{\rm C}(x)$ increases with an increase in
$\Delta P$. In addition, $G_{\rm C}(x)$ also depends on the absolute pressures
at the inlet and outlet, since we can see that the curves are non-linear
functions of $x$ for $P_{L}=1\,{\rm kPa}$, whereas for $P_{L}=100\,{\rm kPa}$,
they show linear behavior. Furthermore, $G_{\rm C}(x)$ approaches $\Delta
P/P_{L}$ at $x=L$ for all the data sets.

To explain the dependency of $G_{\rm C}(x)$ on $\Delta P$ and $P_{L}$, we
recall Equation 1 and rewrite Equation 12 as,

$\displaystyle G_{\rm C}(x)=\frac{P_{0}-P(x)}{P(x)}\;,$ (13)

where $P(x)$ is the pressure inside a bubble at $x$. For $x=L$, $P(x)=P_{L}$
and therefore $G_{\rm C}(L)=\Delta P/P_{L}$ as observed. In Figure 5, we plot
$P(x)$, averaged over different time steps in the steady state, for the two
outlet pressures, $P_{L}=1\,{\rm kPa}$ and $100\,{\rm kPa}$. Both of the plots
show linear variation along $x$ with the slope $-\Delta P$. We therefore have
$P(x)=-x\Delta P/L+P_{L}+\Delta P$ and thus,

$\displaystyle\frac{G_{\rm C}(x)}{n_{P}}=\frac{x/L}{1+n_{P}(1-x/L)}\;,$ (14)

where $n_{P}=\Delta P/P_{L}$. This leads to

$\displaystyle\frac{G_{\rm
C}(x)}{n_{P}}\sim\begin{cases}\displaystyle\frac{1}{n_{P}}\left(\frac{1}{1-x/L}-1\right)&{\rm
for}\quad n_{P}\gg 1\;,\\\ x/L&{\rm for}\quad n_{P}\ll 1\;,\end{cases}$ (15)

which explains the concave and linear variation of $G_{C}$ as function of
$x/L$ observed in Figure 4 (a) and (b) respectively. The growth of the bubbles
along the tube is therefore a function of $n_{P}=\Delta P/P_{L}$.

(a) (b)

Figure 6: Variation of the bubble growth $G_{\rm C}$, scaled with
$n_{P}=\Delta P/P_{L}$, with $x/L$. Results are plotted for the same sets of
$n_{P}$ for two different values of $P_{L}$. The left and right figures
correspond to $\gamma=0$ and $0.03\,{\rm N/m}$ respectively. In each plot, the
line corresponds to $P_{L}=1\,{\rm kPa}$ and the symbols correspond to
$P_{L}=100\,{\rm kPa}$.

In Figure 6 we plot $G_{\rm C}/n_{P}$ for the two outlet pressures $P_{L}$
with the same sets of values of $n_{P}$ for (a) $\gamma=0$ and (b)
$\gamma=0.3\,{\rm N/m}$. The plots show that the results for the same values
of $n_{P}$ follow the same curves, irrespective of the outlet pressures
$P_{L}$. Furthermore, for the non-zero surface tension case in Figure 6 (b),
$G_{\rm C}$ also shows a periodic oscillation along $x$ when both the $n_{P}$
and $P_{L}$ are small, that is, for $P_{L}=1\,{\rm kPa}$ and $n_{P}\leq 1$. In
addition, there is no data point for $n_{P}\leq 0.3$ with $P_{L}=1\,{\rm
kPa}$, as the movement of the bubbles stopped due to high capillary barriers.
This suggests the existence of an effective threshold pressure, below which
there will be no flow through the tube. This threshold depends on both
$\gamma$ and $P_{L}$, which we will explore more in the following section. We
show the different characteristics of flow in the videos provided in
electronic supplementary material.

###  3.3 Effective rheology

Equations 1 and 2 resist analytical solutions even in the case when there is
only a single compressible bubble in the tube. This is due to the pressure in
the compressible bubble being inversely proportional to the difference in
position of the two menisci surrounding it, whereas the motion of the two
surrounding incompressible fluids is determined by the cosine of the positions
of the same menisci. These equations, even in this simplest case, are
therefore highly non-linear with an essential singularity lurking in the very
neighborhood where we seek solutions. We therefore stick to numerical analysis
in the following.

Figure 7:  Plot of the flow rates for the total ($Q_{\rm T}^{\rm i,o}$),
compressible ($Q_{\rm C}^{\rm i,o}$) and incompressible ($Q_{\rm I}^{\rm
i,o}$) fluids at the inlet (left column) and at the outlet (right column) for
$P_{L}=1\,{\rm kPa}$ as a function of $\Delta P$. The different sets in each
plot correspond to different values of the surface tension indicated in the
legends. The dashed line in each plot has a slope of $1$.

Due to the volumetric growth of the compressible bubbles during their flow
towards the outlet, the volumetric flow rate varies along the tube. In
addition, this volumetric growth is a function of the pressures, making the
average saturation and the effective viscosity of the two fluids inside the
tube pressure dependent. These two mechanisms together control the effective
rheological behavior of the steady-state flow. In Figure 7 we show the
variation of the volumetric flow rates ($Q_{\rm T}^{\rm i,o}$, $Q_{\rm C}^{\rm
i,o}$, $Q_{\rm I}^{\rm i,o}$) as functions of the pressure drop $\Delta P$ for
the outlet pressure $P_{L}=1\,{\rm kPa}$ and for different values of the
surface tension ($\gamma$). Note the differences between the inlet and outlet
flow rates for the total and for the each component of flow. For the
incompressible fluid, there is no significant increase in the outlet flow rate
compared to its inlet flow rate (third row in Figure 7) whereas there is a
noticeable increase in the outlet flow rate of the compressible fluid (second
row in Figure 7). This increase in $Q_{\rm C}^{\rm o}$ effectively increases
the total flow rate at the outlet (first row in Figure 7). The dashed line in
Figure 7 has a slope equal to $1$. The total flow rates show deviations from
this dashed line. For the inlet, $Q_{\rm T}^{i}$ shows small deviations from
the dashed line for $\gamma>0$ at small $\Delta P$. Whereas at the outlet, the
deviations are significantly higher due to the increase in the volumetric
growth of the compressible fluid.

Another point to note in Figure 7 is that there is a minimum value of $\Delta
P$, below which there is no data point available. This is due to the existence
of a threshold pressure below which the flow stops. In the supplementary
material we show a simulation video in this regime where one can observe that
the flow of the bubbles stops at a certain time step. The threshold is due to
the capillary forces at the menisci between the two fluids that create
capillary barriers at the narrowest points along the tube. Such threshold was
also observed in the case of two-phase flow of two incompressible fluids in a
tube with variable radius shb13 . There, it was shown analytically that the
average flow rate $Q$ in the steady state varies with the applied pressure
drop $\Delta P$ as, $Q\sim\sqrt{\Delta P^{2}-P_{\rm t}^{2}}$ where $P_{\rm t}$
is the effective threshold pressure. When $\mathinner{\\!\left\lvert\Delta
P\right\rvert}-P_{\rm t}\ll P_{\rm t}$, this relationship translates to
$Q\sim\sqrt{\mathinner{\\!\left\lvert\Delta P\right\rvert}-P_{\rm t}}$, that
is, the flow rate varies with the excess pressure drop to the power of $0.5$.
The threshold pressure depends on the surface tension and on the configuration
of the menisci positions inside the tube. If the total capillary barrier is
higher than the applied pressure drop, the flow stops. This is similar here
for the two-phase flow with one of the fluids being compressible.

(a) (b) (c)

Figure 8: Plot of the volumetric inlet flow rate $Q_{\rm T}^{\rm i}$ as a
function of the excess pressure drop $(\Delta P-P_{\rm t})$ for $P_{L}=1\,{\rm
kPa}$ and $100\,{\rm kPa}$, where the values of $P_{\rm t}$ are obtained from
a minimization of the least square fit error $\epsilon$. The minimization is
illustrated in the insets of (a) and (b) for $\gamma=0.05\,{\rm N/m}$ (green)
and $0.09\,{\rm N/m}$ (purple). The dashed lines in (a) and (b) have a slope
$1$ whereas in (c), the lower and upper dashed lines have slopes $1$ and $1.3$
respectively.

We assume a general relation between the average volumetric flow rates $Q_{\rm
T}^{\rm i,o}$ and the pressure drop $\Delta P$ as,

$\displaystyle Q_{\rm T}^{\rm i,o}\sim(\Delta P-P_{\rm t})^{\beta_{\rm i,o}}$
(16)

where $\beta_{\rm i,o}$ is the corresponding exponent. In order to find both
the effective threshold pressure $P_{T}$ and the exponent $\beta_{\rm i,o}$
from the measurements of $Q_{\rm T}^{\rm i,o}$, we adopt an error minimization
technique that was used in earlier studies sh12 ; fsr21 . There we choose a
series of trial values for $P_{\rm t}$ and calculate the mean square error
$\epsilon$ for the linear least square fit by fitting the data points with
$\log(Q)\sim\log(\Delta P-P_{\rm t})$. Then we select the value of $P_{\rm t}$
that corresponds to the minimum value of $\epsilon$, implying the best fit of
the data points with Equation 16. This is illustrated in the insets of Figure
8 (a) and (b). The slope for the selected threshold $P_{\rm t}$ provides the
exponent $\beta_{\rm i,o}$. The variation of the total inlet and outlet flow
rates $Q_{\rm T}^{\rm i,o}$ with the excess pressure drop $(\Delta P-P_{\rm
t})$ are plotted in Figure 8 for the two outlet pressures $P_{L}=1$ and
$100\,{\rm kPa}$. The data sets show agreement with Equation 16 with the
selected values of $P_{\rm t}$ and $\beta$. There is a noticeable difference
between the slopes for the inlet and outlet flow rates for $P_{L}=1\,{\rm
kPa}$ whereas for $P_{L}=100\,{\rm kPa}$ they are similar. For
$P_{L}=100\,{\rm kPa}$ the data points for both $Q_{\rm T}^{\rm i}$ and
$Q_{\rm T}^{\rm o}$ follow a slope of $\approx 1.0$ whereas for $P_{L}=1\,{\rm
kPa}$, the data points for $Q_{\rm T}^{\rm i}$ and $Q_{\rm T}^{\rm o}$ follow
the slopes of $\approx 1.0$ and $1.3$ respectively. These are indicated by the
dashed lines in the figures.

(a) (b) (c)

Figure 9: Variation of the threshold pressure $P_{\rm t}$ and the exponents
$\beta_{\rm i,o}$ as functions of the effective surface tension $\gamma$ for
$P_{L}=1$ and $100\,{\rm kPa}$. $P_{\rm t}$ increases with the increase of
$\gamma$ and the values are much higher for $P_{L}=1\,{\rm kPa}$ compared to
$P_{L}=100\,{\rm kPa}$. The exponent $\beta_{\rm i}$ for the inlet flow rate
are close to $1$ for both the values of $P_{L}$ whereas for the outlet flow
rate $\beta_{\rm o}\approx 1.3$ for $P_{L}=1\,{\rm kPa}$. For $P_{L}=100\,{\rm
kPa}$, $\beta_{\rm o}$ remains close to $\beta_{\rm i}$. The dashed horizontal
lines indicate the value $1.0$ of the $y$ axis.

The variations of $P_{\rm t}$ and $\beta_{\rm i,o}$ with the surface tension
$\gamma$ are plotted in Figure 9. The data points were calculated by
considering different ranges of $\Delta P$ and taking averages over the
ranges, and the corresponding standard deviations are plotted as error bars.
The threshold pressure $P_{\rm t}$ is zero at $\gamma=0$ and then increases
gradually with $\gamma$ which shows that the threshold appears due to
capillary forces. The increase in $P_{\rm t}$ with $\gamma$ appears to be
linear here which is similar to the case of two incompressible fluids, where
the linear dependence of $P_{\rm t}$ on the surface tension was shown
analytically shb13 . Additionally for the compressible flow here, the
thresholds also depend on the outlet pressure $P_{L}$. For the lower outlet
pressure $P_{L}=1\,{\rm kPa}$, the thresholds are systematically higher
compared to those for $P_{L}=100\,{\rm kPa}$ for the whole range of $\gamma$.
Furthermore, the exponents $\beta_{\rm i,o}$ also depend on the outlet
pressure as seen from the Figures 9 (b) and (c). The difference is more
visible for the exponents related to the outlet flow rates than the inlet. For
the inlet flow rate, $\beta_{\rm i}$ has values around $\approx 0.95$ and
$1.02$ for $P_{L}=1\,{\rm kPa}$ and $100\,{\rm kPa}$ respectively, showing
almost linear dependence for both the cases. For the outlet flow rates,
$\beta_{\rm o}$ remains close to $\beta_{\rm i}$ for $P_{L}=100\,{\rm kPa}$
whereas for $P_{L}=1\,{\rm kPa}$, $\beta_{\rm o}$ increases to $\approx 1.3$.
This increase in $\beta_{\rm o}$ compared to $\beta_{\rm i}$ reflects the
dependence of the volumetric growth $G_{\rm C}(x)$ of the bubbles on $P_{L}$,
indicating an underlying dependence of the rheological behavior on the
absolute inlet or outlet pressures. However, at this point we are unable to
describe how the two parameters $P_{\rm t}$ and $\beta$ scale with $P_{L}$,
which needs further study.

Existing studies of the power-law volumetric flow rate-pressure drop relation
for porous media have shown the existence of different regimes characterized
by different exponents. These studies involve experiments tkr09 ; tlk09 ;
rcs11 ; sbd17 ; glb20 ; zbg21 , Lattice Boltzmann simulations yts13 , pore-
network modeling sh12 ; sbd17  and analytical calculations tlk09 ; sh12 ;
rhs19 ; zbg21 . There are three regimes, an intermediate non-linear regime
where the flow rate $Q$ increases at a rate much faster than the applied
pressure drop $\Delta P$ with a power law exponent larger than one and up to a
value around 2.5. There are in addition two linear regimes for either smaller
yts13 ; glb20 ; zbg21  or larger yts13 ; sh12 ; sbd17  volumetric flow rates
than the non-linear regime. This allows the definition of a lower and upper
crossover pressure drop. A simple explanation for these three regimes may be
drawn from the study of the conductivity of a disordered network of threshold
resistors rh87  which is based on the following idea. Each resistor has a
threshold voltage to start conducting the current and then the current
increases linearly. In a network with links of such flow properties, there
will be a regime when new conducting paths will appear when increasing the
global pressure drop. The increase in the flow rate through each path together
with the increase in the number of paths leads to an effective increase of $Q$
faster than $\Delta P$. This results in the non-linear exponent being higher
than $1$, the value of which depends on the distribution of the thresholds in
each link rhs19 . The linear regime above this non-linear regime appears from
all the available paths being conducting whereas the linear regime below
appears from the flow being flow in single percolating channels without
interfaces. With this idea, the experimental glb20 ; zbg21  and numerical
yts13  observations of two-phase flow in porous media showing linear variation
of flow rate in this low pressure regime therefore indicate that the flow in
single channels consisting of many pores are linear, which is similar to what
we have found for the lower outlet pressure in the present
compressible/incompressible flow case.

##  4 Conclusions

We have studied the flow of alternating compressible bubbles and
incompressible droplets through a capillary tube with variable radius. The
motion of the bubbles was given by the model Equations 1 and 2, thus assuming
the compressible fluid to be an ideal gas with zero viscosity, whereas the
incompressible fluid is Newtonian. The incompressible fluid is more wetting
than the compressible gas, but not to a degree that films form.

We switch between injecting the compressible and incompressible fluid at
intervals so that the fractional flow rate is essentially constant at the
inlet. We fix a pressure drop along the tube in addition to an ambient
pressure. This creates steady-state flow conditions in the tube.

The compressible bubbles expand as they move from the higher pressure region
at the inlet towards the lower pressure at the outlet. This expansion
accelerates the incompressible fluid, thus making the volumetric flow rate
larger at the outlet than at the inlet. The lower the ambient pressure is, the
stronger this effect is.

We measure volumetric flow rate at the inlet, finding essentially a linear
relationship between the volumetric flow rate and the pressure drop. However,
there is a threshold pressure that needs to be overcome in order to have flow
through the tube.

At the outlet, we find that the volumetric flow rate is still linear in the
excess pressure drop when the ambient pressure is low. However, when the
ambient pressure is high, the volumetric flow rate at the outlet becomes
proportional to the excess pressure to a power of around $1.3$.

This behavior is very different from that of two incompressible fluids moving
through a corresponding tube: Here the volumetric flow rate, being the same at
the inlet and the outlet, is proportional to the square root of the excess
pressure.

We expected the flow rate-pressure drop constitutive relations to be different
in this compressible/incompressible case than that of two incompressible
fluids. However, that we should find linearity was a big surprise. A precise
explanation as to why this is so, is still lacking.

Besides these surprising results, this work makes a first step in implementing
the modeling of compressible/incompressible fluid mixtures in dynamic network
models. We may then envision using more sophisticated equations of state for
the compressible fluid beyond the ideal gas law. This allows the consideration
of e.g. phase transitions such as boiling and condensation in porous media.

##  5 Supplementary material

The electronic supplementary material contains videos showing different flow
characteristics. The videos can be found in the list of ancillary files in the
arXiv abstract page of this article. We considered a tube with $L=10\,{\rm
cm}$, $w=1\,{\rm cm}$, $a=0.25\,{\rm cm}$ and $h=5$ (Equation 4). The
simulations were performed for $P_{L}=1\,{\rm kPa}$ and $\gamma=0.2\,{\rm
N/m}$. The compressible bubbles are colored with magenta whereas the
incompressible droplets are colored with black. The videos are not in real
time. We show four different simulations with different values of $\Delta P$:

  1. (a)

Flow of a single bubble of compressible gas in incompressible fluid. Here
$\Delta P=5\,{\rm kPa}$. The video shows the increase in the volume of the
bubble as it approaches towards the outlet.

  2. (b)

Injection of multiple compressible bubbles and incompressible droplets at a
very low pressure drop, $\Delta P=0.3\,{\rm kPa}$. The flow stops after a
certain time when several interfaces appeared in the tube. This shows the
existence of a total capillary barrier, which is higher than the applied
pressure drop here.

  3. (c)

Two-phase flow of multiple compressible bubbles and incompressible droplets at
a low pressure drop, $\Delta P=0.4\,{\rm kPa}$. Here the bubbles speed up and
slow down as they flow, showing the combined effect of the surface tension and
the shape of the tube. The bubbles also grow in volume towards the outlet.

  4. (d)

Two-phase flow of multiple compressible bubbles and incompressible droplets at
a higher pressure drop, $\Delta P=3\,{\rm kPa}$. The bubbles do not show any
significant slowing down in this case, indicating the capillary forces being
negligible compared to the viscous pressure drop. The volumetric expansion of
the compressible bubbles can also be observed here.

## Acknowledgments

We thank Federico Lanza, Marcel Moura and Håkon Pedersen for helpful
discussions. This work is supported by the Research Council of Norway through
its Centers of Excellence funding scheme, project number 262644.

## References

  *   * (1) Bear, J.: Dynamics of Fluids in Porous Media. Dover, Mineola, New York (1988)
  * (2) Dullien, F.A.L.: Porous Media: Fluid, Transport and Pore Structure. Academic Press, San Diego (1992)
  * (3) Blunt, M.J.: Multiphase Flow in Permeable Media. Cambridge University Press, Cambridge (2017)
  * (4) Feder, J., Flekkøy, E.G., Hansen, A.: Physics of Flow in Porous Media. Cambridge University Press, Cambridge (2022)
  * (5) Chen, J.D., Wilkinson, D.: Pore-scale viscous fingering in porous media. Phys. Rev. Lett. 55, 1892 (1985). https://doi.org/10.1103/PhysRevLett.55.1892
  * (6) Lenormand, R., Zarcone, C.: Invasion percolation in an etched network: Measurement of a fractal dimension. Phys. Rev. Lett. 54, 2226 (1985). https://doi.org/10.1103/PhysRevLett.54.2226
  * (7) Måløy, K.J., Feder, J., Jøssang, T.: Viscous fingering fractals in porous media. Phys. Rev. Lett. 55, 2688 (1985). https://doi.org/10.1103/PhysRevLett.55.2688
  * (8) Løvoll, G., Méheust, Y., Toussaint, R., Schmittbuhl, J., Måløy, K.J.: Growth activity during fingering in a porous hele-shaw cell. Phys. Rev. E 70, 026301 (2004). https://doi.org/10.1103/PhysRevE.70.026301
  * (9) Zhao, B., MacMinn, C.W., Primkulov, B.K., Chen, Y., Valocchi, A.J., Zhao, J., et al: Comprehensive comparison of pore-scale models for multiphase flow in porous media. Proc. Natl. Acad. Sci. USA. 116, 13799 (2019). https://doi.org/10.1073/pnas.1901619116
  * (10) Lenormand, R., Touboul, E., Zarcone: Numerical models and experiments on immiscible displacements in porous media. J. Fluid Mech. 189, 165 (1988). https://doi.org/10.1017/S0022112088000953
  * (11) Darcy, H. Les Fontaines publiques de la ville de Dijon 647, (1856)
  * (12) Tallakstad, K.T., Knudsen, H.A., Ramstad, T., Løvoll, G., Måløy, K.J., Toussaint, R., Flekkøy, E.G.: Steady-state two-phase flow in porous media: Statistics and transport properties. Phys. Rev. Lett. 102, 074502 (2009). https://doi.org/10.1103/PhysRevLett.102.074502
  * (13) Rassi, E.M., Codd, S.L., Seymour, J.D.: Nuclear magnetic resonance characterization of the stationary dynamics of partially saturated media during steady-state infiltration flow. New J. Phys. 13, 015007 (2011). https://doi.org/10.1088/1367-2630/13/1/015007
  * (14) Sinha, S., Bender, A.T., Danczyk, M., Keepseagle, K., Prather, C.A., Bray, J.M., Thrane, L.W., Seymour, J.D., Codd, S.L., Hansen, A.: Effective rheology of two-phase flow in three-dimensional porous media: experiment and simulation. Transp. Porous Med. 119, 77 (2017). https://doi.org/10.1007/s11242-017-0874-4
  * (15) Gao, Y., Lin, Q., Bijeljic, B., Blunt, M.J.: Pore-scale dynamics and the multiphase darcy law. Phys. Rev. Fluids. 5, 013801 (2020). https://doi.org/10.1103/PhysRevFluids.5.013801
  * (16) Zhang, Y., Bijeljic, B., Gao, Y., Lin, Q., Blunt, M.J.: Quantification of non-linear multiphase flow in porous media. Geophys. Res. Lett. 48, 2020–090477 (2021). https://doi.org/10.1029/2020GL090477
  * (17) Tallakstad, K.T., Løvoll, G., Knudsen, H.A., Ramstad, T., Flekkøy, E.G., Måløy, K.J.: Steady-state, simultaneous two-phase flow in porous media: An experimental study. Phys. Rev. E. 80, 036308 (2009). https://doi.org/10.1103/PhysRevE.80.036308
  * (18) Sinha, S., Hansen, A.: Effective rheology of immiscible two-phase flow in porous media. Europhys. Lett. 99, 44004 (2012). https://doi.org/10.1209/0295-5075/99/44004
  * (19) Yiotis, A.G., Talon, L., Salin, D.: Blob population dynamics during immiscible two-phase flows in reconstructed porous media. Phys. Rev. E. 87, 033001 (2013). https://doi.org/10.1103/PhysRevE.87.033001
  * (20) Sinha, S., Gjennestad, M.A., Vassvik, M., Hansen, A.: Fluid meniscus algorithms for dynamic pore-network modeling of immiscible two-phase flow in porous media. Front. Phys. 9, 548497 (2021). https://doi.org/10.3389/fphy.2020.548497
  * (21) Roux, S., Herrmann, H.J.: Disorder-induced nonlinear conductivity. Europhys. Lett. 1227, 4 (1987). https://doi.org/10.1209/0295-5075/4/11/003
  * (22) Roy, S., Sinha, S., Hansen, A.: Role of pore-size distribution on effective rheology of two-phase flow in porous media. Front. Water 3, 709833 (2021). https://doi.org/10.3389/frwa.2021.709833
  * (23) Fyhn, H., Sinha, S., Roy, S., Hansen, A.: Rheology of immiscible two-phase flow in mixed wet porous media: Dynamic pore network model and capillary fiber bundle model results. Transp. Porous Med. 139, 491 (2021). https://doi.org/%****␣compress.bbl␣Line␣400␣****10.1007/s11242-021-01674-3
  * (24) Leverett, M.C.: Capillary behavior in porous solids. Trans. AIME 142, 152 (1941). https://doi.org/10.2118/941152-G
  * (25) Li, X., Yortsos, Y.C.: Bubble growth and stability in an effective porous medium. Phys. Fluids 6, 1663 (1994). https://doi.org/10.1063/1.868229
  * (26) Reynolds, C.A., Krevor, S.: Characterizing flow behavior for gas injection: Relative permeability of co2-brine and n2-water in heterogeneous rocks. Water Resources Res. 51, 9464 (2015). https://doi.org/10.1002/2015WR018046
  * (27) Abidoye, L.K., Khudaida, K.J., Das, D.B.: Geological carbon sequestration in the context of two-phase flow in porous media: A review. Crit. Rev. Env. Sc. Tech. 45, 1105 (2015). https://doi.org/10.1080/10643389.2014.924184
  * (28) Iglauer, S., Paluszny, A., Rahman, T., Zhang, Y., Wülling, W., Lebedev, M.: Residual trapping of co2 in an oil-filled, oil-wet sandstone core: Results of three-phase pore-scale imaging. Geophys. Res. Lett. 46, 11146 (2019). https://doi.org/10.1029/2019GL083401
  * (29) Niblett, D., Mularczyk, A., Niasar, V., Eller, J., Holmes, S.: Two-phase flow dynamics in a gas diffusion layer - gas channel - microporous layer system. J. Power Sources 471, 228427 (2020). https://doi.org/10.1016/j.jpowsour.2020.228427
  * (30) Huang, G., Zhu, Y., Liao, Z., Ouyang, X.L., Jiang, P.X.: Experimental investigation of transpiration cooling with phase change for sintered porous plates. Int. J. Heat Mass Transfer 114, 1201 (2017). https://doi.org/10.1016/j.ijheatmasstransfer.2017.05.114
  * (31) Gedupudi, S., Zu, Y.Q., Karayiannis, T.G., Kenning, D.B.R., Yan, Y.Y.: Confined bubble growth during flow boiling in a mini/micro-channel of rectangular cross-section part i: Experiments and 1-d modelling. Inr. J. Therm. Sci. 50, 250 (2011). https://doi.org/10.1016/j.ijthermalsci.2010.09.001
  * (32) Li, D., Wu, G.S., Wang, W., Wang, Y.D., Liu, D., Zhang, D.C., Chen, Y.F., Peterson, G.P., Yang, R.: Enhancing flow boiling heat transfer in microchannels for thermal management with monolithically-integrated silicon nanowires. Nano Lett. 12, 3385 (2012). https://doi.org/10.1021/nl300049f
  * (33) Li, W., Wang, Z., Yang, F., Alam, T., Jiang, M., Qu, X., Kong, F., Khan, A.S., Liu, M., Alwazzan, M., Tong, Y., Li, C.: Supercapillary architecture-activated two-phase boundary layer structures for highly stable and efficient flow boiling heat transfer. Adv. Matter 32, 1905117 (2020). https://doi.org/10.1002/adma.201905117
  * (34) Bremer, J., Sundmacher, K.: Operation range extension via hot-spot control for catalytic co2 methanation reactors. React. Chem. Eng. 4, 1019 (2019). https://doi.org/10.1039/C9RE00147F
  * (35) Sapin, P., Gourbil, A., Duru, P., Fichot, F., Prat, M., Quintard, M.: Reflooding with internal boiling of a heating model porous medium with mm-scale pores. Int. J. Heat Mass Transfer 99, 512 (2016). https://doi.org/10.1016/j.ijheatmasstransfer.2016.04.013
  * (36) Sun, Y., Zhang, L., Xu, H., Zhong, X.: Subcooled flow boiling heat transfer from microporous surfaces in a small channel. Int. J. Therm. Sci. 50, 881 (2011). https://doi.org/10.1016/j.ijthermalsci.2011.01.019
  * (37) C., R., Nimmo, J.R.: Modeling of soil water retention from saturation to oven dryness. Water Resour. Res. 30, 701 (1994). https://doi.org/10.1029/93WR03238
  * (38) Sinha, S., Hansen, A., Bedeaux, D., Kjelstrup, S.: Effective rheology of bubbles moving in a capillary tube. Phys. Rev. E 87, 025001 (2013). https://doi.org/10.1103/PhysRevE.87.025001
  * (39) Blunt, M.J.: Flow in porous media pore-network models and multiphase flow. Current Opinion Colloid & Interface Science 6, 197 (2001). https://doi.org/10.1016/S1359-0294(01)00084-X
  * (40) Meakin, P., Tartakovsky, A.M.: Modeling and simulation of pore-scale multiphase fluid flow and reactive transport in fractured and porous media. Reviews of Geophysics 47, 3002 (2009). https://doi.org/10.1029/2008RG000263
  * (41) Joekar-Niasar, V., M., H.S.: Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: A review. Crit. Rev. Environ. Sci. Technol. 42, 1895 (2012). https://doi.org/10.1080/10643389.2011.574101
  * (42) Gunstensen, A.K., H., R.D., S., Z., G., Z.: Lattice boltzmann model of immiscible fluids. Phys. Rev. A 43, 4320 (1991). https://doi.org/10.1103/PhysRevA.43.4320
  * (43) Ramstad, N. T. Idowu, Nardi, C., Øren, P.E.: Relative permeability calculations from two-phase flow simulations directly on digital images of porous rocks. Transp. Porous Media 94, 487 (2012). https://doi.org/10.1007/s11242-011-9877-8
  * (44) Qiu, R.F., You, Y.C., Zhu, C.X., Chen, R.Q.: Lattice boltzmann simulation for high-speed compressible viscous flows with a boundary layer. Appl. Math. Model. 48, 567 (2017). https://doi.org/10.1016/j.apm.2017.03.016
  * (45) Guo, S., Feng, Y., Jacob, J., Renard, F., Sagaut, P.: An efficient lattice boltzmann method for compressible aerodynamics on d3q19 lattice. J. Comput. Phys. 418, 109570 (2020). https://doi.org/10.1016/j.jcp.2020.109570
  * (46) Washburn, E.W.: The dynamics of capillary flow. Phys. Rev. 17, 273 (1921). https://doi.org/10.1103/PhysRev.17.273
  * (47) Vazquez, A., Leifer, I., Sánchez, R.M.: Consideration of the dynamic forces during bubble growth in a capillary tube. Chem. Eng. Sc. 65, 4046 (2010). https://doi.org/10.1016/j.ces.2010.03.041
  * (48) Welch, S.W.J.: Direct simulation of vapor bubble growth. Int. J. Heat Mass Trans. 41, 1655 (1998). https://doi.org/10.1016/S0017-9310(97)00285-8
  * (49) Kenning, D.B.R., Wen, D.S., Das, K.S., Wilson, S.K.: Confined growth of a vapour bubble in a capillary tube at initially uniform superheat: experiments and modelling. Int. J. Heat Mass Trans. 49, 4653 (2006). https://doi.org/10.1016/j.ijheatmasstransfer.2006.04.010
  * (50) Roy, S., Hansen, A., Sinha, S.: Effective rheology of two-phase flow in a capillary fiber bundle model. Front. Phys. 7, 92 (2019). https://doi.org/10.3389/fphy.2019.00092


# Demonstration of measurement-only blind quantum computing

Chiara Greganti1, Marie-Christine Roehsner1, Stefanie Barz1,2, Tomoyuki
Morimae3, Philip Walther1 1 University of Vienna, Faculty of Physics, Austria,
2 Present address: University of Oxford, Clarendon Laboratory, UK, 3 ASRLD
Unit, Gunma University, 1-5-1 Tenjin-cho Kiryu-shi Gunma-ken, 376-0052, Japan

###### Abstract

Blind quantum computing allows for secure cloud networks of quasi-classical
clients and a fully-fledged quantum server. Recently, a new protocol has been
proposed, which requires a client to perform only measurements. We demonstrate
a proof-of-principle implementation of this measurement-only blind quantum
computing, exploiting a photonic set-up to generate four-qubit cluster states
for computation and verification. Feasible technological requirements for the
client and the device-independent blindness make this scheme very applicable
for future secure quantum networks.

Introduction. – Quantum physics enables one to enhance security for processing
data over a distributed network. In particular, quantum cloud computing allows
quasi-classical clients (i.e. clients with a limited amount of quantum
resources, such as qubit preparation or detection) to  do calculations beyond
their computational power, namely perform quantum algorithms. The first
proposed and demonstrated two-party secure quantum cloud computation protocol
is known as blind quantum computing (BQC) Broadbent et al. (2009): a client,
Alice, who can generate only single-qubit states, delegates her quantum
computing to a remote server, Bob, who has a fully-fledged quantum computer,
without leaking any of her privacy. Many theoretical studies based on this
have been performed recently Morimae and Fujii (2013); Fitzsimons and Kashefi
(2012); Dunjko et al. (2012); Morimae et al. (2015); Morimae and Fujii (2012);
Morimae (2012); Giovannetti et al. (2013); Mantri et al. (2013); Sueki et al.
(2013); Li et al. (2014), and also experimental demonstrations have been
reported Barz et al. (2012, 2013); Fisher et al. (2014). A simplified and
novel version for secure quantum computing consists of a two-party protocol
Morimae and Fujii (2013); Morimae (2014), where Alice only makes measurements
and Bob’s blindness is proven by the no-signaling principle Popescu and
Rohrlich (1994). Here blindness indicates that whatever Bob does, he cannot
learn any of Alice’s privacy.

In order to underline the feasibility of the measurement-only BQC we
demonstrate the computation protocol in a photonic experiment. Bob generates
four-qubit resource states that are used by Alice to implement generic two-
qubit entangling gates and verification protocols.

Theory. – The idea of measurement-only BQC is shown in Fig.1. Bob generates a
resource state for measurement-based quantum computing (MBQC) Raussendorf and
Briegel (2001) Briegel et al. (2009), and sends the corresponding qubits, one
by one, through a one-way quantum channel to Alice. She measures each qubit
according to her program. For any kind of a malicious Bob, he cannot learn
anything about Alice’s quantum computation, because information is sent only
in one direction. The no-signaling principle then ensures that if Alice and
Bob share a system (classical, quantum, or superquantum) and she measures her
part, this does not transmit any information to Bob. This principle is more
fundamental than quantum physics Popescu and Rohrlich (1994) and consequently
provides security even against superquantum attacks.

Figure 1: Measurement-only blind quantum computing. Alice, a client with a
classical device (laptop) and a quantum detector device (grey box), receives
quantum information via a one-way channel from Bob, a quantum server capable
of generating entangled quantum states as universal resources. Here Bob’s
resource corresponds to a fault-tolerant 3D-cluster state, but it can be any
other measurement-based quantum computation resource.

Figure 2:  Schematic Alice’s (left) and Bob’s (right) optical set-ups. Right:
The four-qubit entangled states are generated via a fs-pulsed laser pumping a
type-II BBO crystal in a double-passage scheme. The compensator, HWPs and PBSs
allow Bob to select the desired four-qubit resource. Left: The detection
requires HWP, QWP and PBS for measuring in different polarization basis, and
single-photon APD per qubit. The laboratories are connected by single-mode
fibers as quantum channels. See main text for acronyms.

We remark that the original BQC examines a different approach. A quasi-
classical Alice must be able to: generate randomly-rotated single qubit
states, send these via a quantum channel to Bob, and interact via classical
channels in order to control and receive the computation results. Recently,
the single-qubit generation requirement is extended to the coherent state
generation Dunjko et al. (2012) or single-qubit measurement Hajdusek et al.
(2015).  The measurement-only BQC scheme centers on a quasi-classical Alice,
who now only receives qubits for measuring those. This alternative concept
produces a practical computing protocol. In particular in optical systems the
technological demand for a client is readily available. An additional feature
of measurement-only BQC is the device-independent blindness: even if Alice
owns a malicious device, probably bought from another company, no information
is transmitted to Bob because of the no-signaling principle.  In our case, we
obtain device independence regarding blindness; recently, the concept of
device independence has also been generalized to verifiability Hajdusek et al.
(2015); Gheorghiu et al. (2015).  Furthermore the full-fledged quantum
computer can be based on  any model of MBQC (see Ref.Morimae and Fujii (2013)
for more detailed discussion).

The concept of secure quantum computing opened-up feasible verification
methods Fitzsimons and Kashefi (2012); Barz et al. (2013); Morimae (2014)
where Alice can test whether Bob is performing the computation correctly. It
was shown that the verification is possible for the original protocol
Fitzsimons and Kashefi (2012); Barz et al. (2013) as well as for the
measurement-only protocol Morimae (2014). The central idea in these protocols
is that Alice secretly hides some ”trap” qubits in the resource state. This
fundamentally reduces to the situation where Alice tries to verify, Bob’s
quantum resource with a minimal set of measurements Tóth and Gühne (2005);
Greganti et al. (2015); Knips et al. (2014). If Bob deviates from the correct
protocol, he changes the states of the traps, and if Alice detects the change
of any trap, she can detect Bob’s malicious behaviour and abort the
computation. The security corresponds to the probability that Alice does not
accept the results received by a cheating Bob.

In measurement-only BQC the trap qubits are randomly prepared and placed via
measurements by Alice within the computation resource and are associated to
qubit states in $Z$ and $X$ basis (corresponding to Pauli operators
$\sigma_{z}$ and $\sigma_{x}$ respectively). When Alice receives the traps,
she measures the qubits in the respective basis.

For the special case of a four-qubit linear cluster state a verification
protocol with only two different trap measurements exists. This protocol runs
as follows:

  * 1.

Bob prepares the four-qubit linear cluster state, and sends each qubit one by
one to Alice.

  * 2.

Alice chooses one of the two tests randomly below:

    * (a)

Alice measures qubits 1 and 3 in the $Z$ basis, and qubits 2 and 4 in the $X$
basis.

    * (b)

Alice measures qubits 2 and 4 in the $Z$ basis, and qubits 1 and 3 in the $X$
basis.

If she chooses the option (a), qubits 2 and 4 becomes trap qubits. If any trap
qubit is changed, i.e. she gets not the expected result, then she detects
Bob’s malicious behavior. We call this option (1,3) test. On the other hand,
if she chooses the option (b), qubits 1 and 3 become trap qubits, and she can
check those. We call this option (2,4) test. We show now that Bob has to
prepare the exact four-qubit linear cluster state in order to pass all Alice’s
trap tests in the limit of $n$ repetitions where $n$ tends toward infinity.
In the original verification protocols Morimae (2014); Fitzsimons and Kashefi
(2012), it is shown that the probability that Alice is fooled by Bob can be
exponentially small, by using quantum error correcting codes.

Here, let us show the case without quantum error correcting code, which leads
to a probability of accepting a wrong outcome to be polynomially small. We
want to point out that this probability can be minimized to become
exponentially small by exploiting standard error amplification techniques
Kitaev et al. (2002) via repeating the computation a number of times
proportional to $n$. Bob can generate any state, but in order to pass the
(1,3) test, Bob has to prepare the state

$\begin{split}|\Psi\rangle\equiv\frac{1}{2}\Big{(}|0+0+\rangle|a_{1}\rangle+e^{i\theta_{2}}|0-1-\rangle|a_{2}\rangle+\\\
+e^{i\theta_{3}}|1-0+\rangle|a_{3}\rangle+e^{i\theta
4}|1+1-\rangle|a_{4}\rangle\Big{)}\end{split}$ (1)

where $\\{|a_{j}\rangle\\}$ are certain states of Bob’s ancilla, which are
normalized $|\langle a_{j}|a_{j}\rangle|^{2}=1$, but not necesarily mutually
orthogonal. We consider ancilla states, since Bob could prepare a larger
system and keep a subsystem. Since Bob does not know which option Alice takes,
this state also has to pass the other Alice’s test, i.e., (2,4) test. Assume
that Alice gets the result $(0,0)$ when she measures qubits 2 and 4. Then the
state after the measurement becomes

$\begin{split}|\Psi^{\prime}\rangle\equiv\frac{1}{2}\Big{(}|0\rangle_{1}|0\rangle_{3}|a_{1}\rangle+e^{i\theta_{2}}|0\rangle_{1}|1\rangle_{3}|a_{2}\rangle+\\\
+e^{i\theta_{3}}|1\rangle_{1}|0\rangle_{3}|a_{3}\rangle+e^{i\theta
4}|1\rangle_{1}|1\rangle_{3}|a_{4}\rangle\Big{)}.\end{split}$ (2)

In order to pass the (2,4) test (again in the limit of $n$ repetitions), this
state must be $|+\rangle_{1}|+\rangle_{3}|b\rangle$ for a certain state
$|b\rangle$. This means that, first, the reduced density operator of
$|\Psi^{\prime}\rangle$ for Bob’s ancilla

$\displaystyle\rho=\frac{1}{4}\sum_{j=1}^{4}|a_{j}\rangle\langle a_{j}|$

must have rank 1, which leads to $|a_{j}\rangle=|a_{k}\rangle$ (up to a phase
factor). Now the state is

$\displaystyle\frac{1}{2}\Big{(}|0\rangle_{1}|0\rangle_{3}+e^{i\theta_{2}^{\prime}}|0\rangle_{1}|1\rangle_{3}+e^{i\theta_{3}^{\prime}}|1\rangle_{1}|0\rangle_{3}+e^{i\theta
4^{\prime}}|1\rangle_{1}|1\rangle_{3}\Big{)}|a_{1}\rangle.$

Next, in order for this state to be $|+\rangle_{1}|+\rangle_{3}|b\rangle$,
$\theta_{j}^{\prime}=0$ for all $j=2,3,4$.
Therefore, repeating both tests Alice verifies that Bob has the exact four-
qubit cluster state, except for a small probability of undetected cheating.

In a general case, Alice can choose to use the resource state for either
verification or computation. Increasing the number of verifications per
computation provides a higher level of security at the cost of efficiency. The
probability of undetected errors, i. e. Bob cheats in the computation and not
in the verification, is linearly bounded as in Barz et al. (2013). In Ref.
Barz et al. (2013) Alice is assumed to send qubits, whereas here Alice
measures the received qubits. Remarkably the same technique of verification
can be applied in both schemes: Alice generates trap qubits via either
choosing Bob’s measurement settings or by directly performing the measurements
at her side. This allows for the same analysis, discussed already in Ref. Barz
et al. (2013). In our case, Alice can randomly choose between the two
verification options and a regular computation on a four-qubit linear cluster.
Ref. Hayashi and Morimae (2015) describes the asymptotic behavior of the
scaling for linear cluster states of increasing length.

It is worth to note that in Ref. Broadbent et al. (2009), a random-number
generator is necessary for the blindness, whereas in measurement-only BQC
Morimae and Fujii (2013), no random-number generator is required for Alice to
guarantee the blindness. If we add the option of the verification, both
protocols require random-number generators, since Alice has to randomly place
trap qubits. Nevertheless the use of quantum random numbers is nowadays
accessible at the consumer grade Sanguinetti et al. (2014).

Experiment and results. – We practically realize a proof-of-principle
implementation of the protocol using photons, computing two-qubit entangling
gates and verifying two single trap qubits. In contrast to the proposed
theoretical scheme Morimae (2014), where traps are hidden within the
computation resource, our experiment exploits a four-qubit cluster state
either for a computation or a verification run, due to the number of available
qubits.

The four-qubit resource for measurement-only BQC, (Fig. 3 a, b), is produced
in Bob’s laboratory via a photonic set-up in a so-called ”railway-crossing
configuration” (see Fig. 2). A double spontaneous parametric down conversion
process allows to generate two pairs of polarization entangled photons.
Interferometers with polarizing beam splitters (PBSs) entangle the four
photons. Additional half-wave plates (HWPs) on both pairs directions enable
the generation of different graph states. The scheme has been already
exploited in several other works to create four-qubit linear cluster states
and states that can be obtained from them via local complementations (see e.g.
Barz et al. (2012, 2014)).  Here we focus on the generation of a four-qubit
star cluster state and a four-qubit linear cluster state, respectively
(Fig.3.a and b):

$\displaystyle|C_{star}\rangle=\frac{1}{\sqrt{2}}(|++0+\rangle+|--1-\rangle)_{1234},$
(3)
$\displaystyle|C_{lin}\rangle=\frac{1}{2}(|0+0+\rangle+|0-1-\rangle+|1-0+\rangle+|1+1-\rangle)_{1234},$
(4)

where $|\pm\rangle=(|0\rangle\pm|1\rangle)/\sqrt{2}$ are the eigenstates of
the Pauli operator $\sigma_{x}=X$. Remarkably these cluster states belong to
different classes of entanglement.

The four-qubit star cluster was generated within this setup only recently (see
Greganti et al. (2015) for details) and now exploited for quantum information
computing. Switching between the two entangled classes involves: preparing
specific Bell states at each SPDC process, different photonic interferences
between the two pairs of photons, precise wavelength-scale alignment, and,
therefore, high stability.

Alice’s laboratory consists of four HWPs, four quarter-wave plates (QWPs),
four PBSs and eight single photon counting detectors (APDs) in order to speed
up the data acquisition (four APDs are sufficient to measure all possible
polarization-basis of four-qubit state) and proceed to a complete analysis
using quantum state tomography (QST) James et al. (2001). The connecting
quantum channel from Bob to Alice is achieved by four single mode fibers,
which carry the photonic qubits. We reconstructed through over-complete QST
the density matrix of the two four-qubit resources, obtaining fidelities of
the state with respect to the ideal star cluster and linear cluster of
$F=0.731\pm 0.008$ and $F=0.676\pm 0.007$ (under local unitary operations),
respectively (see SI for the density matrix histograms).

Figure 3:  Bob’s quantum resources: four-qubit (a) star and (b) linear cluster
states. (c) Implementation of an entangling gate on a four-qubit star cluster.
$C$ ($C^{\prime}$) and $T$ ($T^{\prime}$) are the control and target input
(output) logical qubits. (d) Verification (1,3)-test on a linear cluster
state. The dashed frame corresponds to a single-qubit measurement, whereas the
red frame represents a trap. Verification (2,4)-test is equivalent under
exchange of trap and measured qubits.

Computation – The four-qubit star and linear cluster states are the minimal
resources for one-way computation, since the full universal set of gates can
be reproduced Raussendorf and Briegel (2001). This has been already
demonstrated in few works Walther et al. (2005); Prevedel et al. (2007); Chen
et al. (2007); Vallone et al. (2008); Tokunaga et al. (2008); Bell et al.
(2013). In this work we reproduce different two-qubit entangling gates using
the star cluster, in order to validate the computation from Alice’s device. An
entangling gate is performed using a star cluster, where qubits 1 and 2 are
acting as input control ($C$) and target ($T$) qubits, respectively, whereas
qubits 1 and 4 present the output control ($C^{\prime}$) and target
($T^{\prime}$) qubits, respectively, as shown in Fig.3c. Different
combinations of measurement bases for qubit 2 and 3 enable to create
entanglement between the output qubits. Detailed analysis for some entangling
gates are reported in the SI. In table 1 we present, as example, the results
related to measuring qubit 2 and 3 in the $Y_{2}X_{3}$ basis (where $Y$ is the
Pauli $\sigma_{y}$), which corresponds to implementing a Controlled NOT (CNOT)
gate for input state $|++_{i}\rangle$ and $|+-_{i}\rangle$ (where
$|+_{i}\rangle=(|0\rangle+i|1\rangle)/\sqrt{2}$ and
$|-_{i}\rangle=(|0\rangle-i|1\rangle)/\sqrt{2}$ are the eigenstates of $Y$) up
to local unitary operations. The two-qubit output states are analysed through
two-qubit QST with acquisition-time of $600s$ per measurement setting. The
corresponding uncertanties are due to Poissonian counting statistics and
represent only a lower bound for the errors.

$\mathbf{s_{2}}\mathbf{s_{3}}$ |  Ideal Output State |  Fidelity
---|---|---
$00$ |  $(|0+_{i}\rangle+i|1-_{i}\rangle)_{14}/\sqrt{2}$ |  F=0.87 $\pm$ 0.03
$01$ |  $(|0+_{i}\rangle+i|1-_{i}\rangle)_{14}/\sqrt{2}$ |  F=0.74 $\pm$ 0.04
$10$ |  $(|0-_{i}\rangle-i|1+_{i}\rangle)_{14}/\sqrt{2}$ |  F=0.77 $\pm$ 0.03
$11$ |  $(|0-_{i}\rangle-i|1+_{i}\rangle)_{14}/\sqrt{2}$ |  F=0.77 $\pm$ 0.04
Table 1: Results from measuring qubit 2 and qubit 3 of the star cluster onto
$Y_{2}X_{3}$, which corresponds to a CNOT gate on states $|++_{i}\rangle$ and
$|+-_{i}\rangle$ up to a $(Z_{1}Z_{4})^{s_{3}+1}$, where $s_{2}$ and $s_{3}$
are the measurement outcomes. The fidelities of the tomographic reconstructed
two-qubit state with respect of the ideal state are reported.

Verification – The four-qubit linear cluster allows verifying computation with
only two different trap measurements and is especially suited for the
verification protocol as described in the theory section. We present the
results for this state along with the results for a four-qubit star cluster
state. As has been shown in Barz et al. (2013), the probability that Alice is
fooled by Bob is bounded in such a setting.

For the case of the four-qubit linear cluster state we implement the (1,3) and
(2,4) tests by having Alice choosing the respective basis. Per trap we analyze
the measurement outcomes in order to quantify the probabilities that Alice
obtains the correct state, see Fig.4. For single trap the results are within
the values $[0.74\pm 0.03,0.98\pm 0.01]$, where the range is due to unbalanced
phase noises in the set-up. Each of Alice’s measures has an acquisition-time
of $1h$ to decrease the uncertainty. Alice verified the resource with non-
ideal probability, due to experimental imperfections of the set-up, which are
present during the generation of the four-qubit resource as already seen from
the full QST fidelity of the state.

Additionally, we performed the verification protocol on a four-qubit star
cluster state, which we used before to implement entangling gates. We report
the two trap tests performed on the star cluster state, equivalently to the
linear cluster case, see Fig.5. In this case in order to get two trap qubits
each time Alice measures $Z_{1}Z_{3}$ (expecting trap qubits in $X$ basis) for
(1,3)-test, and $Z_{2}X_{4}$ (expecting one trap qubit in $X$ basis and one in
$Z$) for (2,4)-test. The single probabilities of individual trap qubits,
corresponding to Alice’s expected results, and according to respective
measurement outcomes of the non-trap qubits are in the range $[0.90\pm
0.04,1.00-0.16]$ with an acquisition time of $600s$ per single measurement
(see SI for details). The imbalance of the obtained probabilities with respect
to the quantum state fidelity is due to asymmetric noise.

Here we want to point out that the small increased value of fidelity of the
star cluster, with respect to the linear cluster, leads to significant
improvements for the verification results.

Discussion. – The use of four-qubit photonic cluster states allows us to prove
the feasibility of two-party measurement-only BQC in current quantum optics
laboratories. Nevertheless the demonstration can be expanded to several
quantum systems and other MBQC models. In the photonic case, we want to
emphasize that just one HWP, one QWP, one PBS, and one APD would be sufficient
for Alice to measure every qubit received from Bob and consequently to
implement computation and verification. The only additional requirement in
Bob’s laboratory would be a time-delay multiplexer (such as the one used in
Collins et al. (2013)) or a delay line in each photon’s path (such as in
Megidish et al. (2013)). The tomographic reconstruction of Bob’s resource
state, as we did in our experiment, is in fact superfluous for Alice’s
computation, since already from the single-qubit measurement she can verify
Bob’s state. The quantum power required for Alice is then restricted to
measuring the state of the qubits. It is important to note the high losses,
either due to low detection efficiencies or imperfect quantum channels, would
break Alice’s computation. However, the threshold for losses can be increased
by using fault tolerant MBQC models, which are robust against errors and
losses Barrett and Stace (2010), and besides, detection devices with almost
unit efficiencies are now available Lita et al. (2008); Marsili et al. (2013).
In conclusion the demonstrated protocol constitutes a step further to more
realistic secure quantum computing models.

Figure 4: Probability that Alice receives the expected outcomes for the
(1,3)-test and (2,4)-test on a linear cluster state. According to the
measurement outcomes of the non-trap qubits (shown on the abscissa) we report
the probability that measurements on each of Alice’s trap qubits return the
expected result (i.e. Alice trusts the state Bob sent).  Figure 5:
Probability that Alice receives the expected outcomes for the two verification
tests on a star cluster. According to the measurement outcomes of the non-trap
qubits (shown on the abscissa) we report the probability that measurements on
each of Alice’s trap qubits return an expected result (i.e. Alice trusts the
state Bob sent).

###### Acknowledgements.

This work was supported by the European Commission, Initial Training Network
PICQUE (No. 608062), QUILMI (No.295293), EQUAM (No. 323714), GRASP (No.
613024), and the Vienna Center for Quantum Science and Technology (VCQ), the
Austrian Science Fund (FWF) through START (No. Y585-N20) and the doctoral
programme CoQuS, the Vienna Science and Technology Fund (WWTF) under grant
ICT12-041, and the Air Force Office of Scientific Research, Air Force Material
Command, United States Air Force, under grant number FA8655-11-1-3004, and the
Tenure Track System by MEXT Japan and KAKENHI 26730003 and 15H00850 by JSPS.

## References

  * Broadbent et al. (2009) A. Broadbent, J. Fitzsimons, and E. Kashefi, in _Proceedings of the 50th Annual Symposium on Foundations of Computer Science_ (2009), pp. 517–526.
  * Morimae and Fujii (2013) T. Morimae and K. Fujii, Phys. Rev. A 87, 050301 (2013).
  * Fitzsimons and Kashefi (2012) J. Fitzsimons and E. Kashefi, arXiv:1203.5217 (2012).
  * Dunjko et al. (2012) V. Dunjko, E. Kashefi, and A. Leverrier, Phys. Rev. Lett. 108, 200502 (2012).
  * Morimae et al. (2015) T. Morimae, V. Dunjko, and E. Kashefi, Quantum Information and Computation 15, 0200 (2015).
  * Morimae and Fujii (2012) T. Morimae and K. Fujii, Nat. Commun. 3, 1036 (2012).
  * Morimae (2012) T. Morimae, Phys. Rev. Lett. 109, 230502 (2012).
  * Giovannetti et al. (2013) V. Giovannetti, L. Maccone, T. Morimae, and T. Rudolph, Phys. Rev. Lett. 111, 230501 (2013).
  * Mantri et al. (2013) A. Mantri, C. Pérez-Delgado, and J. Fitzsimons, Phys. Rev. Lett. 111, 230502 (2013).
  * Sueki et al. (2013) T. Sueki, T. Koshiba, and T. Morimae, Phys. Rev. A 87, 060301 (2013).
  * Li et al. (2014) Q. Li, W. H. Chan, C. Wu, and Z. Wen, Phys. Rev. A 89, 040302 (2014).
  * Barz et al. (2012) S. Barz, E. Kashefi, A. Broadbent, J. Fitzsimons, A. Zeilinger, and P. Walther, Science 335, 303 (2012).
  * Barz et al. (2013) S. Barz, J. F. Fitzsimons, E. Kashefi, and P. Walther, Nature Physics 9, 727 (2013).
  * Fisher et al. (2014) K. Fisher, A. Broadbent, L. Shalm, Z. Yan, J. Lavoie, R. Prevedel, T. Jennewein, and K. Resch, Nature Communications 5, 3074 (2014).
  * Morimae (2014) T. Morimae, Phys. Rev. A 89, 060302 (2014).
  * Popescu and Rohrlich (1994) S. Popescu and D. Rohrlich, Foundations of Physics 24, 379 (1994), ISSN 0015-9018.
  * Raussendorf and Briegel (2001) R. Raussendorf and H. Briegel, Phys. Rev. Lett. 86, 5188 (2001).
  * Briegel et al. (2009) H. Briegel, D. E. Browne, W. Dür, R. Raussendorf, and M. Van den Nest, Nature Phys. 5, 19 (2009), ISSN 1745-2473.
  * Hajdusek et al. (2015) M. Hajdusek, C. Pérez-Delgado, and J. Fitzsimons, arXiv:1502.02563v1 (2015).
  * Gheorghiu et al. (2015) A. Gheorghiu, E. Kashefi, and P. Wallden, arXiv:1502.02571v2 (2015).
  * Tóth and Gühne (2005) G. Tóth and O. Gühne, Phys. Rev. A 72, 022340 (2005).
  * Greganti et al. (2015) C. Greganti, M. C. Roehsner, S. Barz, M. Waegell, and P. Walther, Phys. Rev. A, 91 022325 (2015).
  * Knips et al. (2014) L. Knips, C. Schwemmer, N. Klein, M. Wieśniak, and H. Weinfurter, arXiv:1412.5881 (2014).
  * Kitaev et al. (2002) A. Y. Kitaev, A. H. Shen, and M. N. Vyalyi, _Classical and Quantum Computation_ (American Mathematical Society, Boston, MA, USA, 2002), ISBN 0821832298\.
  * Hayashi and Morimae (2015) M. Hayashi and T. Morimae, arXiv:1505.07535 (2015).
  * Sanguinetti et al. (2014) B. Sanguinetti, A. Martin, H. Zbinden, and N. Gisin, Phys. Rev. X 4, 031056 (2014).
  * Barz et al. (2014) S. Barz, R. Vasconcelos, C. Greganti, M. Zwerger, W. Duer, H. Briegel, and P. Walther, arXiv. 1308.5209 (2014).
  * James et al. (2001) D. James, P. Kwiat, W. Munro, and A. White, Phys. Rev. A 64, 52312 (2001).
  * Walther et al. (2005) P. Walther, K. Resch, T. Rudolph, E. Schenck, H. Weinfurter, V. Vedral, M. Aspelmeyer, and A. Zeilinger, Nature 434, 169 (2005).
  * Prevedel et al. (2007) R. Prevedel, P. Walther, F. Tiefenbacher, P. Böhi, R. Kaltenbaek, T. Jennewein, and A. Zeilinger, Nature 445, 65 (2007).
  * Chen et al. (2007) K. Chen, C.-M. Li, Q. Zhang, Y.-A. Chen, A. Goebel, S. Chen, A. Mair, and J.-W. Pan, Phys. Rev. Lett. 99, 120503 (pages 4) (2007).
  * Vallone et al. (2008) G. Vallone, E. Pomarico, F. D. Martini, and P. Mataloni, Phys. Rev. A 78, 042335 (2008).
  * Tokunaga et al. (2008) Y. Tokunaga, S. Kuwashiro, T. Yamamoto, M. Koashi, and N. Imoto, Phys. Rev. Lett. 100, 210501 (2008), ISSN 1079-7114.
  * Bell et al. (2013) B. Bell, M. Tame, A. Clark, R. Nock, W. Wadsworth, and J. Rarity, New J. Phys. 15, 053030 (2013).
  * Collins et al. (2013) M. Collins, C. Xiong, I. Rey, T. Vo, J. He, S. Shahnia, T. K. C. Reardon, M. Steel, A. Clark, and B. Eggleton, Nat. Comm. 4, 2582 (2013).
  * Megidish et al. (2013) E. Megidish, A. Halevy, T. Shacham, T. Dvir, L. Dovrat, and H. S. Eisenberg, Phys. Rev. Lett. 110, 210403 (2013).
  * Barrett and Stace (2010) S. Barrett and T. M. Stace, Phys. Rev. Lett. 105, 200502 (2010).
  * Lita et al. (2008) A. E. Lita, A. J. Miller, and S. W. Nam, Opt. Express 16, 3032 (2008).
  * Marsili et al. (2013) F. Marsili, V. B. Verma, J. A. Stern, S. Harrington, T. G. A. E. Lita, I. Vayshenker, B. Baek, M. D. Shaw, R. P. Mirin, and S. W. Nam, Nat. Phot. 7, 210 (2013).