# Rheology of dense fiber suspensions: Origin of yield stress, shear thinning
and normal stress differences

Monsurul Khan  Department of Mechanical Engineering, Purdue University, IN
47905, USA  Rishabh V. More  Department of Mechanical Engineering, Purdue
University, IN 47905, USA  Luca Brandt  Linn´e Flow Centre and SeRC (Swedish
e-Science Research Centre), KTH Mechanics, SE 100 44 Stockholm, Sweden  Arezoo
M. Ardekani  Department of Mechanical Engineering, Purdue University, IN
47905, USA

(August 27, 2024)

###### Abstract

We explain the origins of yield stress, shear-thinning, and normal stress
differences in rigid fiber suspensions. We investigate the interplay between
the hydrodynamic, colloidal attractive and repulsive, and inter-fiber contact
interactions. The shear-thinning viscosity and finite yield stress obtained
from the computational model are in quantitative agreement with experiential
results from the literature. In this study, we show that attractive
interactions result in yield stress and shear-thinning rheology in the
suspensions of rigid fibers. This is an important finding, given the ongoing
discussion regarding the origin of the yield stress for suspensions of fibers.
The ability of the proposed model to quantitatively predict the rheology is
not limited to only shear thinning and yield stress but also extends to normal
stresses.

Suggested keywords

††preprint: APS/123-QED

##  I Introduction

Fiber suspensions are widely encountered in natural and industrial
applications, with examples in paper and pulp production, biomass solutions,
and chemical processing [1, 2, 3, 4, 5]. Under shear, these suspensions
display several non-Newtonian properties such as the Weissenberg effect [6,
7], shear thinning [8, 9, 10], non-zero normal stresses [11, 12], yield stress
[10]. Specifically, shear thinning in fiber suspensions has been an active
area of research, and consequently, the literature provides many
phenomenological explanations, such as the increase in the effective particle
size due to the presence of electric double layer [13], excluded volume
interactions between rigid fibers [14], elastic bending of flexible micron-
sized or nano-sized fibers [15, 16], and fiber aggregation [17] and nonlinear
lubrication force [18]. However, numerical models incorporating these
phenomenological explanations are missing from the literature.

Non-Newtonian rheology arises from a variety of interactions such as
hydrodynamics, cohesive [19, 20], contact [21, 22], and is influenced by the
physical properties of fibers, e.g., roughness, shape, size distribution, etc.
[23, 24] Each of these interactions leads to a corresponding stress scale as
${F}_{A}/\pi d^{2}$ (d is the fiber characteristic scale) that is competing
with hydrodynamic interactions that scale as $\pi\eta\dot{\gamma}d^{2}$, where
$\dot{\gamma}$ is the imposed shear rate and $\eta$ is the viscosity of the
suspending fluid. The competition between these stress scales could lead to a
rate-dependent rheological behavior [25, 26]. Among these different
interactions between fibers, the influence of short-range hydrodynamic forces
and direct mechanical contact has been understood through theoretical modeling
[27, 28] and numerical simulations [21, 29]. However, there is little
documentation on the study of fiber suspensions with colloidal interactions,
even though van der Waals forces [30], depletion forces due to dissolved non-
interacting polymers, [31], presence of external fields [32] can lead to
attractive forces between fibers. Therefore, it is crucial to understand the
effect of attractive and repulsive forces in fiber suspensions as they produce
nonlinear scaling of the shear stress with shear rate. Earlier efforts,
including experiments on an attractive system with nanofibers, show that the
interplay between electrostatic repulsion and the van der Waals attraction
governs the degree of fiber flocculation and yielding behavior [33, 34].
Moreover, shear thinning was observed due to adhesive interactions in the
suspension of rigid micro-sized rod-like particles such as polyamide (PA) [35]
and ceramic fibers [36]. However, the role of colloidal interactions remained
unclear. Lately, a theoretical model of aggregated fiber suspensions,
considering the adhesive force between fibers, shows a good match with
experiments on the rigid fiber suspensions [10, 37]. This model included
fitting parameters, which were obtained from experimental data. None of these
explanations describe the behavior of normal stresses, and so their
applicability is limited.

The apparent yield stress, the minimum stress required to begin the flow, is
considered one of the most important rheological properties of fiber
suspensions. The yield stresses in the suspensions rise as the volume fraction
increases, and they are more noticeable for greater aspect ratios [38, 12]. A
recent theoretical model that considered attractive interactions between
fibers in the dilute regime predicted the Bingham law for the shear stress,
with apparent yield stress proportional to the square of volume fraction
$(\phi^{2})$ [37]. The shear stress for suspensions of larger-sized rigid
fibers also follows the Bingham law, but the yield stress rises with higher
power laws in $\phi$ than predicted in previous studies [38]. The origin of
yield stress has been attributed to adhesive contacts even though the fiber
size was large [35]. As pointed out in [38], recent experimental studies
failed to come to a conclusion on whether the fiber size or the attractive
forces are responsible for the yield stress and mentioned the necessity of
further investigation to pinpoint the origin of the yield stress.

As there is no clear explanation of the shear thinning and yield stress in
fiber suspensions, we provide a predictive model in this paper. Our proposed
model incorporates short-range interactions via attractive and repulsive
interactions and contact interactions that quantitatively capture shear
thinning rheology and yield stress in the fiber suspensions, and elucidate the
effect of colloidal interactions. Moreover, our model not only explains the
yield stress and shear thinning behavior but also accurately predicts the
normal stresses, which further strengthens the validity of our model. Lastly,
we demonstrate the versatility of the proposed model by capturing the effect
of changing surface properties by accurately modeling the contact dynamics. To
the best of the authors’ knowledge, this is the first computational study to
quantitatively capture the experimentally observable rheological behavior
(relative viscosity $\eta_{r}$, first normal stress coefficient $\alpha_{1}$,
second normal stress coefficient $\alpha_{2}$, and yield stress $\sigma_{y}$)
for fiber suspensions and demonstrate the underlying physical mechanism.

##  II Simulation methodology

We perform direct numerical simulations of neutrally buoyant fibers with
aspect ratio $AR=l/d$, where $l$ is the length and $d$ is the diameter of the
fiber, in a shear flow generated by the top and bottom walls moving in
opposite directions with a velocity $U_{\infty}=\dot{\gamma}L$ generating an
imposed shear rate $\dot{\gamma}$, as shown in figure 1. Here, $L$ is the
distance between the walls.

Figure 1: Simulation setup of the shear flow of a fiber suspension. The top
and bottom walls move with velocities $U_{\infty}=\dot{\gamma}L$ in the
directions shown by the arrows. $\dot{\gamma}$ is the imposed shear rate, and
$L$ is the distance between the walls.

The suspending fluid is incompressible Newtonian with a viscosity $\eta$, and
its flow is governed by the Navier-Stokes equations.

$\frac{{\partial\mathbf{u}}}{{\partial
t}}+\mathbf{\nabla}\cdot(\mathbf{u}\otimes\mathbf{u})=-\nabla
p+\frac{1}{{{\mathop{\rm\textit{Re}}\nolimits}}}{\mathbf{\nabla}^{2}}\mathbf{u}+\mathbf{f},$
(1) $\mathbf{\nabla}\cdot\mathbf{u}=0,$ (2)

where $\mathbf{u}$ is the velocity field, $p$ is the pressure, $\mathbf{f}$ is
the volume force to account for the suspending fibers, and
$Re=\rho\dot{\gamma}l^{2}/\eta$ is the Reynolds number, where $\rho$ is the
fluid density, $l$ is the characteristic length scale which is also the fiber
length. We model the fibers as in-extensible slender bodies. So, their motion
for the neutrally buoyant case is described by the Euler-Bernoulli beam
equation as [24]:

$\frac{{{\partial^{2}}\mathbf{X}}}{{\partial{t^{2}}}}=\frac{{{\partial^{2}}{\mathbf{X}_{fluid}}}}{{\partial{t^{2}}}}+\frac{\partial}{{\partial
s}}(T\frac{{\partial\mathbf{X}}}{{\partial
s}})-B\frac{{{\partial^{4}}\mathbf{X}}}{{\partial{s^{4}}}}-\mathbf{F}+\mathbf{F}^{f},$
(3)

where $s$ is the curvilinear coordinate along the fiber,
$\mathbf{X}=(x(s,t),y(s,t),z(s,t))$ is the position of the Lagrangian points
on the fiber axis, $T$ is the tension, $B$ is the bending rigidity,
$\mathbf{F}$ is the fluid-solid interaction force, $\mathbf{F}^{f}$ is the net
inter-fiber interaction. The fibers are considered inextensible, expressed as
[39, 40] :

$\frac{{\partial\mathbf{X}}}{{\partial
s}}.\frac{{\partial\mathbf{X}}}{{\partial s}}=1.$ (4)

We use the immersed boundary method (IBM) [41] to couple the motion of fluid
and solid fibers. For the details of the numerical method, readers are
referred to the supplemental material.

In the numerical simulation, the hydrodynamic interactions are well resolved
with the IBM. However, a fine Eulerian mesh is required to capture the short-
range interactions that increase the computational cost. So, we use the
proposed model to calculate the short-range interactions. The short-range
interaction,
$\mathbf{F}^{f}=\mathbf{F}^{lc}+\mathbf{F}^{c}+\mathbf{F}^{cons}$, is split
into the lubrication correction $\mathbf{F}^{lc}$, contact force
$\mathbf{F}^{c}$, conservative force
$\mathbf{F}^{cons}=\mathbf{F}^{A}+\mathbf{F}^{R}$, where $\mathbf{F}^{A}$ is
the van der Waals attractive force, and $\mathbf{F}^{R}$ is the repulsive
force of electrostatic origin. The implementation of the lubrication
correction $\mathbf{F}^{lc}$ can be found in the supplemental material. The
expressions for attractive and repulsive interactions are readily available
from theoretical analyses and previous experimental data [36, 42]. The
attractive force of van der Waals origin acts in the normal directions toward
the fibers and is modeled as $\mathbf{F}^{A}=F_{A}/(h^{2}+H^{2})$, where $h$
is the inter-fiber surface separation. Moreover, H is fixed to 0.01 to prevent
the divergence in $\mathbf{F}^{A}$ when $h\rightarrow 0$ (during contact). The
strength of the attraction is controlled by $F_{A}$, which determines the
values of the attractive force in contact. The repulsive force
$\mathbf{F}^{R}$ also acts in the normal direction to the fibers but is
opposite to the attractive force. $\mathbf{F}^{R}$ decays with the inter-fiber
separation $h$ over the Debye length $\kappa$ as
$\mathbf{F}^{R}=F_{R}\exp(-h/\kappa)$ [20, 25]. The contact between the fibers
occurs when the inter-fiber separation distance $h$ becomes smaller than the
height of surface asperity $h_{R}$ as shown in figure 2. Specifically, the
single-asperity model of the surface roughness has been widely used owing to
its simplicity and effectiveness [43, 22, 44, 45, 46]. Hence, we take the same
approach and model the asperity as a hemispherical bump on the fiber surface.
Actual asperities might not be just hemispherical and can come in various
geometries [47]. However, on average, we can model their behavior by
approximately assuming them hemispherical as routinely done in the tribology
literature [48]. Finally, we split the contact force $\mathbf{F}^{c}$ into the
tangential ($\mathbf{F}_{t}^{C}$) and normal ($\mathbf{F}_{n}^{C}$)
components. The normal contact force is modelled using a Hertz law,
$|\mathbf{F}_{n}^{C}|=-k_{n}|\delta|^{3/2}$, where $\delta=h_{R}-h$ is the
asperity deformation, and $k_{n}$ is the normal stiffness, which is a function
of the fiber material properties [22, 46]. The Coulomb’s friction law gives
the tangential force, $|\mathbf{F}_{t}|=\mu|\mathbf{F}_{n}|$, where $\mu=0.4$
in the current work unless mentioned otherwise. The typical value of $\mu$ is
$0.3-0.5$ as measured experimentally for the polymer fibers [49, 50]. Even
though little is known about the friction mechanics at the nanoscale, the
Coulomb friction allows a correct prediction for flows of colloidal
nanoparticles [51] or hard-sphere suspensions [43]. We use a repulsive force
magnitude as the characteristic force to scale the various forces.
Consequently, the characteristics shear rate scale is
$\dot{\gamma}_{0}={F_{R}/\pi\eta d^{2}}$ and the corresponding shear stress
scale is $\sigma_{0}=F_{R}/\pi d^{2}$.

Figure 2: A sketch of the roughness model, $l$ and $d$ are the length and
diameter of the fiber, respectively, $h_{r}$ is the asperity height, and
$\delta=h-h_{r}$ is the surface overlap. Contact occurs when $\delta\leq 0$.
Dots along the axes of the fibers indicate Lagrangian points.

We calculate the bulk stress ($\Sigma_{ij}$) by volume averaging the viscous
fluid stress and stress generated by the presence of fibers and inter-fiber
interactions. There are three main contributions to the bulk stress: 1) the
hydrodynamic contribution, $\Sigma_{ij}^{h}$ 2) the contact contribution,
$\Sigma_{ij}^{c}$ and 3) the non-contact contribution, $\Sigma_{ij}^{nc}$. The
calculation of bulk stress, including different contributions, is described in
the supplemental material. Rheological properties can be quantified from the
bulk stress, e.g., the relative viscosity, $\eta_{r}=\sigma/\eta\dot{\gamma}$,
where $\sigma$ is the total shear stress in the suspension. The first and
second normal stress coefficients are defined as
$\alpha_{1}(\phi)=(\sigma_{xx}-\sigma_{yy})/\sigma$ and
$\alpha_{2}(\phi)=(\sigma_{yy}-\sigma_{zz})/\sigma$, where $\phi$ is the
volume fraction of fibers defined as $\phi=\frac{{N\pi}}{{4V{(AR)}}^{2}}$,
where $V$ is the volume $(5l\times 8l\times 5l)$ of the simulation cell and
$N$ is the total number of fibers. We present the time-averaged steady-state
values of rheological properties after discarding the initial transient
behavior.

##  III Results and discussion

###  III.1 Shear-rate dependent rheology

We start our analysis by demonstrating the accuracy of the proposed model by
directly comparing the calculated relative viscosity and shear stress with
experimental measurements of neutrally buoyant Polyamide (PA) fiber
suspensions, which exhibit a yield stress and shear thinning viscosity [10,
37]. PA fibers were suspended in a Newtonian fluid (a mixture of UCON oil
75H90000 in distilled water) at different volume fractions ranging from 1% to
17%, depending on the fiber aspect ratio. The roughness ($h_{R}$) of the
fibers has been measured using atomic force microscopy, having values
$h_{R}=5\pm 2nm$ and $h_{R}=14\pm 4nm$ for the aspect ratio, $AR$ = 18 and
$AR$ = 33, respectively. The dimensionless roughness calculated as
$\epsilon_{r}=h_{R}/d$ is 0.0003 and 0.001 for aspect ratios 18 and 33,
respectively. Moreover, the fibers used in the experiment are almost rigid.
The parameters chosen in our simulations mimic these experimental conditions.

The numerical and experimental comparison of the relative viscosity and flow
curves in figures 3a and 3b for the aspect ratios $AR=18$ and $AR=33$ shows
that the proposed computational model does an excellent job capturing the
rate-dependent rheological behavior. We notice that the high shear plateaus of
the relative viscosity shift to higher shear rate values as the volume
fraction increases. The origin of the shear thinning rheology will be
explained in the next section.

Figure 3: Experimental and numerical comparison of relative viscosity
$\eta_{r}$ (red symbol) and shear stress (blue symbol) versus non-dimensional
shear rate at different volume fractions: fibers with aspect ratio (a)
$AR=18$; (b) $AR=33$. Filled symbols with dashed lines show the numerical
data. Hollow symbols correspond to experiments, and solid lines in the shear
stress curve denote the best fit with the Herschel–Bulkley model (equation 5)
for numerical data. Experimental shear rate and shear stress has been scaled
by $\dot{\gamma_{0}}=2.5$ and ${\sigma_{0}}=0.7$, respectively. Inclusion of
attractive and repulsive interactions along with inter-fiber contact
interactions quantitatively reproduce the experimentally observed shear-
thinning viscosity and shear stress [10].

From the computational flow curve (blue symbols) in figure 3, we observe that
our numerical model does reproduce the low shear plateau,
($\frac{\dot{\gamma}}{\dot{\gamma}_{0}}<0.10$), reminiscent of the yield
stress and the linear segments of the flow curve at high dimensionless shear
rates, ($\frac{\dot{\gamma}}{\dot{\gamma}_{0}}>0.10$), observed in
experiments. Due to this trend in the flow curve, the shear stress of the
fiber suspensions can be described by the Herschel–Bulkley model as:

$\frac{\sigma}{\sigma_{0}}=\sigma_{y}+k(\frac{\dot{\gamma}}{\dot{\gamma}_{0}})^{n}$
(5)

where $\sigma_{y}$ is the dimensionless yield stress, $k$ is the consistency
index, and $n<1$ is the shear thinning index. The solid lines in the flow
curve (figure 3) demonstrate the best fit to the Herschel-Bulkley model for
the computational data. The yield stress $\sigma_{y}$ and the fitting
parameter are provided in table 1. The yield stress for a given volume
fraction is higher for fibers with a higher aspect ratio. Furthermore, as the
volume fraction rises, the yield stress increases for a fixed aspect ratio. To
analyze the concentration dependence of the yield stress, we plot the yield
stress for different volume fractions in figure 4a. It is clear that
$\sigma_{y}$ has a weak dependence on $\phi$ in dilute and semi-dilute regimes
($\phi\leq 0.10$ for $AR=10$ and $\phi\leq 0.05$ for $AR$= 18); the dependence
is seen to increase more rapidly with $\phi$ for concentrated regimes. The
dependence of yield stress $\sigma_{y}$ on $\phi$ can be well described using
the relation proposed by [52]:

${\sigma_{y}}={\sigma_{c}}{\left[{\frac{{\frac{\phi}{{\phi_{y}^{0}}}-1}}{{1-\frac{\phi}{{\phi_{y}^{m}}}}}}\right]^{\frac{1}{\beta}}}\
$ (6)

where $\phi_{y}^{0}$ is a threshold volume fraction above which the suspension
exhibits yield stress. In addition, $\phi_{y}^{m}$ denotes the maximum volume
fraction at which the yield stress diverges and beyond which the suspension
does not move at any applied stress. The fitting parameters are provided in
table 2. The threshold volume fraction $\phi_{y}^{0}$ at which the suspension
first exhibits yield stress and maximum volume fraction $\phi_{y}^{m}$ at
which the yield stress diverges decreases with the aspect ratio. Figure 4b
shows the data after rescaling the volume fractions of figure 4a by
$\phi_{y}^{m}$, which leads to the collapse of the data. Hence, the ratio of
$\phi$ to the maximum volume fraction $\phi_{y}^{m}$ at which the yield stress
diverges determines the suspension rheology. Hence, $\phi/\phi_{y}^{m}$ can be
used as a design parameter for tuning the fiber suspension rheology, as it is
easier to measure and control compared to other parameters like size
distribution, roughness, and friction in real-world suspensions. In the next
section, we focus on the origin of shear thinning and yield stress to
understand the underlying mechanism and control the rheological response.

$AR$ | $\phi$ | $\sigma_{y}$ | $k$ | $n$
---|---|---|---|---
$18$ |  1.25 $0.05$1.25 |  1.25$0.1898$1.25 |  1.25$1.409$1.25 |  1.25$0.9975$ 1.25
1.25 $0.11$1.25 |  1.25$0.8571$1.25 |  1.25$2.301$1.25 |  1.25$0.9865$1.25
1.25 $0.17$1.25 |  1.25$2.928$1.25 |  1.25$7.608$1.25 |  1.25$0.9851$1.25
33 |  1.25 $0.05$1.25 |  1.25$0.1903$1.25 |  1.25$1.882$1.25 |  1.25$0.9935$1.25
1.25 $0.09$1.25 |  1.25$1.587$1.25 |  1.25$3.028$1.25 |  1.25$0.9255$1.25
1.25 $0.11$1.25 |  1.25$6.5$1.25 |  1.25$4.58$1.25 |  1.25$0.9218$1.25
Table 1: Herschel–Bulkley parameters (equation 5) for different volume fractions for aspect ratios 18 and 33. $AR$ | $\phi_{y}^{0}$ | $\phi_{y}^{m}$ | $\sigma_{c}$ | $\beta$
---|---|---|---|---
$10$ |  1.25 $0.05$1.25 |  1.25$0.40$1.25 |  1.25$0.17$1.25 |  1.25$0.82$1.25
$18$ |  1.25 $0.03$1.25 |  1.25$0.30$1.25 |  1.25$0.1369$1.25 |  1.25$0.8015$1.25
$33$ |  1.25 $0.017$1.25 |  1.25$0.16$1.25 |  1.25$0.09$1.25 |  1.25$0.7473$1.25
Table 2: Fitting parameters of equation 6 for aspect ratios, $AR$ = 10, 18 and
33.

Figure 4: (a) Yield stress $\sigma_{y}$ as a function of volume fraction. The
solid line shows the best fit of the data to equation 6. (b) Yield stress
$\sigma_{y}$ as a function of re-scaled volume fraction. Re-scaling the volume
fraction $\phi$ by $\phi_{y}^{m}$ leads to the collapse of data indicating
that it is the ratio of $\phi$ to $\phi_{y}^{m}$ that determines the
suspension rheology.

###  III.2 Origin of shear thinning and yield stress

We focus on the origin of yield stress, especially with increasing the
attractive force, in an attempt to acquire a more mechanistic understanding of
the shear-thinning behavior. Figure 5 depicts the effect of varying the
attractive force magnitude on the rheology of a fiber suspension. The aspect
ratio and volume fraction are fixed to 18 and 17%, respectively. We do not
observe yield stress in the suspension for $F_{A}=0$. With an increase in the
magnitude of the attractive force, the slope of the shear thinning curve at
low shear rate increases, which leads to a rise in the yield stress, as shown
in figure 5a. With the increase in $F_{A}$, the separation distance $h$ below
which the conservative force is repulsive decreases (supplemental material
figure S3). At a sufficiently high attractive force, the position of zero
force and cut-off separation for lubrication ($d/4$) are the same, bringing
fibers into direct frictional contact. These frictional contacts can withstand
applied shear stress, increasing the yield stress and viscosity. The role of
colloidal interactions can be better understood by separating the
hydrodynamic, contact, and non-contact (attractive+repulsive) contributions to
viscosity (figures 5b, 5c). The hydrodynamic contribution to the overall
viscosity is negligible for the conditions under consideration. The data
plotted here shows an increase in the contact contribution to the relative
viscosity as the attractive force increases. At low attractive force values,
the non-contact forces contribute the most to the overall viscosity, while
contact contributions take over at higher shear rates. Even though the
attractive force is the cause of the yield stress and must be large enough to
generate the yield stress, we find that as the attractive force increases, the
yield stress shifts from being dominated by the non-contact attraction to
frictional contacts generated by the attractive force. Higher attractive
forces pull the fibers into more direct contact, as shown in figure 5d,
leading to an increase in the yield stress.

Figure 5: (a) Relative viscosity plotted as a function of dimensionless shear
rate for different values of attractive force $F_{A}$. Total relative
viscosity and contributions arising from hydrodynamic, conservative, and
contact forces, plotted as a function of dimensionless shear rate for $AR=18$,
with (b) $F_{A}=0.23F_{R}$ and (c) $F_{A}=0.78F_{R}$. (d) Average number of
fibers in contact as a function of attractive force at
$\dot{\gamma}/\dot{\gamma_{0}}=0.02$. For all cases, the volume fraction is
fixed to $\phi=0.17$. As the attractive force increases, the relative
viscosity shifts from being dominated by the non-contact interactions to the
frictional contacts induced by the attractive force.

Figure 6: The flow-state diagram shown for the volume fraction - attractive
force plane showing no yield stress (light blue), finite yield stress, and
divergence of the yield stress (dark blue). $\phi_{y}^{0}$ and $\phi_{y}^{m}$
lines show the limit below which there is no yield stress and above which the
yield stress diverges, respectively.

The results of this study show that as we change the volume fraction and the
attractive force, the material undergoes different rheological states, which
can be demonstrated in a flow-state diagram. Using the numerical results in
the current work, we generate a flow-state diagram in the $\phi-F_{A}$ plane
for $AR=18$, as shown in figure 6. For the range of attractive forces, the
suspension is in different states for $\phi<\phi_{y}^{0}$,
$\phi_{y}^{0}<\phi<\phi_{y}^{m}$, and $\phi>\phi_{y}^{m}$. Below an attractive
force threshold ($F_{A}^{0}$), no yield stress is observed at any volume
fraction. Above the threshold of attractive force, the system is in the finite
yield stress state for $\phi_{y}^{0}<\phi<\phi_{y}^{m}$. Above $\phi_{y}^{0}$,
the yield stresses increase with $\phi$ and diverge at $\phi_{y}^{m}$.

###  III.3 Normal stress coefficients

In fiber suspensions, normal stress differences inevitably arise due to the
presence of hydrodynamics and inter-fiber interactions [12]. We compute the
first and second normal stress coefficients for $AR$ = 10, 18, and 33, and
compare them with experiments [19] as shown in figure 7. The materials and
experimental conditions are kept the same as the ones used for the viscosity
measurement [10]. Specifically, the experiments are carried out at a higher
shear rate when the viscosity reached a plateau. That is why we fix the shear
rates to $\dot{\gamma}/\dot{\gamma_{0}}=30$ when the suspension viscosity
reaches a plateau in the numerical simulations. All the simulation parameters
are the same as in the last section. While the first normal stress coefficient
$\alpha_{1}$ is positive, the second normal stress coefficient $\alpha_{2}$ is
negative; in quantitative agreement with the experiments [37, 11, 12].
Moreover, the magnitude of $\alpha_{2}$ is smaller than that of $\alpha_{1}$.

In addition, both normal stress coefficients decrease in magnitude with
decreasing the aspect ratio, in agreement with experimental studies [11, 37].
The origin of the normal stress coefficients at a high shear rate can be
better understood by separating the hydrodynamic, contact, and non-contact
(attractive+repulsive) contributions to the normal stress coefficient as shown
in figures 7c and 7d. The hydrodynamic contribution to the normal stress
coefficient is negligible for the conditions under consideration. The data
presented here clearly show that normal stress coefficients are primarily
affected by contact. The normal stress coefficients of suspension increase in
magnitude as the volume fraction rise due to an increase in both the contact
and non-contact contributions. While the non-contact contribution increases
weakly, the contact contribution becomes more significant at a higher volume
fraction.

Figure 7: Numerical and experimental comparison of (a) the first normal stress
coefficient $\alpha_{1}$ and (b) the second normal stress coefficient
$\alpha_{2}$ as a function of volume fractions for $AR$ = 10, 18 and 33.
Hollow circles correspond to experimental data. Contributions arising from
hydrodynamic, conservative, and contact forces to total (c) $\alpha_{1}$, (d)
$\alpha_{2}$ are plotted as a function volume fractions. $\alpha_{1}$ is
positive and $\alpha_{2}$ is negative in agreement with the experiments [12,
11, 19]. The shear rate was fixed to $\dot{\gamma}/\dot{\gamma_{0}}=30$ when
the suspension viscosity reached a plateau. The data plotted here demonstrates
that the contact contribution is the dominant contributor to the normal stress
coefficients at higher shear rates.

Moreover, to investigate the shear rate dependencies, we measure the first and
second normal stress coefficients for $AR$ = 18 and 33, as a function of f
shear rate. We observe a similar shear rate dependence for the normal stresses
as we do for the viscosity, and the data can be found in supplemental material
figure S5.

###  III.4 Effect of coefficient of friction on the normal stress coefficients

Figure 8: (a) First normal stress coefficients $\alpha_{1}$, and (b) second
normal stress coefficient $\alpha_{2}$ as a function of volume fraction $\phi$
for different values of friction coefficient $\mu$. Filled symbols with dashed
lines show the numerical data. The black and red hollow circles show
experimental values for unwashed and washed fibers, respectively. The aspect
ratio of the fibers was fixed to $AR$ = 33. The effect of friction on the
normal stress coefficients is significant, with a rise in $|\alpha_{1}|$ and
$|\alpha_{2}|$. Hydrodynamic, contact and non-contact contributions to (c)
first normal stress coefficient $\alpha_{1}$ and (d) second normal stress
coefficient $\alpha_{2}$ as a function of volume fraction $\phi$ for friction
coefficient $\mu=0.40$ (filled symbols) and $\mu=0.55$ (hollow symbols). As
the volume fraction increases, the contact contribution increases
significantly, leading to an increase in the magnitude of the normal stress
coefficients.

We have provided a quantitative explanation of the yield stress, normal
stresses, and shear thinning in the fiber suspensions. Our model does an
excellent job of reproducing the experimental data. This model can also
capture and explain the effect of changing different parameters that control
rheology. As an example, we will show the capability of the model to quantify
the effect of changing fiber surface properties.

We employed our model to quantify the effect of changing fiber surface
properties and compare the results against the experimental observations by
Bounoua et al. [37]. They changed the surface properties by washing the
fibers. We hypothesize that changing the fiber surface properties will
directly modify the solid contact that can be numerically quantified by
varying the coefficient of friction. We measure the normal stress coefficients
for five different volume functions ($\phi=0.10,0.11,0.13,0.14,0.15$) with
aspect ratio $AR$ = 33, varying the coefficient of friction $\mu$ between 0.30
to 0.55. Numerical results, having $\mu$ = 0.40 and 0.48, match closely with
the experiment for washed and unwashed fibers, respectively, as shown in
figures 8a and 8b . These results demonstrate the applicability of the
proposed model to capture and quantify the effect of modifying fiber surface
properties.

The effect of friction is more pronounced on $|\alpha_{1}|$ compared to
$|\alpha_{2}|$. To understand the underlying mechanism of the effect of
modifying solid contact between fibers, we examine the hydrodynamic, contact,
and non-contact contributions to the normal stress coefficients for $\mu$ =
0.40 (filled symbols) and $\mu$ = 0.55 (hollow symbols) in figures 8c and 8d.
We observe that the hydrodynamics contribution is almost independent of the
friction coefficient. While the non-contact contribution increases weekly, the
contact contribution increases significantly for the first normal stress
coefficient $\alpha_{1}$. The noted increase in the magnitude of $\alpha_{1}$
is, therefore, solely due to the increase in the contact contribution. For the
second normal stress coefficient $\alpha_{2}$ in figure 8d, a similar
conclusion holds: the friction weakly affects the hydrodynamic and non-contact
contributions while it significantly increases the contact contribution.

##  IV Conclusion

Our research provides a fundamental insight into the complex rheological
behavior of fiber suspensions based on balances between hydrodynamic,
conservative, and contact forces. In this work, we provided the first
quantitative explanation of the origin of yield stress, shear thinning, and
normal stress differences in fiber suspensions by focusing on contact and non-
contact contributions. Comparing the relative viscosity, the first normal
stress coefficient $\alpha_{1}$, and the second normal stress coefficient
$\alpha_{2}$ with experimental measurements corroborate the proposed model. We
demonstrated that the attractive interaction of van der Waals origin results
in yield stress. Moreover, we explored the divergence of the yield stress as
the suspension volume fraction, $\phi$ approaches the maximum flowable limit,
$\phi_{y}^{m}$. Re-scaling volume fraction, $\phi$ by the jamming volume
fraction, $\phi_{y}^{m}$ collapsed the yield stress data for varying aspect
ratios on a single curve, denoting that changing the fiber aspect ratio
affects the maximum volume fraction at which the yield stress diverges. In
addition, we demonstrated that the degree of shear thinning and yield stress
depends on the strength of attractive force, which can be regulated in
principle by fiber size, microstructure, chemistry at solid-fluid interfaces,
and fluid and solid phase parameters such as dielectric properties [53, 54].

The first and second normal stress coefficients were compared with experiments
for different aspect ratios and volume fractions to explore the dilute, semi-
dilute, and concentrated regimes. The results showed that the first normal
stress coefficient $\alpha_{1}$ is positive and the second normal stress
coefficient $\alpha_{2}$ is negative, in agreement with available experiments
[12, 19, 11]. As expected, the contact contribution is the dominant
contribution to the normal stress differences at higher shear rates because
non-contact contributions become less and less important as the shear rate
increases.

A direct comparison of normal stress coefficients with the experiment was
performed when the fiber surface was modified. We captured the impact of
changing fiber properties by modifying the frictional contact by changing the
coefficient of friction. Interestingly, the friction appeared to act primarily
through the contact stresses, as the hydrodynamic and non-contact stresses
were unaffected by friction. Our results demonstrated the importance of
accurately modeling the inter-fiber interactions to capture the experimentally
observed shear rate-dependent rheological behavior of fiber suspensions. Due
to the complexity of attractive and friction forces at the microscopic scale
and the number of parameters that are potentially relevant in the surface
physical chemistry, the next step would be to quantitatively determine the
attractive and frictional force from colloidal probe AFM measurements.

ACKNOWLEDGEMENT

AMA would like to acknowledge financial support from the Department of Energy
via grant EE0008910.

## References

  * Bivins  _et al._ [2005] C. H. Bivins, C. Boney, C. Fredd, J. Lassek, P. Sullivan, J. Engels, E. O. Fielder, T. Gorham, T. Judd, A. E. S. Mogollon, _et al._ , New fibers for hydraulic fracturing, Oilfield Review 17, 34 (2005).
  * Elgaddafi  _et al._ [2012] R. Elgaddafi, R. Ahmed, M. George, and F. Growcock, Settling behavior of spherical particles in fiber-containing drilling fluids, Journal of Petroleum Science and Engineering 84, 20 (2012).
  * Hassanpour  _et al._ [2012] M. Hassanpour, P. Shafigh, and H. B. Mahmud, Lightweight aggregate concrete fiber reinforcement–a review, Construction and Building Materials 37, 452 (2012).
  * Lundell  _et al._ [2011] F. Lundell, L. D. Söderberg, and P. H. Alfredsson, Fluid mechanics of papermaking, Annual Review of Fluid Mechanics 43, 195 (2011).
  * Lindström and Uesaka [2008] S. B. Lindström and T. Uesaka, Simulation of semidilute suspensions of non-brownian fibers in shear flow, The Journal of chemical physics 128, 024901 (2008).
  * Nawab and Mason [1958] M. Nawab and S. Mason, Viscosity of dilute suspensions of thread-like particles, The Journal of Physical Chemistry 62, 1248 (1958).
  * Mewis and Metzner [1974] J. Mewis and A. Metzner, The rheological properties of suspensions of fibres in newtonian fluids subjected to extensional deformations, Journal of Fluid Mechanics 62, 593 (1974).
  * Kitano and Kataoka [1981] T. Kitano and T. Kataoka, The rheology of suspensions of vinylon fibers in polymer liquids. i. suspensions in silicone oil, Rheologica Acta 20, 390 (1981).
  * Goto  _et al._ [1986] S. Goto, H. Nagazono, and H. Kato, The flow behavior of fiber suspensions in newtonian fluids and polymer solutions., Rheologica acta 25, 119 (1986).
  * Bounoua  _et al._ [2016a] S. Bounoua, E. Lemaire, J. Férec, G. Ausias, and P. Kuzhir, Shear-thinning in concentrated rigid fiber suspensions: Aggregation induced by adhesive interactions, Journal of Rheology 60, 1279 (2016a).
  * Snook  _et al._ [2014] B. Snook, L. M. Davidson, J. E. Butler, O. Pouliquen, and E. Guazzelli, Normal stress differences in suspensions of rigid fibres, Journal of fluid mechanics 758, 486 (2014).
  * Keshtkar  _et al._ [2009] M. Keshtkar, M. Heuzey, and P. Carreau, Rheological behavior of fiber-filled model suspensions: Effect of fiber flexibility, Journal of Rheology 53, 631 (2009).
  * Quemada and Berli [2002] D. Quemada and C. Berli, Energy of interaction in colloids and its implications in rheological modeling, Advances in colloid and interface science 98, 51 (2002).
  * Raghavan and Douglas [2012] S. R. Raghavan and J. F. Douglas, The conundrum of gel formation by molecular nanofibers, wormlike micelles, and filamentous proteins: gelation without cross-links?, Soft Matter 8, 8539 (2012).
  * Bennington  _et al._ [1990] C. Bennington, R. Kerekes, and J. Grace, The yield stress of fibre suspensions, The Canadian Journal of Chemical Engineering 68, 748 (1990).
  * Song and Youn [2005] Y. S. Song and J. R. Youn, Influence of dispersion states of carbon nanotubes on physical properties of epoxy nanocomposites, Carbon 43, 1378 (2005).
  * Ma  _et al._ [2008] W. Ma, F. Chinesta, A. Ammar, and M. Mackley, Rheological modeling of carbon nanotube aggregate suspensions, Journal of Rheology 52, 1311 (2008).
  * Natale  _et al._ [2014] G. Natale, M. Heuzey, P. Carreau, G. Ausias, and J. Férec, Rheological modeling of carbon nanotube suspensions with rod–rod interactions, AIChE Journal 60, 1476 (2014).
  * Bounoua  _et al._ [2016b] S. Bounoua, P. Kuzhir, and E. Lemaire, Normal stress differences in non-brownian fiber suspensions, Journal of rheology 60, 661 (2016b).
  * Singh  _et al._ [2019] A. Singh, S. Pednekar, J. Chun, M. M. Denn, and J. F. Morris, From yielding to shear jamming in a cohesive frictional suspension, Physical review letters 122, 098004 (2019).
  * Khan  _et al._ [2021] M. Khan, R. V. More, A. A. Banaei, L. Brandt, and A. M. Ardekani, Rheology of concentrated suspension of fibers with load dependent friction coefficient (2021), arXiv:2106.07702 [physics.flu-dyn] .
  * Lobry  _et al._ [2019] L. Lobry, E. Lemaire, F. Blanc, S. Gallier, and F. Peters, Shear thinning in non-brownian suspensions explained by variable friction between particles, Journal of Fluid Mechanics 860, 682 (2019).
  * Wu and Aidun [2010] J. Wu and C. K. Aidun, A numerical study of the effect of fibre stiffness on the rheology of sheared flexible fibre suspensions, Journal of fluid mechanics 662, 123 (2010).
  * Banaei  _et al._ [2020] A. A. Banaei, M. E. Rosti, and L. Brandt, Numerical study of filament suspensions at finite inertia, Journal of Fluid Mechanics 882 (2020).
  * More and Ardekani [2020a] R. V. More and A. M. Ardekani, A unifying mechanism to explain the rate dependent rheological behavior in non-brownian suspensions: One curve to unify them all, arXiv preprint arXiv:2011.07038 (2020a).
  * Guazzelli and Pouliquen [2018] É. Guazzelli and O. Pouliquen, Rheology of dense granular suspensions, Journal of Fluid Mechanics 852 (2018).
  * Djalili-Moghaddam and Toll [2006] M. Djalili-Moghaddam and S. Toll, Fibre suspension rheology: effect of concentration, aspect ratio and fibre size., Rheologica acta 45, 315 (2006).
  * Férec  _et al._ [2009] J. Férec, G. Ausias, M. Heuzey, and P. Carreau, Modeling fiber interactions in semiconcentrated fiber suspensions, Journal of Rheology 53, 49 (2009).
  * Switzer III and Klingenberg [2003] L. H. Switzer III and D. J. Klingenberg, Rheology of sheared flexible fiber suspensions via fiber-level simulations, Journal of Rheology 47, 759 (2003).
  * Maranzano and Wagner [2001] B. J. Maranzano and N. J. Wagner, The effects of interparticle interactions and particle size on reversible shear thickening: Hard-sphere colloidal dispersions, Journal of Rheology 45, 1205 (2001).
  * Gopalakrishnan and Zukoski [2004] V. Gopalakrishnan and C. Zukoski, Effect of attractions on shear thickening in dense suspensions, Journal of Rheology 48, 1321 (2004).
  * Brown  _et al._ [2010] E. Brown, N. A. Forman, C. S. Orellana, H. Zhang, B. W. Maynor, D. E. Betts, J. M. DeSimone, and H. M. Jaeger, Generality of shear thickening in dense suspensions, Nature materials 9, 220 (2010).
  * Solomon and Boger [1998] M. Solomon and D. Boger, The rheology of aqueous dispersions of spindle-type colloidal hematite rods, Journal of Rheology 42, 929 (1998).
  * Michot  _et al._ [2009] L. J. Michot, C. Baravian, I. Bihannic, S. Maddi, C. Moyne, J. F. Duval, P. Levitz, and P. Davidson, Sol- gel and isotropic/nematic transitions in aqueous suspensions of natural nontronite clay. influence of particle anisotropy. 2. gel structure and mechanical properties, Langmuir 25, 127 (2009).
  * Chaouche and Koch [2001] M. Chaouche and D. L. Koch, Rheology of non-brownian rigid fiber suspensions with adhesive contacts, Journal of Rheology 45, 369 (2001).
  * Bergström [1998] L. Bergström, Shear thinning and shear thickening of concentrated ceramic suspensions, Colloids and Surfaces A: Physicochemical and Engineering Aspects 133, 151 (1998).
  * Bounoua  _et al._ [2016c] S. Bounoua, E. Lemaire, J. Férec, G. Ausias, A. Zubarev, and P. Kuzhir, Apparent yield stress in rigid fibre suspensions: the role of attractive colloidal interactions, Journal of Fluid Mechanics 802, 611 (2016c).
  * Tapia  _et al._ [2017] F. Tapia, S. Shaikh, J. E. Butler, O. Pouliquen, and E. Guazzelli, Rheology of concentrated suspensions of non-colloidal rigid fibres, Journal of Fluid Mechanics 827 (2017).
  * Huang  _et al._ [2007] W.-X. Huang, S. J. Shin, and H. J. Sung, Simulation of flexible filaments in a uniform flow by the immersed boundary method, Journal of computational physics 226, 2206 (2007).
  * Pinelli  _et al._ [2017] A. Pinelli, M. Omidyeganeh, C. Brücker, A. Revell, A. Sarkar, and E. Alinovi, The pelskin project: part iv—control of bluff body wakes using hairy filaments, Meccanica 52, 1503 (2017).
  * Peskin [1972] C. S. Peskin, Flow patterns around heart valves: a numerical method, Journal of computational physics 10, 252 (1972).
  * Israelachvili [2011] J. N. Israelachvili,  _Intermolecular and surface forces_ (Academic press, 2011).
  * Gallier  _et al._ [2014] S. Gallier, E. Lemaire, F. Peters, and L. Lobry, Rheology of sheared suspensions of rough frictional particles, Journal of Fluid Mechanics 757, 514 (2014).
  * More and Ardekani [2020b] R. More and A. Ardekani, A constitutive model for sheared dense suspensions of rough particles, Journal of Rheology 64, 1107 (2020b).
  * More and Ardekani [2020c] R. More and A. Ardekani, Roughness induced shear thickening in frictional non-brownian suspensions: A numerical study, Journal of Rheology 64, 283 (2020c).
  * More and Ardekani [2020d] R. V. More and A. M. Ardekani, Effect of roughness on the rheology of concentrated non-brownian suspensions: A numerical study, Journal of Rheology 64, 67 (2020d).
  * Tanner and Dai [2016] R. I. Tanner and S. Dai, Particle roughness and rheology in noncolloidal suspensions, Journal of Rheology 60, 809 (2016).
  * Broedersz and MacKintosh [2014] C. P. Broedersz and F. C. MacKintosh, Modeling semiflexible polymer networks, Reviews of Modern Physics 86, 995 (2014).
  * Bowden and Tabor [1964] F. P. Bowden and D. Tabor,  _The Friction and Lubrication of Solids-Part II._ (Oxford, England, University Press, 1964).
  * Petrich and Koch [1998] M. P. Petrich and D. L. Koch, Interactions between contacting fibers, Physics of Fluids 10, 2111 (1998).
  * Fujita and Yamaguchi [2008] M. Fujita and Y. Yamaguchi, Simulation model of concentrated colloidal nanoparticulate flows, Physical Review E 77, 026706 (2008).
  * Zhou  _et al._ [1995] J. Z. Zhou, P. H. Uhlherr, and F. T. Luo, Yield stress and maximum packing fraction of concentrated suspensions, Rheologica acta 34, 544 (1995).
  * Scales  _et al._ [1998] P. J. Scales, S. B. Johnson, T. W. Healy, and P. C. Kapur, Shear yield stress of partially flocculated colloidal suspensions, AIChE Journal 44, 538 (1998).
  * Gálvez  _et al._ [2017] L. O. Gálvez, S. de Beer, D. van der Meer, and A. Pons, Dramatic effect of fluid chemistry on cornstarch suspensions: Linking particle interactions to macroscopic rheology, Physical Review E 95, 030602 (2017).


with the two sides switched everywhere.

If S makes an $X$-move, it can either be a return to $u$ or to a vertex $s$
below $u$. Assume that the return is made to $u$ and that $X$ is a
$\mu$-variable. Now the models in $\mathrm{left}(v)$ are moved to
$\mathrm{left}^{\prime}(u)$ so any that were already in $\mathbb{L}(u)$ stay
there and more may be added so $\mathbb{L}(u)\subseteq\mathbb{L}^{\prime}(u)$.
The models in $\mathrm{right}(v)$ cease to be relevant in
$\mathrm{right}^{\prime}(u)$ but they remain as old models in
$\mathrm{right}^{\prime}(v)$ and are still counted for
$\mathbb{R}^{\prime}(v)$ in $P^{\prime}$ as they were in $P$ so
$\mathbb{R}^{\prime}(u)=\mathbb{R}(u)$. The case of a $\nu$ is symmetric with
the two sides switched.

Finally assume that an $X$-move is made returning to a vertex $s$ below $u$.
Assume again that $X$ is a $\mu$-variable. The models in $\mathrm{left}(v)$
are moved to $\mathrm{left}^{\prime}(s)$ so
$\mathbb{L}(s)\subseteq\mathbb{L}^{\prime}(s)$. Everything not below $s$
remains unchanged so $\mathbb{L}(u)\subseteq\mathbb{L}^{\prime}(u)$. The
models in $\mathrm{right}(v)$ cease to be relevant in
$\mathrm{right}^{\prime}(s)$ but they remain as old models in
$\mathrm{right}^{\prime}(v)$ and are still counted for
$\mathbb{R}^{\prime}(s)$ and therefore also for $\mathbb{R}^{\prime}(u)$ just
like in $P$. Thus $\mathbb{R}^{\prime}(u)=\mathbb{R}(u)$. The case of a $\nu$
is again symmetric. ∎

We finally have all of the required notation and lemmas to show that the non-
elementary succinctness gap is present also between $\mathrm{FO}$ and modal
$\mu$-calculus.

######  Theorem 6.6.

First-order logic is non-elementarily more succinct than modal $\mu$-calculus.

###### Proof.

We prove an analogous result to Lemma 4.9 for the $\mathrm{L}_{\mu}$ game. We
use the notation $\mathcal{G}(\mathbb{V},\mathbb{E})$ for the same graph as in
subsection 4.3. The precise statement we prove is as follows:

Let $n\in\mathbb{N}$ and let $k_{0}\in\mathbb{N}$. If
$k_{0}<\log(\chi(\mathcal{G}(\mathbb{C}_{n},\mathbb{D}_{n})))$, then D has a
winning strategy in the game
$\mu\textnormal{-}\mathrm{FS}^{\Phi}_{k_{0}}(\mathbb{C}_{n},\mathbb{D}_{n})$.

In this proof we only consider relevant models. Many positions in the game
$\mu\textnormal{-}\mathrm{FS}^{\Phi}_{k_{0}}(\mathbb{C}_{n},\mathbb{D}_{n})$
have also non-relevant models but they are not needed for the strategy of D we
describe here and can safely be ignored. We will assume all models are
relevant and occasionally comment on why models remain or cease to be relevant
after some moves of the game.

We show by induction that D has a strategy to maintain the following condition
in any position
$P=(V,E,B,\mathrm{lab},\mathrm{res},\mathrm{left},\mathrm{right},v)$:

$\mathrm{res}(v)<\log(\chi(\mathcal{G}(\mathbb{L}(v),\mathbb{R}(v)))$ ($\ast$)

At the start of the game, $\ast$ ‣ 6 holds by assumption. Since the sets
$\mathbb{C}_{n}$ and $\mathbb{D}_{n}$ are already finite, D can keep the full
sets for the first move of the game.

We first show that if S ever makes a modal move while $\ast$ ‣ 6 holds, D gets
a winning strategy for the game. We assume $v\notin\mathrm{dom}(\mathrm{lab})$
since the first modal move in a game must always be made in an unlabelled
vertex. In this case there are no other vertices below $v$ so
$\mathbb{L}(v)=\mathrm{left}(v)$ and $\mathbb{R}(v)=\mathrm{right}(v)$. We
assume $\mathrm{res}(v)>1$ so that S can make a modal move and not lose
immediately due to the resource running out. From $\ast$ ‣ 6 we obtain
$\chi(\mathcal{G}(\mathrm{left}(v),\mathrm{right}(v)))>2$ so there are
relevant clocked models
$(\textstyle{\bigtriangleup}(\mathcal{M}_{1},w_{1}),c_{1},\mathsf{new}),(\textstyle{\bigtriangleup}(\mathcal{M}_{2},w_{2}),c_{2},\mathsf{new})\in\mathrm{left}(v)$
and
$(\textstyle{\bigtriangleup}\\{(\mathcal{M}_{1},w_{1}),(\mathcal{M}_{2},w_{2})\\},c_{3},\mathsf{new})\in\mathrm{right}(v)$.
Now if S makes a $\Diamond$\- or $\square$-move, then in the following
position $P^{\prime}$ there is $i\in\\{1,2\\}$ s.t.
$(\mathcal{M}_{i},w_{i},c_{i},\mathsf{new})\in\mathrm{left}(v^{\prime})$ and
$(\mathcal{M}_{i},w_{i},c_{3},\mathsf{new})\in\mathrm{right}(v^{\prime})$. As
these two share the same underlying pointed model they are bisimilar and
moreover, since the depth has decreased by at least $1$ from the previous
position, the models are strictly relevant. By Lemma 6.3, D now has a winning
strategy from position $P^{\prime}$.

If S makes a $\lor$-move, let $v_{1}$ and $v_{2}$ be the successors of $v$. In
the following position, whichever it may be, we have
$\mathbb{L}(v)=\mathbb{L}(v_{1})\cup\mathbb{L}(v_{2})$ and
$\mathbb{R}(v)=\mathbb{R}(v_{1})\cap\mathbb{R}(v_{2})$. Let
$\mathcal{G}_{s}=\mathcal{G}(\mathbb{L}(s),\mathbb{R}(s))=(V_{s},E_{s})$ for
$s\in\\{v,v_{1},v_{2}\\}$. We obtain $V_{v}=V_{v_{1}}\cup V_{v_{2}}$ and
$E_{v}\upharpoonright V_{v_{i}}\subseteq E_{v_{i}}$ for $i\in\\{1,2\\}$. By
Lemma 4.8,

$\chi(\mathcal{G}_{v})\leq\chi(V_{v_{1}},E_{v}\upharpoonright
V_{v_{1}})+\chi(V_{v_{2}},E_{v}\upharpoonright
V_{v_{2}})\leq\chi(\mathcal{G}_{v_{1}})+\chi(\mathcal{G}_{v_{2}}).$

Thus (just like in the proof of Theorem 4.9) we obtain
$\mathrm{res}(v_{i})<\log(\chi(\mathcal{G}_{v_{i}})$ for some $i\in\\{1,2\\}$
so $\ast$ ‣ 6 holds in the following position after D chooses this $i$.

If S makes a $\land$-move, let $v_{1}$ and $v_{2}$ be the successors of $v$.
In the following position we have
$\mathbb{L}(v)=\mathbb{L}(v_{1})\cap\mathbb{L}(v_{2})$ and
$\mathbb{R}(v)=\mathbb{R}(v_{1})\cup\mathbb{R}(v_{2})$. We use the notation
$\mathcal{G}_{s}=(V_{s},E_{s})$ from the previous case and obtain
$V_{v}=V_{v_{1}}\cap V_{v_{2}}$ and $E_{v}=(E_{v_{1}}\upharpoonright
V_{v})\cup(E_{v_{2}}\upharpoonright V_{v})$. By Lemma 4.8,

$\chi(G_{v})\leq\chi(V_{v},E_{v_{1}}\upharpoonright
V_{v})\chi(V_{v},E_{v_{2}}\upharpoonright
V_{v})\leq\chi(\mathcal{G}_{v_{1}})\chi(\mathcal{G}_{v_{2}}).$

Thus we again obtain $\mathrm{res}(v_{i})<\log(\chi(\mathcal{G}_{v_{i}})$ for
some $i\in\\{1,2\\}$ so $\ast$ ‣ 6 holds in the following position after D
chooses this $i$.

If S makes a $\eta X$-move, where $\eta\in\\{\mu,\nu\\}$, then D sets her
clock for each model at the same value as the depth of the model. All relevant
models remain relevant and $\ast$ ‣ 6 is maintained.

If S makes an $X$-move, by $\ast$ ‣ 6, S does not immediately win the game.
Assume that $u$ is the vertex returned to and $P^{\circ}$ is the previous
position when $u$ was the current vertex. Let $P^{\prime}$ be the position
after this $X$-move. By Lemma 6.5 we obtain
$\mathbb{L}^{\circ}(u)\subseteq\mathbb{L}^{\prime}(u)$ and
$\mathbb{R}^{\circ}(u)\subseteq\mathbb{R}^{\prime}(u)$. By induction
hypothesis, $\ast$ ‣ 6 held in $P^{\circ}$ and clearly
$\mathrm{res}^{\prime}(u)=\mathrm{res}^{\circ}(u)$, so $\ast$ ‣ 6 still holds
in $P^{\prime}$.

If S makes a $\mathit{Lit}$-move, by $\ast$ ‣ 6,
$\mathrm{left}(v)\neq\emptyset$ and $\mathrm{right}(v)\neq\emptyset$. Since
all the models are propositionally equivalent, D wins the game.

By Theorem 5.6, we obtain that there is no sentence
$\varphi\in\mathrm{L}_{\mu}(\emptyset)$ that separates $\mathbb{C}_{n}$ from
$\mathbb{D}_{n}$ with
$\mathrm{sz}(\varphi)\leq\log(\chi(\mathcal{G}(\mathbb{C}_{n},\mathbb{D}_{n})))=\mathrm{twr}(n-1)$.
Thus $\mathrm{FO}$ is non-elementarily more succinct than $\mathrm{L}_{\mu}$.

∎

######  Remark 6.7.

Just like in the case of $\mathrm{ML}$, we remark that the result of Theorem
6.6 also holds for DAG-size. This is again because the difference between the
size of an $\mathrm{L}_{\mu}$ formula in our sense and the DAG-size of the
same formula is at most exponential.

##  7 Conclusion

We have defined formula size games for basic modal logic and modal
$\mu$-calculus. The games utilize resource parameters to achieve a truly two-
player game. In the case of modal logic the players only construct one branch
of the game tree. This is in contrast with the original Adler-Immerman game,
where the players form the whole tree in a single play. For modal
$\mu$-calculus the recursive nature of fixed point operators necessitates
returning to previously visited nodes. However, the game still traverses only
one possible path through the game tree in a single play and some branches can
remain unvisited for the entire play. The $\mu$-calculus game has infinite
branching but the use of decreasing ordinal clocks, as in [16], makes each
play of the game finite.

We used the games to show that the property “all successor models are
$n$-bisimilar with each other” cannot be defined in basic modal logic or modal
$\mu$-calculus with a formula of size less than the exponential tower of
height $n-1$. On the other hand this property can be defined in $\mathrm{FO}$
with a formula of size linear in $n$. This means that $\mathrm{FO}$ is non-
elementarily more succinct than both $\mathrm{ML}$ and $\mathrm{L}_{\mu}$. We
also show that the same property can be defined in two-dimensional modal logic
$\mathrm{ML}^{2}$ with a formula of size exponential in $n$. Thus the non-
elementary succinctness gap is also present between $\mathrm{ML}^{2}$ and both
$\mathrm{ML}$ and $\mathrm{L}_{\mu}$

We find the $\mathrm{ML}$-game to be a useful tool for proving lower bounds on
the size of $\mathrm{ML}$-formulas. Depending on the desired result, the game
can also be modified to count a more specific parameter instead of formula
size, such as the number or nesting depth of a specific operator.

The $\mu$-calculus game is also functional for proving lower bounds, but with
some caveats. The main theorem stating the usefulness of the game, Theorem
5.6, requires uniform strategies for S. This means that we assume S has a
single formula in mind and always plays according to that formula. It may be
that this restriction could be removed but we have been unable to prove this.
However, to show succinctness results we only need one direction of the
equivalence so the issue is usually not relevant in practice. The greater
concern is whether the game can be successfully used to prove succinctness
results for $\mu$-calculus in more complicated contexts. Here we only
generalize a result already obtained with the $\mathrm{ML}$ game and we have
so far failed to produce any other results with the $\mathrm{L}_{\mu}$ game
due to its sheer complexity. A question related to this difficulty would be
whether the game could be simplified significantly while still preserving its
functionality. It would be especially interesting to apply the game to open
problems related to $\mu$-calculus and succinctness, such as whether there is
a polynomial transformation from $\mathrm{L}_{\mu}$ to the guarded fragment or
from vectorial form to regular $\mathrm{L}_{\mu}$ [4].

## References

  * [1] M. Adler and N. Immerman.  An _n!_ lower bound on formula size.  ACM Trans. Comput. Log., 4(3):296–314, 2003.
  * [2] P. Blackburn, M. de Rijke, and Y. Venema.  Modal Logic.  Cambridge University Press, New York, NY, USA, 2001.
  * [3] J. Bradfield and C. Stirling.  Chapter 12 - modal mu-calculi.  In P. Blackburn, J. V. Benthem, and F. Wolter, editors, Handbook of Modal Logic, volume 3 of Studies in Logic and Practical Reasoning, pages 721 – 756. Elsevier, 2007.
  * [4] F. Bruse, O. Friedmann, and M. Lange.  On guarded transformation in the modal mu-calculus.  Logic Journal of the IGPL, 23(2):194 – 216, 2015.
  * [5] K. Etessami, M. Y. Vardi, and T. Wilke.  First-order logic with two variables and unary temporal logic.  Inf. Comput., 179(2):279–295, 2002.
  * [6] K. Etessami and T. Wilke.  An until hierarchy and other applications of an ehrenfeucht-fraïssé game for temporal logic.  Inf. Comput., 160(1-2):88–108, 2000.
  * [7] D. Fernández-Duque and P. Iliev.  Succinctness in subsystems of the spatial mu-calculus.  ArXiv e-prints, Aug. 2017.
  * [8] T. French, W. van der Hoek, P. Iliev, and B. P. Kooi.  Succinctness of epistemic languages.  In IJCAI 2011, Proceedings of the 22nd International Joint Conference on Artificial Intelligence, Barcelona, Catalonia, Spain, July 16-22, 2011, pages 881–886, 2011.
  * [9] T. French, W. van der Hoek, P. Iliev, and B. P. Kooi.  On the succinctness of some modal logics.  Artif. Intell., 197:56–85, 2013.
  * [10] D. Gabbay and V. Shehtman.  Products of modal logics, part 1.  Logic Journal of the IGPL, 6(1):73–146, 1998.
  * [11] D. Gabbay and V. Shehtman.  Products of modal logics. part 2: Relativised quantifiers in classical logic.  Logic Journal of the IGPL, 8(2):165–210, 2000.
  * [12] D. Gabbay and V. Shehtman.  Products of modal logics. part 3: Products of modal and temporal logics.  Studia Logica, 72(2):157–183, 2002.
  * [13] G. Gogic, H. Kautz, C. Papadimitriou, and B. Selman.  The comparative linguistics of knowledge representation.  In Proceedings of the 14th International Joint Conference on Artificial Intelligence - Volume 1, IJCAI’95, pages 862–869, San Francisco, CA, USA, 1995. Morgan Kaufmann Publishers Inc.
  * [14] M. Grohe and N. Schweikardt.  The succinctness of first-order logic on linear orders.  Logical Methods in Computer Science, 1(1), 2005.
  * [15] Grohe, Martin and Schweikardt, Nicole.  Comparing the succinctness of monadic query languages over finite trees.  RAIRO-Theor. Inf. Appl., 38(4):343–373, 2004.
  * [16] L. Hella, A. Kuusisto, and R. Rönnholm.  Bounded game-theoretic semantics for modal mu-calculus.  ArXiv e-prints, June 2017.
  * [17] L. Hella and J. Väänänen.  The size of a formula as a measure of complexity.  In Logic Without Borders - Essays on Set Theory, Model Theory, Philosophical Logic and Philosophy of Mathematics, pages 193–214. De Gruyter, 2015.
  * [18] L. Hella and M. Vilander.  The succinctness of first-order logic over modal logic via a formula size game.  In Advances in Modal Logic 11, Budapest, Hungary, August 30 - September 2, 2016, pages 401–419, 2016.
  * [19] G. Lenzi.  The modal mu-calculus: A survey.  Task Q, 9:293–316, 2005.
  * [20] C. Lutz.  Complexity and succinctness of public announcement logic.  In 5th International Joint Conference on Autonomous Agents and Multiagent Systems (AAMAS 2006), Hakodate, Japan, May 8-12, 2006, pages 137–143, 2006.
  * [21] C. Lutz, U. Sattler, and F. Wolter.  Modal logic and the two-variable fragment.  In Computer Science Logic, Paris, France, September 10-13, 2001, Proceedings, pages 247–261, 2001.
  * [22] N. Markey.  Temporal logic with past is exponentially more succinct, concurrency column.  Bulletin of the EATCS, 79:122–128, 2003.
  * [23] M. Marx and Y. Venema.  Multi-dimensional modal logic, volume 4 of Applied Logic Series.  Kluwer Academic Publishers, Dordrecht, 1997.
  * [24] J. C. McCabe-Dansted, T. French, S. Pinchinat, and M. Reynolds.  Expressiveness and succinctness of a logic of robustness.  Journal of Applied Non-Classical Logics, 25(3):193–228, 2015.
  * [25] M. Otto.  Bisimulation invariance and finite models.  In Z. Chatzidakis, P. Koepke, and W. Pohlers, editors, Logic Colloquium ’02, Lecture Notes in Logic, pages 276–298. Cambridge University Press, 2006.
  * [26] A. P. Sistla and E. M. Clarke.  The complexity of propositional linear temporal logics.  J. ACM, 32(3):733–749, 1985.
  * [27] L. J. Stockmeyer.  The Complexity of Decision Problems in Automata Theory and Logic.  PhD thesis, Massachusetts Institute of Technology, 1974.
  * [28] W. van der Hoek and P. Iliev.  On the relative succinctness of modal logics with union, intersection and quantification.  In International conference on Autonomous Agents and Multi-Agent Systems, AAMAS ’14, Paris, France, May 5-9, 2014, pages 341–348, 2014.
  * [29] W. van der Hoek, P. Iliev, and B. P. Kooi.  On the relative succinctness of two extensions by definitions of multimodal logic.  In How the World Computes - Turing Centenary Conference and 8th Conference on Computability in Europe, CiE 2012, Cambridge, UK, June 18-23, 2012\. Proceedings, pages 323–333, 2012.
  * [30] H. van Ditmarsch, J. Fan, W. van der Hoek, and P. Iliev.  Some exponential lower bounds on formula-size in modal logic.  In Advances in Modal Logic 10, held in Groningen, The Netherlands, August 5-8, 2014, pages 139–157, 2014.
  * [31] H. van Ditmarsch and P. Iliev.  The succinctness of the cover modality.  Journal of Applied Non-Classical Logics, 25(4):373–405, 2015.
  * [32] T. Wilke.  CTL${}^{\mbox{+}}$ is exponentially more succinct than CTL.  In Foundations of Software Technology and Theoretical Computer Science, 19th Conference, Chennai, India, December 13-15, 1999, Proceedings, pages 110–121, 1999.


# Variational Quantum Algorithms for Chemical Simulation and Drug Discovery

Hasan Mustafa  AI@Scale, Fractal Analytics
Mumbai, India
<EMAIL_ADDRESS>Sai Nandan Morapakula  AI@Scale, Fractal Analytics
Bangalore, India
<EMAIL_ADDRESS>Prateek Jain  Srinjoy Ganguly  AI@Scale, Fractal
Analytics
Gurugram, India
<EMAIL_ADDRESS>AI@Scale, Fractal Analytics
Gurugram, India
<EMAIL_ADDRESS>
###### Abstract

Quantum computing has gained a lot of attention recently, and scientists have
seen potential applications in this field using quantum computing for
Cryptography and Communication to Machine Learning and Healthcare. Protein
folding has been one of the most interesting areas to study, and it is also
one of the biggest problems of biochemistry. Each protein folds distinctively,
and the difficulty of finding its stable shape rapidly increases with an
increase in the number of amino acids in the chain. A moderate protein has
about 100 amino acids, and the number of combinations one needs to verify to
find the stable structure is enormous. At some point, the number of these
combinations will be so vast that classical computers cannot even attempt to
solve them. In this paper, we examine how this problem can be solved with the
help of quantum computing using two different algorithms, Variational Quantum
Eigensolver (VQE) and Quantum Approximate Optimization Algorithm (QAOA), using
Qiskit Nature. We compare the results of different quantum hardware and
simulators and check how error mitigation affects the performance. Further, we
make comparisons with SoTA algorithms and evaluate the reliability of the
method.

###### Index Terms:

Protein Folding, Quantum Computing, Variational Algorithms, VQE, QAOA

##  I Introduction

Proteins are chains of smaller molecules called amino acids. Amino acids bond
with each other by making a so-called Peptide bond. The median length of
proteins in terms of the number of amino acids is 375 for Humans. These
incredibly long protein chains fold and form complex 3-D structures that are
stable. Finding this stable 3-D shape of a particular protein given the amino
acid sequence is known as the protein folding problem. Understanding how
proteins fold is essential as the function of a protein is correlated with its
shape. The 3-D structure of the protein is entirely determined by the amino
acid sequence. Thus, mutations which can lead to a change in the amino acid
sequence of a protein cause the protein to fold differently and render it
incapable of performing the desired function. The effects of a misfolded
protein can range from minor diseases to severe illnesses. Some of the
diseases which are caused by improperly folded proteins are, Type 2 Diabetes,
Alzheimer’s, Parkinson’s, Amyotrophic Lateral Sclerosis(ALS). The amino acid
sequence considered resembles a protein that interacts with Alzheimer’s
disease-causing enzyme.

Drug development involves simulating the behaviour of the proteins and their
interaction with drug molecules. Being able to simulate the 3-D structure of
protein reliably would considerably speed up the drug discovery process, which
currently relies on experimental methods that are expensive and time-
consuming. The protein folding problem is computationally very expensive and
scales exponentially in the number of amino acids. This complexity has led to
the development of new techniques such as those using machine learning. The
most successful of them being AlphaFold, which outperformed other techniques
in benchmarking competition [1].

Through the paper, we provide a comprehensive analysis of the VQE in contrast
to the SoTA AlphaFold model. This is one of the first of its kind benchmarking
result for a particular application on different hardware. We also identify
the caveats in the existing approach and quantify them.

In this paper, we examine the application of Variational Quantum Algorithms to
the protein folding problem. We use the protein folding module developed by
Qiskit [2] to model the problem. Variational algorithms are well suited for
the Noisy intermediate-scale quantum(NISQ) era, which is characterized by
noisy hardware capable of running low-depth circuits.

The paper is organized as follows. Section II outlines the protein folding
model and the quantum algorithms we use. In section III we describe the
specific setting we have used such as the optimizers and ansatz. We describe
all the hardware and simulator results of various types of quantum computers
and compare them in section IV. Section V comments on the results. Finally we
give additional comments and summary in Section VI.

##  II Literature Review

###  II-A Protein Folding Problem

The protein folding problem seeks to answer the question of how a protein’s
3-D structure can be determined from the amino acid sequence [3]. This can be
modelled as a search problem in the space of all possible conformations. The
protein folding problem is characterized by an infinitely large search space.
The solution to this problem represents the functional structure of the
protein that matches the experimentally determined structure. The number of
possible conformers can go as high as 100300. The enormity of the search space
makes it NP-Hard. It is well established that protein folds into a
thermodynamically stable conformation [3]. Thus, the goal of the search
problem can be equivalently set to finding the protein structure with the
minimum energy[4].

###  II-B Protein Folding Model

The peptide is abstracted as a sequence of beads on a tetrahedral lattice.
Each point in a tetrahedral lattice is connected to 4 other points; thus the
chain can ”turn” in one of the four directions when the following amino acid
is added. The direction in which the chain turns with the addition of each
amino acid from the sequence is used to describe a protein conformation. The
amino acid chain is encoded by assigning 4 qubits per bead, each qubit
representing one of the four possible turns. Without loss of generality, the
first two turns can be chosen arbitrarily. The remaining chain can then be
encoded into 4(N-3) qubits. The representation was further condensed to reduce
the number of qubits used to 2(N-3). The interaction between the amino acids
is calculated based on inter-residue contact energies evaluated by Miyazawa-
Jernigan [5]. Interaction qubits are separately used to encode the interaction
energy. Three penalty terms, penalty chiral, penalty back, penalty 1, are also
introduced to prevent the chain from folding back to itself, impose chirality
constraints, and penalize local overlap between beads within a nearest
neighbour contact. A Hamiltonian encoding the interaction energies and the
penalties for which the ground state is to be determined is generated.

$H(\mathbf{q})=H_{gc}(\mathbf{q}_{cf})+H_{ch}(\mathbf{q}_{cf})+H_{in}(\mathbf{q}_{cf},\mathbf{q}_{in})$

  * •

$H_{gc}$ is the geometrical constraint term (governing the growth of the
primary sequence of amino acids without bifurcations)

  * •

$H_{ch}$ is the chirality constraint (enforcing the right stereochemistry for
the system)

  * •

$H_{in}$ is the interaction energy terms of the system. In our case we
consider only nearest neighbor interactions.

Two of the most common variational algorithms that achieve the objective are
VQE and QAOA. A subsequent modification of these algorithms is to use CVaR [6]
as the objective function instead of simply the expectation value.

###  II-C VQE

Variational Quantum Eigensolver(VQE) is a hybrid classical-quantum algorithm
that finds the eigenvalue of the ground state of a Hamiltonian. VQEs are based
on the variational principle. VQEs involve estimating the expectation value of
the Hamiltonian for a parametric ansatz state and then minimizing the
expectation value by using a classical optimizer to tune the parameters. The
ansatz state might be chosen based on the physical properties of the system
[7]. Quantum Phase Estimation(QPE) algorithm can also be used for the same
purpose, but its implementation is not feasible in the NISQ Era and is thus
not considered in this paper.

Figure 1: Working of VQE

Figure 1 illustrates the workflow of VQE as split between the classical and
quantum computers. We use classical computers/optimizers to adjust and tune
the parameters. VQE has three main components, Ansatz, Hamiltonian, and
Measurement. An ansatz is the parametrized quantum circuit that generates the
quantum state we need in order to apply the variational algorithm. It should
be noted that an ansatz is not fixed; the design of the ansatz depends on the
problem and hardware; the better the range and its shallowness the better is
the ansatz. Hamiltonian, denoted as a matrix or a sum of Pauli operators.

###  II-D QAOA

Quantum Approximate optimization algorithm(QAOA) was proposed to solve
combinatorial optimizations. The optimization problem can be remodelled into
an equivalent problem of finding the ground state of the corresponding
Hamiltonian. The variational form chosen for QAOA is[8]

$U(C,\gamma_{1})U(M,\beta_{1})..ptimes...U(C,\gamma_{p})U(M,\beta_{p})\ket{\Psi}$
(1) $U(C,\gamma)=e^{-iH_{C}\gamma}$ (2) $U(M,\beta)=e^{-iH_{B}\beta}$ (3)

Where $H_{C}$ and $H_{B}$ are the cost and mixing Hamiltonian respectively.
The former encodes the optimization problem to be solved and the latter is
required in case the initial state is an eigenstate of $H_{C}$. In the absence
of a mixing Hamiltonian, the algorithm might get stuck in one of the
eigenstates of $H_{C}$. The optimal parameters
$\gamma_{1}^{*},\beta_{1}^{*},....\gamma_{p}^{*},\beta_{p}^{*}$ are found
using a classical optimization scheme so as to minimize the cost function.

###  II-E CVaR VQE and QAOA

Conditional Value at Risk(CVaR), frequently used in finance, is a measure of
risk associated with an investment. More precisely it quantifies the expected
loss associated with a certain level of risk. It was shown that using CVaR
expectation value as the loss function leads to a greater chance of getting a
state that overlaps with the ground state when used with VQE and QAOA
algorithms.[9] CVaR expectation value is calculated by averaging the energies
of the lowest $\alpha$ % of sampled energies. Suppose, the energies $E_{i}$
are arranged in ascending order then CVaR expectation value is given by,

$CVaR=\frac{1}{\frac{\alpha N}{100}}\sum_{i=0}^{i=\alpha N}E_{i}$ (4)

where $n$ is the number of shots i.e. the number of energies measured.

###  II-F Significance of protein selected

The amino acid sequence chosen ”YPYFIP” is a subsequence of the peptide chain
that binds with the exosite of BACE1 Alzheimer’s enzyme [10]. BACE1 has been
shown to produce beta-amyloid which is responsible for symptoms like memory
loss in Alzheimer’s. We have removed the Acetylene and NH2 chains and only
considered the main peptide.

###  II-G Alphafold

Alphafold is a state of the art deep learning model developed by Deepmind that
predicts the 3-D structure of a protein based on the amino acid sequence. The
model uses biological, evolutionary and geometric constraint of protein
structures. This approch is quite different from the one used by Qiskit
nature, the subject of this study where the approach is to minimize the
interaction energy. AlphaFold completely outperformed other techniques in
CASP14, an annual competition that is considered the gold standard for protein
folding [1].

##  III Methodology and Setup

###  III-A Hardware

We run our experiments on quantum devices provided by IBM and Oxford Quantum
Computing(OQC). Further, we tried using Rigetti and IONQ hardware but ran into
compatibility issues and thus had to abandon them. The IBM, Rigetti, and OQC
devices are based on superconducting qubits, while the IONQ device is based on
Trapped Ion technology [11]. Another notable feature of IONQ QPU is that it
has a fully connected topology thus any two-qubit gate can be directly applied
to any qubit pair.

We alternatively also run the experiments on simulators to compare and
contrast the results obtained on real hardware. We use the SV1 simulator
provided by AWS braket, and IBM FakeBrooklyn provided by Qiskit. SV1 is a
state vector simulator while FakeBrooklyn is a mock hardware with noise.

|
---|---
|

OQC LUCY

(a)

(b)

IBM Perth

(c)

(d)

Figure 2: Hardware Architecture of OQC and IBM Perth

###  III-B Ansatz

For VQE we use Real Amplitudes Ansatz with a single repetition provided by
Qiskit. Real Amplitudes ansatz has alternate layers of RY gates and CX gates.
Default full entanglement map was used which applies the CX gate to all the
qubit pairs.

For QAOA and CVaR QAOA, the ansatz is as discussed in II. We experiment with
different numbers of repetitions. Theoretically, as the number of repetitions
increases, the solution should converge to the exact solution. However, it was
found that the best results were obtained for p=2 repetitions. For all cases
except OQC Lucy(see Section V) p=2 was chosen.

Figure 3: QAOA performance vs Ansatz repetitions. p represents the depth of
the QAOA circuit.

###  III-C Optimizer

We use the COBYLA optimizer with 50 iterations. COBYLA is a classical
gradient-free optimizer and calls the cost function only once per iteration.
This is important from a cost point of view. Other optimizers such as gradient
descent or SPSA could have been possibly chosen; however, in the interest of
uniformity across different experiments and relatively erratic behaviour of
other optimizers COBYLA was chosen.

###  III-D Error Mitigation

NISQ devices are typically noisy and there exist several methods to
mitigate/correct for these effects. There are two current approaches, the
first being error correction, which involves methods to remove error during
computation. On the other hand the error mitigation methods try to compensate
for the errors post measurement. We have used the error mitigation technique
provided by qiskit simulating experiments on IBM Perth and FakeBrooklyn -
noisy backend. A calibration matrix is calculated by fitting the measured
state and the idea state data. The calibration matrix is then used to correct
the states used in optimization.

##  IV Results

###  IV-A Ground State Energies on Hardware and Simulators

The tables I through V summarize the ground state energies(in Hartrees) as
achieved on different quantum hardware and simulators across the four
algorithm types. Typically CVaR Expectation by definition is less than the
expectation value. $\alpha$=0.01 was set for all cases except for IBM
FakeBrooklyn, where $\alpha$=0.1 was chosen to compensate for fewer shots. The
reduced number of shots for IBM FakeBrooklyn stems from memory and run-time
constraints.

CVaR VQE and CVaR QAOA seem to converge to the same value accross hardwares
and simulators for simulations without the side chain as can be seen in Tables
1, II, and III. Error mitigation has no effect in this case. For QAOA and VQE
significant variations were noted. VQE in general performed better than QAOA,
with VQE showing improvements when error mitigation was applied.

For simulations with the sidechain, VQE and CVaR VQE seem to outperform their
QAOA counterparts as is evident in Tables IV and V.

IBM PERTH
---
Algorithm |  Ground State Energy |  Error Mitigated Energy
VQE |  391.166 |  86.89
CVaR VQE |  -1.019 |  -1.019
QAOA |  986.73 |  832.81
CVaR QAOA |  -1.019 |  -1.019
TABLE I: Table 1 shows the results of experiments ran on IBM Perth Hardware.
The main chain is YPYFIP and there is no side chain used for this experiment.
It took 6 qubits and the number of shots used are 8192. Fake Brooklyn
---
Algorithm |  Ground State Energy |  Error Mitigated Energy
VQE |  114.91 |  84.22
CVaR VQE |  -1.019 |  -1.019
QAOA |  623.62 |  666.46
CVaR QAOA |  -1.019 |  -1.019
TABLE II: Table 2 shows the results of experiments ran on Fake Brooklyn
Simulator. The main chain is YPYFIP and there is no side chain used. It took 6
qubits and the number of shots used are 128.  Lucy
---
Algorithm |  Ground State Energy
VQE |  768.7303
CVaR VQE |  -1.019
QAOA |  746.8159
CVaR QAOA |  -1.019
TABLE III: Table 3 shows the results of experiments ran on Lucy. There is no
error mitigation technique used in running this experiment. The main chain is
YPYFIP and the there is no side chain used. It took 6 qubits and the number of
shots used are 8192.  Fake Brooklyn
---
Algorithm |  Ground State Energy
VQE |  1601.5051
CVaR VQE |  17.1838
QAOA |  2123.2193
CVaR QAOA |  27.562968
TABLE IV: Table 4 shows the results of experiments run on Fake Brooklyn
Simulator. There is no error mitigation technique used in running this
experiment. The main chain is YPYFIP and the side chain is IPFY. It took 17
qubits and the number of shots used are 128.  SV1 Simulator
---
Algorithm |  Ground State Energy
VQE |  226.27
CVaR VQE |  -1.019
QAOA |  1741.854
CVaR QAOA |  8.33775
TABLE V: Table 5 shows the results of experiments run on SV1 Simulator. There
is no error mitigation technique used in running this experiment. The main
chain is YPYFIP, and the side chain is IPFY. It took 17 qubits and the number
of shots used were 8192.  Figure 4: Sample protein structure with sidechain

###  IV-B Distribution of Algorithm Predictions

We study the repeatability of each of the four algorithms. Thus to quantify
the variation in results, we ran experiments on a noiseless simulator with the
protein chain without any side chains. We characterize the distribution of
predicted ground state energy(or CVaR expectation in the case of CVaR
variants) and the predicted protein structure across multiple trials to
quantify the efficacy of the algorithms. The predicted structures show
considerable variability as is evident in the figures 5 6 7 8. The results
thus, not only help in quantifying the efficacy of the method in its ultimate
goal which is to predict the 3-D structure but also provide a more robust way
of comparing the predicted structures, as will be discussed in the next
section. VQE and CVaR QAOA were statistically the most consistent in terms of
the returning the same structure. However, both of them predicted physically
impossible structures for a significant number of trials.

|
---|---
|

(a) Distribution of predicted eignevalue

(b) Distribution of predicted peptide folding, the x-axis labels denote the
turns as defined in II

(c) Corresponding predicted structure

(d) Second ranked prediction

Figure 5: CVaR QAOA: The histogram(a) summarises the predicted results as
obtained on the QASM simulator for CVaR QAOA ($\alpha$=0.01). The prediction
was made 512 times; the eigenvalue converged to the same value 98.8% of the
time. However, the most frequent prediction had only a frequency of 8.8%.
Moreover, Figure (d) indicates that the chain folds back onto itself leading
to overlap, something that is forbidden.

|
---|---
|

(a) Distribution of predicted eigenvalue

(b) Distribution of predicted protein folding, the x-axis labels denote the
turns as defined in II

(c) First ranked prediction

(d) Second ranked prediction

Figure 6: QAOA: The histogram(a) summarises the predicted results as obtained
on the QASM simulator for QAOA. The prediction was made 512 times, the
algorithm converged to the lowest energy state less than 2% of the time. The
most frequent prediction had only a probability of 8.4%. The top two most
frequent predictions do not show signs of overlap.

|
---|---
|

(a) Distribution of predicted eignevalue

(b) Distribution of predicted protein folding, the x-axis labels denote the
turns as defined in II

(c) First ranked prediction

(d) Second ranked prediction

Figure 7: VQE: The histogram(a) summarises the predicted results as obtained
on the QASM simulator for VQE. The prediction was made 512 times, the
algorithm converged to the lowest energy state 30% of the time. The most
frequent prediction had a frequency of 16%. Despite promising numbers, Figure
(d) indicates that the chain folds back onto itself leading to overlap,
something that is forbidden.

|
---|---
|

(a) Distribution of predicted eignevalue

(b) Distribution of predicted protein folding, the x-axis labels denote the
turns as defined in II

(c) First ranked prediction

(d) Second ranked prediction

Figure 8: CVaR VQE: The histogram(a) summarises the predicted results as
obtained on the QASM simulator for VQE. The prediction was made 512 times, the
algorithm converged to the lowest energy state 99% of the time. The most
frequent prediction had a frequency of only 5%. Moreover, numbers, Figure (d)
indicates that the chain folds back onto itself leading to overlap, something
that is forbidden.

###  IV-C RMSD Calculation

Root Mean Square Deviation(RMSD) is the most commonly used distance-based
measurement tool to compare the protein structure similarity between two
atomic coordinates. RMSD values are calculated by the formula

$\sqrt{1/n\sum_{i=1}^{n}d_{i}^{2}}\ $ (5)

It is to be noted that the any RMSD value less that 2 Å is a good value,
nevertheless, the lesser the value, the higher the accuracy. We compare the
top predictions by QAOA class of algorithms and VQE class of algorithms by
computing the RMSD [12]. The RMSD value turned out to be 0.6357Å between the
QAOA class and VQE class (Fig 9).

|
---|---
|

(a) Protein Structure with QAOA

(b) Protein Structure with VQE

(c) Comparing the structures of QAOA and VQE

(d) Comparing AlphaFold Predictions(Red) with VQE top prediction(Yellow)

Figure 9: Comparing most probable structures predicted by QAOA and VQE

Further we compare the predictions made by AlphaFold [1] with our top
predictions using the same method. In order to make the two predictions
comparable, alpha fold prediction (Figure 11) was reduced to a bead like model
by considering the peptide as a chain of geometric centers of each amino acid.
The comparison is illustrated in Figure 9 (d). RMSD value of 4.795Å was
obtained.

There is a noticeable difference in the scale of the two images. We conjecture
this is because of the lattice parameter inbuilt in Qiskit Nature.

###  IV-D Effect of penalty terms

The effects of the three penalty terms described in Section II are studied.
Each of the three terms is varied on a logarithmic scale keeping other terms
constant and the probability of the most frequent structure is noted (Fig 10).

Figure 10: % probability of the most frequent structure plotted against the
penalty term magnitude as defined in Section II

##  V Discussion

The cost of executing the VQE algorithm scales with two factors; the number of
quantum tasks and the number of shots. The former being the costlier of the
two. These two correspond to the number of circuits created and the number of
shots executed on each circuit. If the Hamiltonian terms commute it is
possible to measure them simultaneously thus reducing the number of calls made
to the quantum computer, and reducing the major cost component. However, if
the terms do not commute, one has to measure each hamiltonian term separately
which scales as $O(N^{4})$. In our case, all Hamiltonian terms were composed
of Pauli Z matrices. Our cost thus only scaled with the number of shots.

Another important factor that determines the capability of the hardware to
simulate is the pulse duration. OQC Lucy allows a maximum pulse duration for
120 $\mu$s thus limiting the circuit depth to 1 at 8192 shots for QAOA.

The number of iterations is also important to consider. Optimizers such as
gradient descent or SPSA might make more than one call to the quantum computer
per iteration further increasing the cost multifold(typically three-fold for
gradient-based optimizers as they require the calculation of gradient by
parameter shift.)

To reduce the computational costs, the proteins qubit operator has been
reduced to the minimum amount of qubits required by omitting identity
operations. For the protein with the main chain, that has been used in the
paper, the problem is mapped to 120 qubits but it has been optimized to 17
qubits. Similarly, it has been optimized to 6 qubits in the case without a
chain.

The structures obtained using the Variational Algorithms differ considerably
from similar techniques. To illustrate, we include the structures predicted by
alpha fold and alpha fold 2.

Figure 11: Structure Predicted by AlphaFold Figure 12: Structure Predicted by
AlphaFold 2

The reliability of the model is quite low, a probable cause might be that the
energy difference between different structures is not substantial. Further,
since only a limited number of interactions were considered, this leads to
inaccuracies. A pressing problem with the current model is that chain folds
onto itself on multiple occasions.

##  VI Conclusion

In the paper, we analyzed the performance of the protein folding model
provided by Qiskit as a part of the Qiskit Nature library and compare the
results against the Alphafold model and find that there are significant
deviations which require that the Qiskit Nature model be improved. The
algorithm is one of a kind, especially in the paradigm of gate-based quantum
computing.

We compare VQE and QAOA on various hardware and simulators, with and without
error mitigation. We conclude that VQE class of algorithms performed better
than QAOA for simulations with and without a side chain. We also compare the
structures that were predicted by the algorithms and go on to examine the
repeatability of the model we find that several unphysical predictions were
being made frequently. Further, we also examined the effect of penalty
parameters on the statistical stability of the model and find a weak trend in
the same.

We conclude that there is an extensive scope of improvement on various
metrics, however, we also acknowledge that the current limitations on the
hardware make it difficult to achieve. Continuous work in this direction is a
positive indication. A more robust encoding method that prevents folding of
the chain back onto itself must be developed among other things.

##  VII Future Research

Future developments in this field would closely follow the availability of
hardware capabilities. The primary goal for future works is expected to be
modelling the problem in a way that allows more degrees of freedom and
prohibits unphysical states. Further, the reliability of the model must be
improved to allow for commercial and practical use.

One avenue that we identify is the use of other emerging techniques such as
Density Matrix Embedding Theory to better characterize the energies of the
proteins.

## Acknowledgments

The authors acknowledge using Google Colab Pro and AWS Braket for carrying out
the experiments and the libraries Pennylane and Qiskit Nature from IBM
Quantum.

## References

  * [1] J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Potapenko, A. Bridgland, C. Meyer, S. A. A. Kohl, A. J. Ballard, A. Cowie, B. Romera- Paredes, S. Nikolov, R. Jain, J. Adler, T. Back, S. Petersen, D. Reiman, E. Clancy, M. Zielinski, M. Steinegger, M. Pacholska, T. Berghammer, S. Bodenstein, D. Silver, O. Vinyals, A. W. Senior, K. Kavukcuoglu, P. Kohli, and D. Hassabis, ”Highly accurate protein structure prediction with AlphaFold,” Nature, vol. 596, no. 7873, pp. 583–589, 2021.
  * [2] A. Robert, P. K. Barkoutsos, S. Woerner, and I. Tavernelli, “Resource-efficient quantum algorithm for protein folding,” npj Quantum Information, vol. 7, no. 1, feb 2021. [Online]. Available: https://doi.org/10.1038
  * [3] K. A. Dill, S. B. Ozkan, M. S. Shell, and T. R. Weikl, “The protein folding problem,” Annual review of biophysics, vol. 37, p. 289, 2008.
  * [4] M. S. R. K. J. E. Maria Kieferova, Paul K. Faehrmann, “Randomizing multi-product formulas for improved hamiltonian simulation,” arXiv preprint arXiv:2101.07808, 2021.
  * [5] S. Miyazawa and R. L. Jernigan, “Residue–residue potentials with a favorable contact pair term and an unfavorable high packing density term, for simulation and threading,” Journal of molecular biology, vol. 256, no. 3, pp. 623–644, 1996.
  * [6] P. K. Barkoutsos, G. Nannicini, A. Robert, I. Tavernelli, and S. Woerner, “Improving variational quantum optimization using CVaR,” Quantum, vol. 4, p. 256, apr 2020. [Online]. Available: https://doi.org/10.22331
  * [7] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, “A variational eigenvalue solver on a photonic quantum processor,” Nature communications, vol. 5, no. 1, pp. 1–7, 2014.
  * [8] E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm,” arXiv preprint arXiv:1411.4028, 2014.
  * [9] P. K. Barkoutsos, G. Nannicini, A. Robert, I. Tavernelli, and S. Woerner, “Improving variational quantum optimization using cvar,” Quantum, vol. 4, p. 256, 2020.
  * [10] L. J. Gutierrez, E. Angelina, A. Gyebrovszki, L. F, N. Peruchena, H. A. Baldoni, B. Penke, and R. D. Enriz, “New small-size peptides modulators of the exosite of bace1 obtained from a structure-based design,” Journal of Biomolecular Structure and Dynamics, vol. 35, no. 2, pp. 413–426, 2017.
  * [11] K. Wright, K. M. Beck, S. Debnath, J. Amini, Y. Nam, N. Grzesiak, J.- S. Chen, N. Pisenti, M. Chmielewski, C. Collins et al., “Benchmarking an 11-qubit quantum computer,” Nature communications, vol. 10, no. 1, pp. 1–6, 2019.
  * [12] R. A. Irina Kufareva, “Methods of protein structure comparison,” Meth- ods in molecular biology (Clifton, N.J.), vol. 857, pp. 231–57, 2015.