# Upper bounds on the sizes of variable strength covering arrays using the
Lovász local lemma

Lucia Moura111University of Ottawa, Ottawa, Canada<EMAIL_ADDRESS>Sebastian Raaphorst222Gemini Observatory, La Serena, Chile,
<EMAIL_ADDRESS>and Brett Stevens333Carleton University, Ottawa, Canada,
<EMAIL_ADDRESS>
(January 15, 2019)

###### Abstract

Covering arrays are generalizations of orthogonal arrays that have been widely
studied and are used in software testing. The probabilistic method has been
employed to derive upper bounds on the sizes of minimum covering arrays and
give asymptotic upper bounds that are logarithmic on the number of columns of
the array. This corresponds to test suites with a desired level of coverage of
the parameter space where we guarantee the number of test cases is logarithmic
on the number of parameters of the system. In this paper, we study variable
strength covering arrays, a generalization of covering arrays that uses a
hypergraph to specify the sets of columns where coverage is required;
(standard) covering arrays is the special case where coverage is required for
all sets of columns of a fixed size $t$, its strength. We use the
probabilistic method to obtain upper bounds on the number of rows of a
variable strength covering array, given in terms of parameters of the
hypergraph. We then compare this upper bound with another one given by a
density-based greedy algorithm on different types of hypergraph such as
$t$-designs, cyclic consecutive hypergraphs, planar triangulation hypergraphs,
and a more specific hypergraph given by a clique of higher strength on top of
a “base strength”. The conclusions are dependent on the class of hypergraph,
and we discuss specific characteristics of the hypergraphs which are more
amenable to using different versions of the Lovász local lemma.

Keywords: covering arrays, variable strength, Lovász local lemma, greedy
algorithm, derandomization.

##  1 Introduction

Covering arrays are well studied combinatorial designs [colbourn04,
crchandbook] in part because of their utility in software and network testing
[cohen96, dalal98, kuhn04]. For more information about covering arrays,
including combinatorial constructions and overview on algorithms and
applications, see the survey by Colbourn [colbourn04]. In this paper, we focus
on a recent covering array generalization called variable strength covering
array (VCA). We begin defining VCAs and indicate how covering arrays are a
special case. For more information on VCAs, see [raaphorstphd, vca].

######  Definition 1.1.

Let $H=(V,E)$ be a hypergraph and let $k=|V|$. A variable-strength covering
array, denoted ${\rm VCA}(n;H,v)$, is an $n\times k$ array filled from
$\mathbb{Z}_{v}$ such that for each $e=\\{v_{0},\ldots,v_{t-1}\\}\in E$, the
$n\times t$ subarray of columns indexed by the elements of $e$ is covered,
that is, it has every possible $t$-tuple in $\mathbb{Z}_{v}$ as a row at least
once. The variable-strength covering array number, written ${\rm VCAN}(H,v)$,
is the smallest $n$ such that a ${\rm VCA}(n;H,v)$ exists.

Consider the complete $t$-uniform hypergraph on $k$ vertices, denoted
$K_{k}^{(t)}$, that is, the hypergraph where $|V|=k$ and $E$ is the set of all
$t$-subsets of $V$. A covering array of strength $t$, denoted
$\mathrm{CA}(N;t,k,v)$, is precisely a ${\rm VCA}(N;K_{k}^{(t)},v)$; the
covering array number is denoted $\mathrm{CAN}(t,k,v)$. In this article, we
use the Lovász local lemma to give an upper bound on ${\rm VCAN}(H,v)$.

######  Theorem 1.2 (Lovász local lemma - symmetric case [lovasz75,
spencer77]).

Let $\mathcal{A}=\\{A_{0},\ldots,A_{m-1}\\}$ be a finite set of events in a
probability space $\Omega$ such that $\Pr(A_{i})\leq p<1$, and each event is
independent of all but at most $d$ of the other events. If $ep(d+1)\leq 1$,
where $e$ is the base of the natural logarithm, then the probability that none
of the events occur is nonzero.

For standard covering arrays, the local lemma has been used [godbole96] to
prove

$\mathrm{CAN}(t,k,v)\leq\frac{(t-1)\log
k}{\log\frac{v^{t}}{v^{t}-1}}(1+o(1)).$

More recently, this local lemma technique and the deterministic analogue,
often called entropy compression [moser2010], has improved the coefficient of
the leading term in this bound for all $v$ when $t=3$ and for $v\leq 10$ when
$4\leq t\leq 6$, by using columns that have balanced numbers of symbols
[MR3328867, francetic2016]. Sarkar and Colbourn match and extend these
improvements by exploiting permutation groups with the use of the local lemma
[sarkar_upper_2017]. They also combine the local lemma, permutation groups,
graph colouring techniques and a density-based greedy approach into two-stage
methods which further improve covering array upper bounds [colbourn_two-
stage_nodate_alt]. In [colbourn_asymptotic_nodate] together with Lanus they
use the local lemma to construct covering perfect hash families which are then
used to construct covering arrays; covering perfect hash families are smaller
than covering arrays and thus can be more efficiently generated. These
constructions are related to covering arrays constructed from linear feedback
shift register sequences [raaphorst2012, MR3433921].

Probabilistic methods have been previously used in the context of variable
strength covering arrays. In [godbole2011, godbole2010cca] probabilistic
methods, but not explicitly the local lemma, are used to derive bounds on
variable strength covering arrays for consecutive hypergraphs which are
similar to the cyclic consecutive hypergraphs we discuss in Section 3 but
without the edges wrapping from the end of the vertex set to the beginning.
Sarkar et al. [MR3808633] use the local lemma to construct partial covering
arrays which cover at least $1-\epsilon$ of the possible $t$-sets of columns.
The main difference between these and variable strength covering arrays is
that variable strength covering arrays specify exactly which $t$-sets of
columns will be covered and partial covering arrays only specify that a
certain proportion of $t$-sets must be covered, without constraining which
ones they are.

###  1.1 Main Result

For a hypergraph $H$, let the rank of $H$, denoted $rank(H)$, be the largest
cardinality of an edge in $H$.

######  Theorem 1.3.

Let $H=(V,E)$ be a hypergraph with $\text{rank}(H)=t\geq 1$, and let $d$ be an
integer such that no edge of $H$ intersects more than $d$ other edges of $H$.
Then, for any $v\geq 2$, we have:

$\displaystyle{\rm VCAN}(H,v)$ $\displaystyle\leq$
$\displaystyle\left\lceil\frac{\ln(d+1)+t\ln
v+1}{\ln\frac{v^{t}}{v^{t}-1}}\right\rceil$ (1.1) $\displaystyle\leq$
$\displaystyle\left\lceil v^{t}\left(\ln(d+1)+t\ln v+1\right)\right\rceil.$
(1.2)

###### Proof.

Let $k=|V|$. Consider a randomly generated $N\times k$ array $M$ with entries
chosen from $\mathbb{Z}_{v}$ with uniform probability. For each edge $e\in E$,
write $s=|e|$, and associate an event $A_{e}$ that the $N\times s$ subarray of
$M$ formed by the columns corresponding to $e$ is missing one or more of the
$v^{s}$ $s$-tuples of $\mathbb{Z}_{v}^{s}$ as a row. Define:

$p=v^{t}\left(\frac{v^{t}-1}{v^{t}}\right)^{N}.$

Since $s\leq t$, we have that

$\Pr(A_{e})\leq v^{s}\left(\frac{v^{s}-1}{v^{s}}\right)^{N}\leq
v^{t}\left(\frac{v^{t}-1}{v^{t}}\right)^{N}=p.$

We apply the symmetric Local Lemma as in Theorem 1.2, which states that if
$ep(d+1)\leq 1$, the probability that none of the events occur is positive,
and hence there is some $N\times k$ array that is a ${\rm VCA}(N;H,v)$. This
happens when:

$N\geq\frac{\ln(d+1)+t\ln v+1}{\ln\frac{v^{t}}{v^{t}-1}}.$

Thus, we have:

${\rm VCAN}(H,v)\leq\frac{\ln(d+1)+t\ln v+1}{\ln\frac{v^{t}}{v^{t}-1}}.$

If we use the approximation

$\left(\ln\frac{m}{m-1}\right)^{-1}<m,$

for $m>1$ from the Taylor series expansion, we can rewrite the equation as
follows:

${\rm VCAN}(H,v)\leq v^{t}\left(\ln(d+1)+t\ln v+1\right).$

∎

Theorem 1.3 can be applied to any hypergraph $H$ and thus gives a very general
existence result for VCAs with arbitrary parameters. The rest of the paper
focuses on comparing the upper bound given by Theorem 1.3 with a constructive
upper bound for VCAs obtained by a density-based greedy algorithm called
VarDens introduced in [raaphorstphd, vcadbga] and given next. This method is a
generalization of the density method of Bryce and Colbourn [bryce09] for
dealing with variable strength.

######  Theorem 1.4.

[raaphorstphd, vcadbga] Let $H$ be a hypergraph such that $rank(H)=t$ and
$e=|E(H)|$. Algorithm VarDens returns a ${\rm VCA}(N;H,v)$ where $N$ satisfies

$\displaystyle{\rm VCAN}(H,v)\leq N$ $\displaystyle\leq$
$\displaystyle\left\lceil\frac{\ln e+t\ln
v}{\ln\frac{v^{t}}{v^{t}-1}}\right\rceil$ (1.3) $\displaystyle\leq$
$\displaystyle\left\lceil v^{t}(\ln e+t\ln v)\right\rceil.$ (1.4)

For the remainder of the paper we compare the bounds given by the
probabilistic method and by VarDens algorithm. Let $N_{prob}(H,v)$ be the
upper bound given by the probabilistic method (the right hand side of equation
(1.1) in Theorem 1.3), and let $N_{dens}(H,v)$ be the upper bound given by
VarDens algorithm (the right hand side of equation (1.3) in Theorem 1.4). That
is, we denote

$\displaystyle N_{prob}(H,v)$ $\displaystyle:=$
$\displaystyle\left\lceil\frac{\ln(d+1)+t\ln
v+1}{\ln\frac{v^{t}}{v^{t}-1}}\right\rceil,$ $\displaystyle N_{dens}(H,v)$
$\displaystyle:=$ $\displaystyle\left\lceil\frac{\ln e+t\ln
v}{\ln\frac{v^{t}}{v^{t}-1}}\right\rceil.$

We note that if absolutely nothing is known about the hypergraph $H$ except
the number of edges $e$ and the rank $t$, then we can substitute $d\leq e-1$,
into $N_{prob}(H,v)$ obtaining an upper bound quite close to $N_{dens}(H,v)$,
specially as we fix $v$ and $t$ and let $e$ grow.

In the rest of the article, we show that when we know better estimates on $d$,
$N_{prob}(H,v)$ outperforms $N_{dens}(H,v)$. In Section 2, we look at
hypergraphs that are combinatorial designs. In Section 3, we study two
families of hypergraphs (cyclic consecutive hypergraphs and planar
triangulations) where we know the exact ${\rm VCAN}$ to assess how well each
of these methods perform. In Section 4, we consider other versions of the
Lovász local lemma, namely the asymmetric and the general cases, and exemplify
challenges and benefits to using them. An extended abstract containing the
main results and discussions in Sections 1-3 appeared in [cai-ea].

##  2 VCAs over designs

Combinatorial designs can be used to obtain hypergraphs that have a great deal
of regularity, which can be exploited to determine the number of dependent
events.

######  Definition 2.1.

An $s$-$(k,t,\lambda)$ design is a collection $\mathcal{B}$ of $t$-subsets
(called blocks) of a $k$-set $V$ with the property that any $s$-subset of
points from $V$ appear in exactly $\lambda$ blocks of $\mathcal{B}$.

For a fixed $s$, $t$ and $\lambda$, these designs are known to exist for all
sufficiently large $k$ that satisfy the necessary conditions [1401.3665].

We begin by counting the number of blocks in a $s$-$(k,t,\lambda)$ design that
intersect a fixed block.

######  Lemma 2.2.

Let $\mathcal{B}$ be an $s$-$(k,t,\lambda)$ design. Then for any
$B\in\mathcal{B}$:

$d=|\\{B^{\prime}\in\mathcal{B}\setminus\\{B\\}:B\cap
B^{\prime}\neq\emptyset\\}|\leq\sum_{i=1}^{2\lfloor\frac{s-1}{2}\rfloor+1}(-1)^{i+1}{t\choose
i}\left(\frac{\lambda{k-i\choose s-i}}{{t-i\choose s-i}}-1\right).$

When $s=t-1$ and the design has no repeated blocks,

$d=\sum_{i=1}^{t-1}(-1)^{i+1}{t\choose
i}\left(\frac{\lambda{k-i\choose(t-1)-i}}{t-i}-1\right).$

Furthermore, for any $1\leq m\leq s$, $m$ odd, we can truncate the summation
after $m$ terms to derive an upper bound on the summation.

###### Proof.

Let $K$ be the point set of the design. To count the number of blocks that
intersect $B$ in at least one point we use the inclusion-exclusion principle.
Let $\lambda_{U}$ be the number of blocks that contain a set $U\subseteq K$.
The inclusion-exclusion principle states that the number of blocks that
intersect $B\in\mathcal{B}$ in at least one point is precisely

$d=\sum_{i=1}^{t}(-1)^{i+1}\sum_{\scriptsize\begin{array}[]{c}U\subset B\\\
|U|=i\end{array}}\left(\lambda_{U}-1\right)$

where the 1 is subtracted to not count $B$ itself. Whenever the outer
summation is truncated ending with a positive term an upper bound is achieved.

For $|U|\leq s$, $\lambda_{U}$ can be computed from the parameters of the
design. Each of the ${t\choose i}$ $i$-sets $U\subseteq B$, occurs in exactly
${k-i\choose s-i}$ $s$-sets each of which appears $\lambda$ times amongst the
blocks of $\mathcal{B}$. Each block of $\mathcal{B}$ which contains $U$,
contains exactly ${t-i\choose s-i}$ $s$-sets containing $U$, and hence
$\lambda_{U}=\lambda\frac{{k-i\choose s-i}}{{t-i\choose s-i}}$ blocks of
$\mathcal{B}$. For $s+1\leq|U|\leq t$, $\lambda_{U}$ cannot be derived from
just the parameters of the design. If there are repeated blocks then
$\lambda_{B}$ itself may be more than one. Thus when $s<t-1$, or repeated
blocks are permitted, the value of $d$ will depend on the particular design
and not just on the parameters. So we calculate an upper bound on $d$, by
truncating the inclusion-exclusion after a positive term. We stop at the
largest odd integer no larger than $s$, that is
$2\lfloor\frac{s-1}{2}\rfloor+1$.

When $s=t-1$ and there are no repeated blocks, every term in the summation can
be computed from the parameters of the design and the computation is exact. ∎

######  Theorem 2.3 (Bound Comparison for ${\rm VCA}$ over $s$-$(k,t,\lambda)$
designs).

Let $\mathcal{B}$ be an $s$-$(k,t,\lambda)$ design. Then,

$N_{prob}(\mathcal{B},v)=\left\lceil\frac{\ln\left(1+\sum_{i=1}^{2\lfloor\frac{s-1}{2}\rfloor+1}(-1)^{i+1}{t\choose
i}\left(\frac{\lambda{k-i\choose s-i}}{{t-i\choose s-i}}-1\right)\right)+t\ln
v+1}{\ln\frac{v^{t}}{v^{t}-1}}\right\rceil,$

and $N_{dens}(\mathcal{B},v)=\left\lceil\frac{\ln b+t\ln
v}{\ln\frac{v^{t}}{v^{t}-1}}\right\rceil.$ Furthermore, for $s$, $t$,
$\lambda$, and $v$ fixed, as $k\rightarrow\infty$, we get

$N_{prob}(\mathcal{B},v)=(s-1)v^{t}\ln k+O(1)\;\;\mbox{ and
}\;\;N_{dens}(\mathcal{B},v)=sv^{t}\ln k+O(1).$

###### Proof.

The first equality comes from Lemma 2.2 to bound $d$ in Theorem 1.3, the
second equality comes from Theorem 1.4 and the asymptotic results follow
easily from these using the Taylor series of the denominator.∎

We conclude that for designs, $N_{prob}(\mathcal{B},v)$ is asymptotically
better than $N_{dens}(\mathcal{B},v)$ as it reduces the coefficient for the
leading term $\ln k$ by $v^{t}$.

When $s=t-1$ and the design has no repeated blocks we can take advantage of
the equality in Lemma 2.2. Table 1 gives the bounds from Theorem 2.3 for
$(t-1)$-$(k,t,1)$ designs for $3\leq t\leq 6$. When $t\leq 6$, explicit
designs are known [crchandbook]. Although for small $k$, $N_{dens}$ may
outperform $N_{prob}$, as $k$ grows even modestly, $N_{prob}$ becomes better.

Table 1: VCA upper bounds on $(t-1)$-$(k,t,1)$ designs.
${\tiny\begin{array}[]{|c|l|l|l|l|}\hline\cr&$2$-$(k,3,1)$\mbox{
designs}&$3$-$(k,4,1)$\mbox{ designs}&$4$-$(k,5,1)$\mbox{
designs}&$5$-$(k,6,1)$\mbox{ designs}\\\ \hline\cr\mbox{
$N_{prob}$}&\frac{\ln\left(\frac{3}{2}k-\frac{7}{2}\right)+3\ln
v+1}{\ln\frac{v^{3}}{v^{3}-1}}&\frac{\ln\left(\frac{2}{3}k^{2}-\frac{11}{3}k+9\right)+4\ln
v+1}{\ln\frac{v^{4}}{v^{4}-1}}&\frac{\ln\left(\frac{5}{24}k^{3}-\frac{35}{12}k^{2}+\frac{125}{8}k-\frac{121}{4}\right)+5\ln
v+1}{\ln\frac{v^{5}}{v^{5}-1}}&\frac{\ln\left(\frac{1}{20}k^{4}-\frac{9}{8}k^{3}+\frac{257}{24}k^{2}-\frac{595}{12}k+\frac{456}{5}\right)+6\ln
v+1}{\ln\frac{v^{6}}{v^{6}-1}}\\\
\hline\cr\mbox{$N_{dens}$}&\frac{\ln\left(\frac{1}{6}k^{2}-\frac{1}{6}k\right)+3\ln
v}{\ln\frac{v^{3}}{v^{3}-1}}&\frac{\ln\left(\frac{1}{24}k^{3}-\frac{1}{8}k^{2}+\frac{1}{12}k\right)+4\ln
v}{\ln\frac{v^{4}}{v^{4}-1}}&\frac{\ln\left(\frac{1}{120}k^{4}-\frac{1}{20}k^{3}+\frac{11}{120}k^{2}-\frac{1}{20}k\right)+5\ln
v}{\ln\frac{v^{5}}{v^{5}-1}}&\frac{\ln\left(\frac{1}{720}k^{5}-\frac{1}{72}k^{4}+\frac{7}{144}k^{3}-\frac{5}{72}k^{2}+\frac{1}{30}\right)+6\ln
v}{\ln\frac{v^{6}}{v^{6}-1}}\\\ \hline\cr\end{array}}$

##  3 Ability of bounds to capture global properties of $H$

In this section, we compare the two bounds for two families of hypergraphs for
which VCANs are almost completely known and do not grow with $k$ [vca]. The
cyclic consecutive hypergraph is $H_{c}^{k,t}=(V,E)$ with
$V=\\{0,\ldots,k-1\\}$ and
$E=\\{\\{i,i+1\bmod{k},\ldots,i+t-1\bmod{k}\\}:0\leq i\leq k-1\\}$; for many
cyclic consecutive hypergraphs ${\rm VCAN}(H_{c}^{k,t},v)$ equals $v^{t}$, and
in all cases we know that ${\rm VCAN}(H_{c}^{k,t},v)\leq d_{t,v}$ for some
$d_{t,v}$ that does not grow with $k$ [vca, Theorem 3.8].

In $H_{c}^{k,t}$, each edge intersects exactly $2t-2$ other edges. Theorem 1.3
gives

$N_{prob}(H_{c}^{k,t},v)=\left\lceil\frac{\ln(2t-1)+t\ln
v+1}{\ln\frac{v^{t}}{v^{t}-1}}\right\rceil.$ (3.1)

In agreement with $d_{t,v}$, $N_{prob}$ is independent of $k$ and thus for
fixed $v$ and $t$, it is $O(1)$. On the other hand, the bound from algorithm
VarDens is

$N_{dens}(H_{c}^{k,t},v)=\left\lceil\frac{\ln k+t\ln
v}{\ln\frac{v^{t}}{v^{t}-1}}\right\rceil,$

which does grow with $k$. Thus we have an example of a family of hypergraphs
where $N_{prob}$ is substantially better than $N_{dens}$. The homomorphism
construction (see [raaphorstphd, Chapter 3]) gives
$VCAN(H_{c}^{k,t},v)\leq{\rm VCAN}(H_{c}^{k^{\prime},t},v)$, where
$k^{\prime}\leq 2t-1$ and the $\ln(2t-1)$ term in Equation 3.1 in place of the
usual $\ln k$ term suggests that $N_{prob}$ “recognizes” this homomorphism
while $N_{dens}$ does not. However, experiments show that running VarDens
algorithm for $H_{c}^{k,t}$ does seem to generate arrays where the array size
is independent of $k$ [raaphorstphd].

A triangulation hypergraph of the sphere, $T=(V,E)$ is a rank-3 hypergraph
which corresponds to a planar graph all of whose faces are triangles; the
rank-3 hyperedges are precisely the faces of the planar embedding. From
colourings and homomorphisms we know that ${\rm VCAN}(T,v)=v^{3}$ [vca]. If
$T$ has $k$ vertices then it has $2k-4$ edges (triangles) and a triangle may
intersect up to $3(\Delta-2)$ other triangles where $\Delta$ is the maximum
degree of the planar graph. Thus, $d\leq\min\\{3(\Delta-2),2k-5\\}$ and

$N_{prob}(T,v)\leq\left\lceil\frac{\ln(\min\\{3(\Delta-2)+1,2k-4\\})+3\ln
v+1}{\ln\frac{v^{3}}{v^{3}-1}}\right\rceil,$

while VarDens gives

$N_{dens}(T,v)=\left\lceil\frac{\ln(2k-4)+3\ln
v}{\ln\frac{v^{3}}{v^{3}-1}}\right\rceil.$

Thus, $N_{prob}$ only performs better than $N_{dens}$ when the maximum degree
of the triangulation grows more slowly than the number of vertices. In
addition, unless $\Delta(T)$ is bounded as $k$ grows, $N_{prob}$ grows with
$k$, even though ${\rm VCAN}(T,v)=v^{3}$. Hence, triangulation hypergraphs
represent a case where, unlike cyclic consecutive hypergraphs, our application
of the local lemma is unable to capture the behaviour of the homomorphism
$H\rightarrow K_{4}^{(3)}$. The bound $N_{dens}$ also grows with $k$ and, in
contrast to the previous class of hypergraphs, experiments running VarDens on
random triangulation hypergraphs suggest that the size of the arrays produced
indeed increases as some function of $k$, see Table 2. For more details about
generating the random triangulation hypergraphs see [raaphorstphd].

| $k=4$ | $k=11$ | $k=100$
---|---|---|---
$v$ | $N_{m}$ | $N_{M}$ | $N_{a}$ | $N_{m}$ | $N_{M}$ | $N_{a}$ | $N_{m}$ | $N_{M}$ | $N_{a}$
2 | 8 | 10 | 8.37 | 8 | 11 | 9.50 | 12 | 14 | 13.13
3 | 27 | 34 | 31.17 | 31 | 36 | 34.23 | 39 | 42 | 40.37
5 | 142 | 152 | 146.30 | 150 | 157 | 154.13 | 167 | 172 | 169.60
Table 2:  Trials of VarDens on random triangulations with $k$ vertices.
$N_{m}$, $N_{M}$, and $N_{a}$ are the minimum, maximum, and mean array sizes,
respectively.

##  4 Using the general and asymmetric local lemma

In Theorem 1.3, we use the symmetric version of the local lemma to establish
the existence of covering arrays. In Section 2, we apply this theorem to
hypergraphs that are highly symmetric: the probability that any set of columns
represented by a hyperedge be covered is precisely the same as for any other
hyperedge, because all the hyperedges have the same size. Additionally, for
$(t-1)$-$(k,t,\lambda)$ designs the size of the sets of dependent events also
does not vary; even for $s$-$(k,t,\lambda)$ designs the sizes of sets of
dependent events can only vary in a fixed range. In Section 3, the hyperedges
were again of a fixed size. If the hyperedges themselves vary in size or the
size of sets of dependent events significantly vary, using the symmetric
version of the local lemma requires taking the worst probability of a set of
columns being uncovered and requires taking the largest set of dependent
events. In this section, we explore the benefit of using versions of the local
lemma that can adapt to varying sizes of hyperedges and sets of dependent
events.

The most general statement of the local lemma was given by Lovász in 1975.

######  Theorem 4.1 (Lovász local lemma - general case [lovasz75]).

Let $\mathcal{A}=\\{A_{0},\ldots,A_{m-1}\\}$ be a finite set of events in a
probability space $\Omega$. Define a function
$\Gamma\colon\mathcal{A}\rightarrow\mathcal{P}(\mathcal{A})$ such that for
$A\in\mathcal{A}$, $A$ is independent from all events in
$\mathcal{A}\setminus(\Gamma(A)\cup\\{A\\})$. If there is a map
$x\colon\mathcal{A}\rightarrow(0,1)$ such that for all $A\in\mathcal{A}$:

$\Pr(A)\leq x(A)\prod_{B\in\Gamma(A)}(1-x(B))$

then the probability that none of the events occur is nonzero, and is:

$\Pr(\overline{A_{0}}\wedge\ldots\wedge\overline{A_{m-1}})\geq\prod_{A\in\mathcal{A}}(1-x(A)).$

The Asymmetric local lemma was given by Habib in 1998.

######  Theorem 4.2 (Lovász local lemma - asymmetric case [habib1998]).

Let $\mathcal{A}=\\{A_{0},\ldots,A_{m-1}\\}$ be a finite set of events in a
probability space $\Omega$. Define a function
$\Gamma:\mathcal{A}\rightarrow\mathcal{P}(\mathcal{A})$ such that for
$A\in\mathcal{A}$, $A$ is independent from all events in
$\mathcal{A}\setminus(\Gamma(A)\cup\\{A\\})$. If, for each $0\leq i<m$, we
have that both:

  1. 1.

$\Pr(A_{i})\leq\frac{1}{8}$

  2. 2.

$\sum_{A_{j}\in\Gamma(A_{i})}\Pr(A_{j})\leq\frac{1}{4}$

then the probability that none of the events occur is positive.

For any particular hypergraph, using these more general versions of the local
lemma can require establishing complicated bounds on the probabilities and
size of sets of dependent events which are different when the hypergraph is
changed. We choose to focus on one specific hypergraph to highlight the
challenges and benefits of this approach. Let $H_{15}$ be the hypergraph with
vertex set $\\{a,b,c,d,e,f,g,h,i,j,k,0,1,2,3\\}$ and edge set containing the
four hyperedges of rank 3 from $\\{0,1,2,3\\}$ and all possible edges of rank
2 that are not contained within any hyperedge of rank 3. Cohen et al.
[cohen2003vsi] previously studied this hypergraph and other similar ones and
used simulated annealing to construct variable strength covering arrays over
them.

With $n$ rows, the probability of missing coverage on any pair of columns is

$p(A)\leq v^{2}\left(\frac{v^{2}-1}{v^{2}}\right)^{n}.$

The probability of missing coverage on one of the four sets of three columns
is

$p(A)\leq v^{3}\left(\frac{v^{3}-1}{v^{3}}\right)^{n}.$

To determine the sets of dependent events we classify the hyperedges into
three kinds. $E_{1}$ is the set of pairs from $\\{a,b,c,d,e,f,g,h,i,j,k\\}$,
$E_{2}$ is the set of pairs from
$\\{a,b,c,d,e,f,g,h,i,j,k\\}\times\\{0,1,2,3\\}$ and $E_{3}$ is the set of
triples from $\\{0,1,2,3\\}$. The number of events of each type that are
dependent on other types is summarized in Table 3. The entry in row $E_{i}$
and column $E_{j}$ is the number of $E_{j}$ events that are dependent on an
$E_{i}$ event.

| $E_{1}$ | $E_{2}$ | $E_{3}$
---|---|---|---
$E_{1}$ | $18$ | $8$ | 0
$E_{2}$ | $10$ | $13$ | $3$
$E_{3}$ | 0 | $33$ | $3$
Table 3:  Bad event dependency counts on hypergraph $H_{15}$.

We start by applying the symmetric local lemma for this hypergraph. For
$A_{1}\in E_{1}$, $|\Gamma(A_{1})|=26$; for $A_{2}\in E_{2}$,
$|\Gamma(A_{2})|=26$; and for $A_{3}\in E_{3}$, $|\Gamma(A_{3})|=36$. Thus all
bad events are avoided when:

$n\geq\frac{3\ln v+\ln 37+1}{3\ln v-\ln(v^{3}-1)}.$

These values of $n$ are given in Table 4 in column $n_{s}$ for
$v=2,\ldots,10$.

To use the general form of the local lemma, as in Theorem 4.1, we need to find
a function $x\colon\mathcal{A}\rightarrow(0,1)$ such that, for each
$A\in\mathcal{A}$:

$\Pr(A)\leq x(A)\prod_{B\in\Gamma(A)}(1-x(B)).$ (4.1)

If such a function $x$ exists, then the local lemma guarantees that the
probability that all bad events can be avoided is nonzero. We look for a
function $x$ that is fixed on each $E_{i}$. Equation (4.1) gives the following
system of inequalities:

$\displaystyle v^{2}\left(\frac{v^{2}-1}{v^{2}}\right)^{n}$ $\displaystyle\leq
x_{1}(1-x_{1})^{18}(1-x_{2})^{8},$ $\displaystyle
v^{2}\left(\frac{v^{2}-1}{v^{2}}\right)^{n}$ $\displaystyle\leq
x_{2}(1-x_{1})^{10}(1-x_{2})^{13}(1-x_{3})^{3},$ $\displaystyle
v^{3}\left(\frac{v^{3}-1}{v^{3}}\right)^{n}$ $\displaystyle\leq
x_{3}(1-x_{2})^{33}(1-x_{3})^{3}.$

No closed form is apparent so we solved this system numerically using OpenOpt
[openopt]. We provide the results in column $n_{g}$ of Table 4. For details of
the solution technique see [raaphorstphd].

Since the solutions for the system of equations from the general local lemma
can be quite difficult to produce [habib1998], we also consider the asymmetric
local lemma, as in Theorem 4.2. This gives the following system of
inequalities:

$\displaystyle v^{2}\left(\frac{v^{2}-1}{v^{2}}\right)^{n}$
$\displaystyle\leq\frac{1}{4},$ $\displaystyle
v^{3}\left(\frac{v^{3}-1}{v^{3}}\right)^{n}$ $\displaystyle\leq\frac{1}{4},$
$\displaystyle 26\left[v^{2}\left(\frac{v^{2}-1}{v^{2}}\right)^{n}\right]$
$\displaystyle\leq\frac{1}{8},$ $\displaystyle
33\left[v^{2}\left(\frac{v^{2}-1}{v^{2}}\right)^{n}\right]+3\left[v^{3}\left(\frac{v^{3}-1}{v^{3}}\right)^{n}\right]$
$\displaystyle\leq\frac{1}{8}.$

which can be solved more directly. Details can be found in [raaphorstphd]. The
results are given in column $n_{a}$ of Table 4.

$v$ | $n_{s}$ | $n_{g}$ | $n_{a}$ |  $p_{g,s}$ (%) |  $p_{a,s}$ (%)
---|---|---|---|---|---
2 | 50.10 | 33.79 | 34.38 | 32.55 | 31.37
3 | 209.50 | 148.30 | 153.17 | 29.21 | 26.88
4 | 556.87 | 407.02 | 421.87 | 26.91 | 24.24
5 | 1175.17 | 881.51 | 910.49 | 24.99 | 22.52
6 | 2152.02 | 1643.10 | 1693.86 | 23.65 | 21.28
7 | 3578.65 | 2777.33 | 2850.50 | 22.39 | 20.34
8 | 5549.38 | 4367.67 | 4461.93 | 21.29 | 19.59
9 | 8161.08 | 6440.68 | 6612.28 | 21.08 | 18.97
10 | 11512.91 | 9171.64 | 9387.96 | 20.34 | 18.45
Table 4:  Different values of $n$ obtained for $H_{15}$ using the symmetric,
general and asymmetric local lemmas.

Table 4 gives the values of $n$ obtained for $H_{15}$ using each of the three
versions of the local lemma. Column $p_{g,s}$ gives the percentage improvement
of the general local lemma with respect to the symmetric local lemma. Column
$p_{a,s}$ gives the percentage improvement of the asymmetric local lemma with
respect to the symmetric local lemma. The results obtained from the general
local lemma showed significant improvement over those from the symmetric local
lemma, with an average improvement of 24.71% and a median improvement of
23.65%. We note, however, that considerable work went into finding valid
functions $x\colon\mathcal{A}\rightarrow(0,1)$ satisfying the conditions in
Theorem 4.1. The process we used would be hard to automate and was highly
intensive, requiring significant manual experimentation and interaction. The
asymmetric local lemma gives results that are slightly worse than those given
by the general local lemma. On average, they are within 2.76% (median: 2.66%)
of those given by the general local lemma. While this difference is small, the
improvement that the use of the asymmetric local lemma gives over the use of
the symmetric local lemma is considerably more significant: the asymmetric
local lemma is, on average, an improvement of 22.63% (median: 21.29%) over the
symmetric local lemma. We believe that when the hypergraph for a variable
strength covering array lacks uniformity with respect to the size of edges or
the sizes of their neighbourhoods, there is benefit using the general and
asymmetric forms of the local lemma.

As a final comparison, we compare the size of covering arrays given by the
local lemma to the guaranteed bound of VarDens and also to results of running
VarDens. The results are shown in Figure 1. The best local lemma bounds (from
the general local lemma) are better than the theoretical guarantee of VarDens
($N_{dens}$) but the actual size of arrays given by running VarDens is even
better.

Figure 1:  A comparison of the VarDens bound (upper line), the general local
lemma bounds (middle line), and the actual results from running VarDens on
$H_{15}$ (lower line). The horizontal axis is the number of levels and the
vertical axis represents the natural logarithm of the size of the arrays.

##  5 Conclusion

We used variants of the Lovász local lemma to find upper bounds on the sizes
of variable strength covering arrays and compared to the ones constructed
using VarDens, a derandomized greedy construction, and its guaranteed upper
bound. Our main result, Theorem 1.3, is a general upper bound on the size of a
variable strength covering array in term of the parameters of the associated
hypergraph obtained via the symmetric local lemma. When nothing is known about
the hypergraph these two bounds are very similar.

The bounds obtained from the symmetric local lemma work best when the
hypergraph’s edges are of a fixed size and the number of edges intersecting an
edge is invariant. For example, for $t$-designs the local lemma bound is
better than the VarDens bound by a constant that depends on $t$. A more
extreme example are the cyclic consecutive hypergraphs for which the local
lemma bound, unlike VarDens bound, remains constant when we let $k$ grow. We
suggest that when edges are of varying sizes, the general and asymmetric
versions of the local lemma may work best, and demonstrate this fact with an
example. We note that in some instances the VarDens bound is much worse than
actual runs of VarDens on given inputs, and in several cases these runs
outperform the local lemma bound.

One direction for future research is examining how much of the recent
improvements in the application of both the local lemma and density based
greedy algorithm for standard covering arrays [MR3328867, francetic2016,
sarkar_upper_2017, colbourn_two-stage_nodate_alt, colbourn_asymptotic_nodate]
can be extended to variable strength covering arrays. Another direction is to
continue the exploration of the general and asymmetric local lemmas and their
utility for variable strength covering arrays. Perhaps some families of
hypergraphs are more amenable to the general case and closed form solutions
for the required function exist. The weighted local lemma is also deserving of
attention [habib1998].

Acknowledgements
Lucia Moura and Brett Stevens were supported by NSERC Discovery grants. The
authors would like to thank an anonymous reviewer for various suggestions that
improved the presentation of this paper.

## References


equations. By analyzing the system and its fixed points (when they exist) near
the flow direction in the FGP and performing numerical integration of the
trajectories starting near but not on the FGP in section 5.3, we determined
the orientational dynamics in the region nearby the flow gradient plane. The
particle either spirals away from or towards the FGP when there are no fixed
points in FGP. When the fixed points in FGP arise, the one closer to the flow
direction is either a saddle point (with its stable manifold in the FGP) or a
stable node. Therefore, the particle first migrates towards the flow direction
nearly parallel to the FGP. Then it may either settle near the flow direction
or monotonically drift away from the FGP along the flow-vorticity plane. The
boundaries in the $c\cdot De-\kappa$ space separating these different
qualitative behaviors were shown in figure 2 for $c=0.005$. The qualitative
nature of these boundaries is not sensitive to the exact value of $c$. As we
observed in figures 6, 7, 9, 10 and 11, there are other invariant features
(stable limit cycle, unstable limit cycle, saddle point, and unstable node)
that arise in the orientation space off the FGP and the vorticity axis. These
features may be viewed as a result of the orientation space being restricted
to a unit sphere.

In $R_{vort}^{(1)}$ (figure 2), which extends to arbitrarily small (but
finite) viscoelasticity (small $c\cdot De$), the vorticity axis is a spiral
source, and the flow-gradient plane is an unstable limit cycle. Therefore, at
least one stable attractor must exist between these two regions on the unit
sphere. In the case of a Newtonian fluid, the particle follows an initial
condition dependent Jeffery orbit, i.e., one of a concatenation of neutrally
stable periodic orbits centered at the vorticity axis. Therefore, at small
(but finite) $c\cdot De$ the simplest change in phase-plane dynamics leads to
a stable limit cycle around and near one of these limit cycles, as shown in
figure 6. As shown by linear stability analysis at the vorticity direction in
section 5.2 the rate of deviation of trajectories away from the vorticity axis
is driven by the $\mathcal{O}(1/\kappa^{2})$ flow generated by the force
doublet and $\mathcal{O}(1/\kappa^{2})$ Stokeslet. Thus, a slender particle
starting at any orientation between the vorticity axis and the flow-gradient
plane leads to a final orientation behavior where the particle undergoes
periodic motion very close to and around the vorticity axis.

In $R_{flow-vort}^{(1)}$ (figure 2), the FGP changes from being an unstable to
a stable limit cycle. Therefore, an unstable invariant object or a repeller in
the form of an unstable limit cycle exists between the stable limit cycle near
the vorticity direction and the FGP, as shown in figure 7. Thus in $R_{flow-
vort}^{(1)}$, a particle with an initial orientation between the unstable
limit cycle and FGP spirals to the FGP, where it undergoes tumbling motion. A
particle with an initial orientation between the unstable limit cycle above
the FGP and the stable limit cycle near the vorticity axis spirals towards the
stable limit cycle near the vorticity axis. A particle released at an
arbitrarily small (but finite) angle from the vorticity axis spirals outwards.
In both these cases, eventually, the particle undergoes perpetual periodic
motion on the stable limit cycle close to the vorticity axis.

In $R_{vort}^{(2)}$ (figure 2), the FGP contains two fixed points as shown in
figure 9. One is an unstable node, and the other is a saddle point with its
unstable manifold perpendicular to the FGP. Thus, no other invariant feature
(attractor or repeller) is needed between the FGP and the stable limit cycle
near the vorticity axis, which is the only stable invariant object on the
orientation space. Hence a particle starting at an arbitrarily small angle
from the vorticity axis or the FGP ends up in a periodic motion very close to
the vorticity axis.

Similar to $R_{vort}^{(1)}$, $R_{vort}^{(2)}$ and $R_{flow-vort}^{(1)}$
discussed above, there is a stable limit cycle around the vorticity axis in
$R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$. Therefore, a particle with an
initial orientation arbitrarily close to the vorticity direction (but not on
the vorticity axis) eventually undergoes a periodic motion around the
vorticity axis. This stable limit cycle exists for all regions except for
$R_{flow}^{(1)}$ and $R_{flow}^{(2)}$ and it will be further analyzed in
section 5.4.1.

In $R_{flow-vort}^{(2)}$ (figure 2) the fixed point in the FGP closest to the
flow direction is a stable node and the other fixed point in the FGP is an
unstable node as shown in figure 10. Hence, a saddle node exists near the
flow-vorticity plane to repel the particle orientation trajectories towards
the stable limit cycle near the vorticity axis on one side and the stable
fixed point near the flow direction on the FGP on the other. This saddle point
receives trajectories emanating from the unstable point on the FGP. Thus,
depending upon the initial orientation of the particle, it may eventually
either obtain a stable orientation near the flow direction or undergo a
periodic motion close to the vorticity axis.

The dynamics in $R_{flow-vort}^{(3)}$ (figure 2) are similar to those in
$R_{flow-vort}^{(2)}$, but the unstable fixed point in the FGP in $R_{flow-
vort}^{(2)}$ changes to a saddle-node in $R_{flow-vort}^{(3)}$ with no
qualitative change to the other invariant objects present in $R_{flow-
vort}^{(2)}$. Therefore, to allow smooth and continuous phase flow, an
unstable node occurs at the intersection of the saddle points in the FGP and
near the flow-vorticity plane, as shown in figure 11. At a given $\kappa$, the
value of $c\cdot De$ in $R_{flow-vort}^{(3)}$ is larger than that in $R_{flow-
vort}^{(2)}$ (figure 2) and the separation between two fixed points increases
with increasing $c\cdot De$ as shown by figure 3. Due to this increased
separation and the unstable node above the flow gradient plane, more of the
phase flow is directed toward the FGP in $R_{flow-vort}^{(3)}$ than in
$R_{flow-vort}^{(2)}$. Thus, the basin of attraction for the stable fixed
point near the flow direction is larger in $R_{flow-vort}^{(3)}$ as compared
to $R_{flow-vort}^{(2)}$. In other words more of the initial particle
orientations lead to a stable final orientation near the flow direction in the
$R_{flow-vort}^{(3)}$ than in $R_{flow-vort}^{(2)}$ as observed by comparing
the proportion of solid lines in figures 10 and 11. When the stable limit
cycle near the vorticity axis exists, the saddle point on the flow-vorticity
plane is important in determining the proportions of initial orientations
leading to the flow alignment and towards the periodic motion around the
vorticity axis. Its location will be further analyzed in section 5.4.2.

The $c\cdot De|_{\text{cut-off}}^{\text{Limit}}$ boundary in figure 2, beyond
which the stable limit cycle near the vorticity direction does not exist, may
be viewed as arising due to the interaction between the saddle-node near the
flow-vorticity plane and the stable limit cycle near the vorticity direction
as these move in the orientation space upon increasing $c\cdot De$. We will
discuss this in section 5.4.3.

####  5.4.1 Stable Limit cycle around the vorticity axis

The previous study of Leal (1975) concerning a slender particle in simple
shear flow of a second-order fluid predicts that viscoelasticity (although
incorrect in attributing it to the second instead of the first normal stress
difference) drives the particle orientation towards the vorticity axis. He
only consider the Stokeslet fluid flow (of order $\mathcal{O}(1/\log(\kappa))$
and $\mathcal{O}(1/\log(\kappa)^{2})$) from the slender body theory. This
velocity field is proportional to $p_{2}$ and hence has no effect when the
particle is in the flow-vorticity plane. However, our study shows that a
stable limit cycle exists around the vorticity axis due to the competition
between polymer-driven torque arising from the doublet flow causing the
particle to spiral away from the vorticity axis and that from the Stokeslet
flow at $\mathcal{O}(1/\log(\kappa))$ causing it move towards the region
around vorticity.

The real part of the eigenvalues of the fixed point at the vorticity axis in
$R_{vort}^{(1)}$ is $2cDe/\kappa^{2}$ (equation (111)), and the rate at which
trajectories leave the flow-gradient plane is $cDe\zeta$ (equation (114)).
Figure 12 shows the variation of $2cDe/\kappa^{2}$ and $cDe\zeta$ with
$\kappa$ for a few values of $c\cdot De$ at $c=0.005$. From this figure, we
conclude that for $\kappa\gtrapprox 15$, trajectories leave the flow gradient
plane much faster than the outward spiraling rate from the vorticity axis.
Hence, the stable limit cycle shifts closer to the vorticity axis as $\kappa$
increases. Furthermore, as $\kappa$ increases beyond $\kappa\approx 20$,
$cDe\zeta$ remains almost constant while the rate of spiralling from the
vorticity axis, $\kappa^{-2}$, decreases. Therefore, the stable limit cycle
approaches the vorticity axes at larger $\kappa$.

Figure 12: The variation with $\kappa$ of the rate of deviation of orientation
trajectories away from the vorticity axis, $2cDe/\kappa^{2}$ (orange), and the
flow gradient plane, $cDe\zeta$ (black), in $R_{vort}^{(1)}$ ($b^{2}>0$) for a
few values of $c\cdot De$ at $c=0.005$. On each curve, $b^{2}>0$ is to the
left of the solid markers. An increasing gap between orange and black curves,
with almost horizontal black curves, with $\kappa$ indicates an larger
increase in the departure rate of the orientation trajectories from the flow
gradient plane than that from the vorticity axis, implying a shift of the
stable limit cycle closer to the vorticity axis.

We quantify the location of the limit cycle observed in the numerical
integration of equation (95). For this purpose we transform this equation into
the $C-\tau$ coordinate system (Leal & Hinch, 1971), using,

$C=\frac{\sqrt{p_{2}^{2}+p_{1}^{2}/\kappa^{2}}}{p_{3}},\tau=\frac{p_{1}}{p_{2}},$
(118)

to obtain,

$\displaystyle\frac{dC}{dt}=cDe\tau^{2}\Big{(}\frac{4(1+C^{2})}{C\kappa^{4}}-\frac{C}{(4\log(2\kappa)-3)(\kappa^{2}+\tau^{2})}\Big{)}p_{2}^{2}+cDe\mathcal{O}(p_{2}^{4}),$
(119) $\displaystyle\frac{d\tau}{dt}=1+cDe\tau
p_{2}^{2}\Big{(}\frac{\tau^{2}}{4\log(2\kappa)-3}+4\frac{\kappa^{2}+\tau^{2}}{C^{2}\kappa^{4}}\Big{)}+\mathcal{O}(c\cdot
De)$ (120)

In the Newtonian limit, ${dC}/{dt}=0$ and the particle trajectory is
determined by its initial condition $C(t)=C(0)=C_{N}$. Each $C_{N}=[0,\infty]$
represents a different Jeffery orbit. $C=0$ represents the log-rolling motion,
i.e., the particle rotating about its major axis which is aligned with the
vorticity axis, and $C=\infty$ represents the major axis rotating within the
flow-gradient plane. $\tau=t$ represents the phase on the Jeffery orbit. A
periodic orbit can thus be ascertained using a single parameter $C$ that
repeats periodically. In the range, $10\leq\kappa\leq 200$, for $c\cdot De\geq
0.01$ we evolve the system of equations for an initial condition starting near
the vorticity direction, at $C=10^{-7}$. We evolve the equations for each
choice of $c\cdot De$ and $\kappa$ until the change in the period-averaged $C$
across 100 successive Jeffery periods is less than $10^{-5}$. Finally, we
define $\bar{C}_{\text{Limit}}^{\text{numerical}}$ as the average $C$ over the
last 100 Jeffery periods. $\bar{C}_{\text{Limit}}^{\text{numerical}}$ vs.
$c\cdot De$ for a few values of $\kappa\in[10,200]$ is shown in the left panel
of figure 13.

(a)  (b)

Figure 13: Left: The location of the stable limit cycle at different $c\cdot
De$ and $\kappa$ at $c=0.005$ (nearly identical curves are obtained for other
polymer concentrations in the range $0.01\leq c\leq 0.2$) indicated by the
average $C=\sqrt{p_{2}^{2}+p_{1}^{2}/\kappa^{2}}/{p_{3}}$,
$\bar{C}_{\text{Limit}}^{\text{numerical}}$ on the limit cycle. $C=0$ is the
vorticity axis and $C=\infty$ is the flow-gradient plane. Right:
$\bar{C}_{\text{Limit}}^{\text{numerical}}$ on the limit cycle vs. $\kappa$ at
low $c\cdot De$.

$\bar{C}_{\text{Limit}}^{\text{numerical}}$ is smaller than $0.1$ for
$\kappa>10$, decreases with $\kappa$ and varies slightly with $c\cdot De$ at
each $\kappa$. For $\kappa>50$ $\bar{C}_{\text{Limit}}^{\text{numerical}}$ is
less than $0.01$ for a range of $c\cdot De$. Therefore, the limit cycle is
close to the vorticity axis. Beyond a certain value of $c\cdot De=c\cdot
De|_{\text{cut-off}}^{\text{Limit}}$, the limit cycle does not exist. This
marks the boundary between $R_{flow-vort}^{(2)}$ and $R_{flow}^{(1)}$ for a
given $\kappa$, as shown in figure 2. $c\cdot De|_{\text{cut-
off}}^{\text{Limit}}$ increases with $\kappa$ and for $10\leq\kappa\leq 200$
we find,

$c\cdot De|_{\text{cut-off}}^{\text{Limit}}\approx 0.1056\kappa-0.5183.$ (121)

At larger $\kappa$, the stable limit cycle is closer to the vorticity axes and
exists for a larger range of $c\cdot De$. Hence, marked as $R_{flow-
vort}^{(1)}$ and $R_{flow-vort}^{(2)}$ in figure 2, there is a significant
range of $c\cdot De$ at large $\kappa$ where two possible types of orientation
dynamics are possible: a steady state orientation close to the flow direction
and a periodic orbit near the vorticity direction. For $\kappa\lessapprox 20$,
increasing $c\cdot De$ from very small values does not affect the position of
the limit cycle much. However, as $c\cdot De$ approaches $c\cdot
De|_{\text{cut-off}}^{\text{Limit}}$, the limit cycle moves away from the
vorticity axis. For $\kappa\gtrapprox 20$, the limit cycle moves towards the
vorticity axis upon increasing $c\cdot De$ from very small values. However,
beyond a given $c\cdot De$, up to $c\cdot De|_{\text{cut-off}}^{\text{Limit}}$
the limit cycle’s angular position is not influenced by $c\cdot De$. These
features are reflected by the trends in
$\bar{C}_{\text{Limit}}^{\text{numerical}}$ in the left panel of figure 13. At
low $c\cdot De$, the limit cycle is very close to one of the degenerate
Jeffery orbits of the Newtonian case, i.e., $C$ does not change much on the
limit cycle. The variation of $\bar{C}_{\text{Limit}}^{\text{numerical}}$ with
$\kappa$ at low $c\cdot De=0.01,0.02$ and 0.03 is shown in the right panel of
figure 13, where we can observe $\bar{C}_{\text{Limit}}^{\text{numerical}}$ to
scale as $\kappa^{-1.5}\log(2\kappa-3)$ for $10\leq\kappa\leq 100$ and
approximately $\kappa^{-1}$ for larger $\kappa$.

We find the particle orientation trajectories to approach this stable limit
cycle close to the vorticity axis instead of spiraling towards the vorticity
axis as concluded from the experiments of Gauthier et al. (1971) and Iso et
al. (1996b). Gauthier et al. (1971) claim the approach toward the vorticity
axis by extrapolating the observed data. In some of the experimental
observations of Iso et al. (1996b) at small elasticity after undergoing
initial spiraling, the particle oscillates around the vorticity axis without
settling at a steady orientation. This may indicate the limit cycle discussed
above. The theoretical rate of spiraling away from the vorticity axis
($2cDe/\kappa^{2}$) is small (figure 12), and small imperfections from the
Couette cell might have led to secondary flows in the experiments. Thus, the
existence and size of the stable limit cycle around the vorticity axis must be
further tested in future experiments and numerical simulations.

####  5.4.2 Saddle node near the flow-vorticity plane in $R_{flow-
vort}^{(2)}$, $R_{flow-vort}^{(3)}$, $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$

From the trajectories of particle orientation in $R_{flow-vort}^{(2)}$ and
$R_{flow-vort}^{(3)}$ shown in figures 10 and 11 we find the stable limit
cycle around the vorticity axis and the stable node on the flow gradient plane
near the flow axis to be separated by a saddle point. The saddle point arises
on the boundary between $R_{vort}^{(2)}$ and $R_{flow-vort}^{(2)}$ when the
saddle point near the flow direction on the flow-gradient plane changes into a
stable node. It also arises on the boundary between $R_{flow-vort}^{(3)}$ and
$R_{flow-vort}^{(1)}$ when the unstable limit cycle in the flow-gradient plane
vanishes. The saddle point persists in the region $R_{flow-vort}^{(3)}$ when
approaching it by increasing $c\cdot De$ from $R_{flow-vort}^{(2)}$. In the
regions $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$, the unstable
direction of the saddle point leads to the limit cycle on one side, and the
stable node on the other. Based on numerical evidence from figures 10 and 11
we assume the saddle point,
$\mathbf{p}_{\text{saddle}}=\begin{bmatrix}{p}_{\text{1,saddle}}&{p}_{\text{2,saddle}}&{p}_{\text{3,saddle}}\end{bmatrix}$,
to approximately lie at the same $p_{2}$ coordinate as the stable node in the
flow gradient plane,

${p}_{\text{2,saddle}}={p}_{2,\text{flow}}^{0-}+{p}_{\text{2,saddle}}^{\prime},\hskip
14.45377pt{p}_{\text{2,saddle}}^{\prime}\ll{p}_{2,\text{flow}}^{0-}$ (122)

where ${p}_{2,\text{flow}}^{0-}$ is given by equation (108). From
$\dot{p}_{\text{2,saddle}}=0$, we obtain,

${p}_{\text{2,saddle}}=\frac{\sqrt{-b^{2}}c\cdot
De+(4\log(2\kappa)-6)b^{2}}{c\cdot
De(1-p_{1,saddle})-\sqrt{-b^{2}}(4\log(2\kappa)-6)}.$ (123)

From equation (95), $\dot{p}_{\text{1,saddle}}=\dot{p}_{\text{3,saddle}}=0$ is
equivalent to,

$p_{\text{1,saddle}}\alpha|_{p_{1}=p_{\text{1,saddle}},p_{3}=p_{\text{3,saddle}}}=p_{\text{2,saddle}}/(c\cdot
De).$ (124)

A general solution to the above equation is intractable, but we can solve it
in the limit of small $p_{\text{1,saddle}}$ or small $p_{\text{3,saddle}}$. We
denote the value of $p_{\text{1,saddle}}$ in the limit of small
$p_{\text{1,saddle}}$ as $\tilde{p}_{\text{1,saddle}}$ and the value of
$p_{\text{3,saddle}}$ in the limit of small $p_{\text{3,saddle}}$ as
$\tilde{p}_{\text{3,saddle}}$. We obtain,

$\tilde{p}_{\text{1,saddle}}=\frac{\begin{split}&-8(c\cdot
De)^{2}f^{2}+16c\cdot Def^{3}\sqrt{-b^{2}}\\\ &-\sqrt{\begin{split}&4c\cdot
Def\left(-(c\cdot De)^{2}\kappa^{2}+2(c\cdot
De)f\kappa^{2}\sqrt{-b^{2}}+4f^{2}\right)\\\ &4\left(-3\pi(c\cdot
De)^{2}\kappa^{2}+2c\cdot Def\left(8f+3\pi\kappa^{2}\sqrt{-b^{2}}\right)+12\pi
f^{2}\right)\\\ &+128c\cdot Def^{3}((c\cdot
De)-2f\sqrt{-b^{2}})^{2}\end{split}}\end{split}}{(c\cdot De)\left(3\pi(c\cdot
De)^{2}\kappa^{2}-2(c\cdot
De)f\left(8f+3\pi\kappa^{2}\sqrt{-b^{2}}\right)-12\pi f^{2}\right)},$ (125)

and,

$\tilde{p}_{\text{3,saddle}}=-\frac{\begin{split}&3\pi(c\cdot De)\left((c\cdot
De)^{2}\kappa^{2}-2c\cdot Def\kappa^{2}\sqrt{-b^{2}}-4f^{2}\right)\\\
&+\sqrt{\begin{split}&9\pi^{2}(c\cdot De)^{2}\left(-(c\cdot
De)^{2}\kappa^{2}+2c\cdot Def\kappa^{2}\sqrt{-b^{2}}+4f^{2}\right)^{2}\\\
&+\frac{1}{2}(c\cdot De)^{3}\kappa^{2}\left(3\pi^{2}c\cdot
De-2f\left(3\pi^{2}\sqrt{-b^{2}}+2\right)\right)\\\ &+\frac{1}{2}c\cdot
Def^{2}\left(f(16-64c\cdot De\sqrt{-b^{2}})+c\cdot
De\left(8\kappa^{2}\sqrt{-b^{2}}-12\pi^{2}\right)\right)\\\
&\left(3\pi^{2}(c\cdot De)^{2}\kappa^{2}-2c\cdot
Def\left(16f+3\pi^{2}\kappa^{2}\sqrt{-b^{2}}\right)-4f^{2}\left(16f\sqrt{-b^{2}}+3\pi^{2}\right)\right)\end{split}}\end{split}}{2\left(\frac{3}{4}\pi^{2}c\cdot
De\left((c\cdot De)^{2}\kappa^{2}-2c\cdot
Def\kappa^{2}\sqrt{-b^{2}}-4f^{2}\right)-8c\cdot Def^{2}((c\cdot
De)+2f\sqrt{-b^{2}})\right)},$ (126)

where $f=2\log(2\kappa)-3$ and $\sqrt{-b^{2}}$ is in equation (115). The
variation of $\tilde{p}_{\text{3,saddle}}$ and
$\sqrt{1-\tilde{p}_{\text{1,saddle}}^{2}-{p}_{\text{2,saddle}}^{2}}$ with
$c\cdot De$ for a few $\kappa$ is shown in the left panel of figure 14. From
either of these approximations, we find that the saddle point moves closer to
the vorticity axis as $c\cdot De$ is increased or $\kappa$ is reduced. Both
approximations yield a similar location for the saddle point. The best
approximation for the location of the saddle point for any $c\cdot De$ and
$\kappa$, is
$\mathbf{p}_{\text{saddle}}=\begin{bmatrix}{p}_{\text{1,saddle}}&{p}_{\text{2,saddle}}&{p}_{\text{3,saddle}}\end{bmatrix},$
where,

$\displaystyle\begin{split}{p}_{\text{1,saddle}}=&\tilde{p}_{\text{1,saddle}},\hskip
7.22743pt{p}_{\text{3,saddle}}=\sqrt{1-\tilde{p}_{\text{1,saddle}}^{2}-{p}_{\text{2,saddle}}^{2}},\hskip
14.45377pt\text{if
}|\tilde{p}_{\text{1,saddle}}|\leq|\tilde{p}_{\text{3,saddle}}|\\\
{p}_{\text{1,saddle}}=&\sqrt{1-{p}_{\text{2,saddle}}^{2}-\tilde{p}_{\text{3,saddle}}^{2}},\hskip
7.22743pt{p}_{\text{3,saddle}}=\tilde{p}_{\text{3,saddle}},\hskip
14.45377pt\text{if
}|\tilde{p}_{\text{1,saddle}}|>|\tilde{p}_{\text{3,saddle}}|\end{split},$
(127)

and $p_{2}$ is from equation 123. We show this saddle point’s more accurate
location of $p_{3}$ in the right panel of figure 14. The location of this
analytically predicted saddle point, shown with a green marker in figures 10,
11 and 15, agrees well with that inferred from the of numerically integrated
trajectories.

In $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$, the unstable manifold of
this saddle point leads to two stable attractors, i.e., either a limit cycle
near the vorticity axis or a stable node near the flow direction. Therefore,
the location of the saddle point partially dictates the relative sizes of the
basins of attraction of these stable attractors. As $c\cdot De$ increases, due
to the changing location of the saddle point, $\mathbf{p}_{\text{saddle}}$,
the particle is more likely to attain a final orientation towards the flow
axis as compared to approaching the limit cycle around the vorticity axis (see
figures 14 and compare the proportion of the solid lines in 10 and 11). The
saddle point also persists in the large $c\cdot De$ regions $R_{flow}^{(1)}$
and $R_{flow}^{(2)}$ as shown in the different plots in figure 15.

(a)  (b)

Figure 14: $p_{3}$ coordinate of the saddle point in the flow-vorticity plane
in $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$. Left: Solid lines (the
graph of $\sqrt{1-\tilde{p}_{\text{1,saddle}}^{2}-p_{\text{2,saddle}}^{2}}$
from equation (125) and (123)) provide the more accurate $p_{3}$ location of
the saddle point when it is closer to the vorticity axis
($p_{\text{3,saddle}}\approx 1$) and dashed lines (the graph of equation
(126)) represent the more accurate $p_{3}$ location of the saddle point when
it is closer to the flow axis ($p_{\text{3,saddle}}\approx 0$). Right: most
accurate location of the saddle point constructed from the left figure using
the cross over point between the solid and dashed lines for each $\kappa$ i.e.
using equation (127).

####  5.4.3 $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$: Very large $c\cdot De$
guarantees flow alignment

In section 5.4.1, we found that for $\kappa\gtrapprox 20$, the stable limit
cycle moves closer to the vorticity axis upon increasing $c\cdot De$. The
saddle point near the flow-vorticity plane that exists in regions shown in
figure 2 with $c\cdot De$ equal and larger than that in $R_{flow-vort}^{(2)}$
also moves towards the vorticity axis upon increasing $c\cdot De$ as discussed
in section 5.4.2. The saddle point moves towards the vorticity axis faster
than the limit cycle. Beyond a $\kappa$-dependent value of $c\cdot De$, when
the saddle point is too close to the limit cycle, the latter ceases to exist,
making the stable fixed point near the flow direction the only stable
attractor. This boundary is marked as $c\cdot De|_{\text{cut-
off}}^{\text{Limit}}$ in figure 2. This implies that a particle with any
initial orientation other than the vorticity axis ends up being aligned near
the flow direction. Plots corresponding to $c\cdot De=1.2$ and 1.8 for
$\kappa=20$ in figure 15 show the vanishing of the stable limit cycle near the
vorticity axis upon going from $R_{flow-vort}^{(3)}$ to $R_{flow}^{(1)}$. For
$c\cdot De=1.2$ ($R_{flow-vort}^{(3)}$) the saddle point is just outside the
limit cycle, whereas for $c\cdot De=1.8$ ($R_{flow}^{(1)}$), corresponding to
a location close to the $R_{flow-vort}^{(3)}-R_{flow}^{(1)}$ boundary, the
limit cycle does not exist and instead the stable manifolds of the saddle
point originate from the unstable spiral at the vorticity axis. The unstable
manifold of the saddle point leads to the stable node near the flow direction.
The saddle point persists for larger $c\cdot De$ in $R_{flow}^{(1)}$ and also
$R_{flow}^{(2)}$ as show by the $c\cdot De=4$ and 12 plots at $\kappa=20$
respectively in figure 15. As also discussed in section 5.2 (left panel of
figure 4) in $R_{flow}^{(2)}$, a particle leaving the vorticity axis changes
from spiralling to monotonic trajectories as the vorticity axis changes from
an unstable spiral to an unstable node when going from $R_{flow}^{(1)}$ to
$R_{flow}^{(2)}$. From the $c\cdot De=4$ and 12 panels of figure 15 at
$\kappa=20$, the proximity of the saddle point to the vorticity axis can also
explain the breakdown of spiral motion near the vorticity axis in
$R_{flow}^{(2)}$.

(a)  (b)  (c)  (d)

Figure 15: Trajectories at large $c\cdot De$ for $\kappa=20$ in $R_{flow-
vort}^{(3)}$, $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$ zoomed near the flow-
vorticity plane. At very large $c\cdot De$ (corresponding to regions
$R_{flow}^{(1)}$ and $R_{flow}^{(2)}$), the stable limit cycle does not exist
around the vorticity axis and a particle starting anywhere apart from
vorticity axis ends up being nearly flow aligned.

##  6 Conclusion

Using a regular perturbation expansion in polymer concentration, $c$, the
balance of torques on the particle surface at each order in $c$, and the
velocity field generated by a slender prolate spheroid (obtained from the
slender body theories of Batchelor (1970) and Cox (1971)), we develop a theory
to characterize the orientation dynamics of a freely rotating (torque-free)
particle in simple shear flow of a viscoelastic fluid with small polymer
concentration, $c$, and an arbitrary polymer relaxation time in the absence of
inertia. This theory predicts a wide variety of orientation behaviors that are
qualitatively similar to previous experimental observations of Gauthier et al.
(1971); Bartram et al. (1975); Iso et al. (1996a, b).

In the absence of inertia, for a viscoelastic fluid, where the fluid stress is
a sum of the solvent and the polymer or non-Newtonian stress, the momentum
equation can be decomposed into a Newtonian and a non-Newtonian part. The
Newtonian part comprises the particle’s motion in a Newtonian fluid undergoing
the imposed flow and leads to the Newtonian stress field. The non-Newtonian
part has zero velocity at the boundaries and comprises the balance of the
divergence of the polymer induced solvent stress and the divergence of the
polymer stress. Thus, three physically distinct stress mechanisms impacting
the particle surface are the elastic or polymer stress, the polymer-induced
solvent stress, and the Newtonian stress. The decomposition of torques in this
way can be useful in obtaining further insights into the particle motion in a
viscoelastic or a non-Newtonian fluid where the fluid stress can be decomposed
into a Newtonian and a non-Newtonian part. In addition to the non-Newtonian
momentum equation, different stress components are coupled with the torque-
free boundary condition and the polymer constitutive equation. The balance of
torques generated by the three stresses leads to the appropriate angular
velocity to allow torque-free particle motion. The boundary conditions for the
motion of the particle and the imposed flow are contained in the Newtonian
part of the momentum equation. Hence, the torque generated from the
corresponding stress is termed the motion induced solvent torque (MIST). The
polymer constitutive equation is driven by the sum of the velocity field from
the Newtonian and non-Newtonian momentum components. The elastic torque is a
function of the polymer stress on the particle surface. By definition, the
polymer-induced solvent torque (PIST) is the antisymmetric first moment of the
polymer-induced solvent stress over the particle surface . However, using a
generalized reciprocal theorem and the non-Newtonian momentum equation, we can
express the PIST directly as a volume integral of the polymer stress.

Using a regular perturbation in the polymer concentration, $c$, we obtain the
equations for particle’s orientation dynamics in a small $c$ viscoelastic
fluid subjected to a simple shear flow. At $\mathcal{O}(1)$, the elastic and
the polymer induced solvent stress are zero. Thus, the particle orientation
dynamics are the same as that in the simple shear flow of a Newtonian fluid,
i.e., the particle undergoes Jeffery (1922) rotations. At $\mathcal{O}(c)$,
all three mechanisms mentioned above lead to a finite torque. The
$\mathcal{O}(c)$ Newtonian stress is that on a particle rotating (at
$\mathcal{O}(c)$ velocity) in a quiescent Newtonian fluid. Therefore, the
$\mathcal{O}(c)$ rotation rate is the one that generates enough
$\mathcal{O}(c)$ MIST to balance the PIST and the elastic torque. The total
$\mathcal{O}(c)$ rotation rate may be decomposed as the sum of
$\dot{\mathbf{p}}^{(1)}_{\text{Elastic}}$ and
$\dot{\mathbf{p}}^{(1)}_{\text{PIST}}$ i.e. the $\mathcal{O}(c)$ elastic and
PIST generated rotation rates. The leading order velocity field is taken from
the slender body theory of Batchelor (Batchelor, 1970) and Cox (Cox, 1971),
and this drives the $\mathcal{O}(c)$ polymer stress. As both the elastic
torque and the PIST can be expressed as a function of the polymer stress, the
$\mathcal{O}(c)$ rotation rate of the particle is obtained from the analytical
velocity field from the slender body theory.

The polymer relaxation time ($\lambda$) is non-dimensionalized with the shear
rate ($\dot{\gamma}$) to yield Deborah number, $De(=\lambda\dot{\gamma})$
($c=0$ or $De=0$ implies a Newtonian fluid). We find
$\dot{\mathbf{p}}^{(1)}_{\text{Elastic}}$, to be independent of $De$ and
initially consider $\dot{\mathbf{p}}^{(1)}_{\text{PIST}}$, separately in the
low and high $De$ limits. In the low $De$ limit,
$\dot{\mathbf{p}}^{(1)}_{\text{Elastic}}$ and $\mathcal{O}(De^{0})$
$\dot{\mathbf{p}}^{(1)}_{\text{PIST}}$ are the rotation rates due to the
additional torque from the rate of strain and fluid pressure respectively of
the Newtonian fluid with an additional viscosity $c$. Therefore, in the low
$De$ limit, the sum of $\dot{\mathbf{p}}^{(1)}_{\text{Elastic}}$ and the
$\mathcal{O}(De^{0})$ contribution to $\dot{\mathbf{p}}^{(1)}_{\text{PIST}}$
is identical to the rotation due to the Newtonian torque. Since the Newtonian
torque is already balanced to be zero at the leading order in $c$, the first
viscoelastic effect on the rotation rate for $De\ll 1$ arises at
$\mathcal{O}(c.De)$ and is entirely from the polymer-induced solvent stress.
In the high $De$ limit,
$\dot{\mathbf{p}}^{(1)}_{\text{PIST}}\sim\mathcal{O}(De)$ and is hence
$\mathcal{O}(De)$ larger than $\dot{\mathbf{p}}^{(1)}_{\text{Elastic}}$.
Therefore, the particle-induced polymer stress is the primary agent changing
the particle dynamics in both the low and high $De$ limits. By interpolating
the rotation rate in $De$ between the large and small (up to
$\mathcal{O}(De)$) $De$ limits, we obtain a uniformly valid equation for the
particle orientation. We analyze the effect of particle aspect ratio, $\kappa$
and viscoelastic fluid parameters, polymer concentration ($c$), and relaxation
rate ($De$) on this equation valid for any $De$.

Depending upon $c\cdot De$ and particle aspect ratio ($\kappa$), several
qualitatively different particle orientation behaviors are possible. At
constant $c\cdot De$, $c$ affects the behavior only quantitatively. In a
Newtonian fluid ($c\cdot De=0$), a particle undergoes periodic motion on an
initial condition dependent Jeffery (1922) orbit around the vorticity axis. As
expected from symmetry, the vorticity axis and the flow-gradient plane (FGP)
remain invariant within the orientational space (unit sphere). Viscoelasticity
(for all $c\cdot De$ and $\kappa$) does not affect the rotation rate of the
log-rolling state on the vorticity axis. In the FGP, even very small $c\cdot
De$ creates a bottleneck near the flow direction where the particle slows down
in proportion to $c\cdot De$. Thus viscoelasticity increases the period of
tumbling motion within the FGP. At $c\cdot
De=cDe_{crit}=(4\log(2\kappa)-6)/\kappa$, an infinite period bifurcation
occurs, and two fixed points arise in the FGP near the flow direction.
Therefore, for $c\cdot De>cDe_{crit}$, the tumbling motion within the FGP does
not exist, and a particle initially oriented within FGP ultimately aligns near
the flow direction.

In a Newtonian fluid, the phase flow of the dynamical system of the particle’s
orientational drift is of equal and opposite magnitude in the vorticity
direction about the flow-vorticity plane. Thus the Jeffery (1922) orbits, in
that case, are symmetric about the flow-vorticity plane. For $0<c\cdot
De<cDe_{crit}$ and $\kappa>\kappa_{FV1}$, the phase flow remains qualitatively
similar, but the magnitude of the phase velocity in the vorticity direction is
larger on the downstream side of the phase-flow. Thus, a particle starting
close to (but not in) the flow-gradient plane spirals towards the vorticity
direction. $\kappa_{FV1}\approxeq 17$ for $c=0.005$ and slightly reduces with
$c$. For $\kappa<\kappa_{FV1}$, spiraling away from the FGP occurs for a
$c\cdot De$ up to a value slightly less than $cDe_{crit}$ indicated by the
$\zeta=0$ curve in figure 2. $R_{vort}^{(1)}$ in figure 2 is the region where
the particle spirals away from the FGP. It continues to spiral towards the
vorticity axis as it moves away from the FGP until it comes close to the axis.
A particle with an initial orientation arbitrarily close to, but not on, the
vorticity axis spirals away from the vorticity axis. Thus, in
$R_{vort}^{(1)}$, a particle ultimately undergoes a periodic motion in a limit
cycle near the vorticity axis. The spiraling towards the vorticity axis due to
viscoelasticity is observed in the experiments of Gauthier et al. (1971) and
Iso et al. (1996a, b).

For $\kappa<\kappa_{FV1}$ and between $c\cdot De$ corresponding to $\zeta=0$
and $c\cdot De=cDe_{crit}$, marked as $R_{flow-vort}^{(1)}$ in figure 2, the
phase velocity in the vorticity direction downstream of the flow-vorticity
plane and close to the FGP changes direction as compared to the
$R_{vort}^{(1)}$ region discussed above or the case for a Newtonian fluid.
Thus, a particle with an initial orientation close to the FGP spirals into the
FGP, where it ultimately undergoes tumbling motion (albeit with a larger
period than in Newtonian fluid). A particle starting further away from the FGP
or near the vorticity axis spirals towards the periodic orbit (similar to
$R_{vort}^{(1)}$) near the vorticity axis.

The phase space for the dynamical system of a particle in a Newtonian fluid
consists of a concatenation of neutral centers. The periodic orbit at small
angles away from the vorticity axis for the case of a viscoelastic fluid
($c\cdot De>0$) is a stable limit cycle. It is the only stable attractor in
the phase space in $R_{vort}^{(1)}$. As $c\cdot De$ increases at large
$\kappa\gtrapprox 50$ in the small $c\cdot De$ regime, the limit cycle shrinks
as it goes towards the vorticity axis. The size of the limit cycle is not
affected by $c\cdot De$ in the small $c\cdot De$ regime for $\kappa\lessapprox
20$. The previous studies by Leal (1975) using the slender body theory at low
$De$ have considered only the $\mathcal{O}(1/(2\log(2\kappa)-3))$ velocity
disturbance and found viscoelasticity to cause a slender particle to spiral
towards the vorticity axis. However, the polymeric torque due to the force
doublet and $\mathcal{O}(1/\kappa^{2})$ Stokeslet, $\mathcal{O}(1/\kappa^{2})$
velocity disturbance from slender body theory of Cox (1971) that accounts for
the fluid velocity disturbance the particle makes when in the flow-vorticity
plane ($p_{2}=0$), makes the vorticity axis an unstable spiral. Thus, a stable
limit cycle occurs from the competition of polymeric torques from the
disturbance at $\mathcal{O}(1/\kappa^{2})$ and the Stokeslet flow at order
$\mathcal{O}(p_{2}/(2\log(2\kappa)-3))$. Unlike Leal (1975) who attributes it
to second normal stress difference, our theory predicts the effect of
viscoelasticity in a second order fluid to arise from the first normal stress
difference of the fluid. In $R_{vort}^{(1)}$, the FGP is an unstable limit
cycle and, in $R_{flow-vort}^{(1)}$, it is a stable limit cycle. In both
cases, the stable limit cycle exists near the vorticity axis. Thus, in
$R_{flow-vort}^{(1)}$, there is an unstable limit cycle between the two stable
ones, and there are two possible final orientation behaviors, tumbling within
the FGP and periodic motion close to the vorticity axis.

For $c\cdot De>cDe_{crit}=(4\log(2\kappa)-6)/\kappa$, the fixed points in the
FGP leads to monotonic particle drift toward or away from the FGP instead of
spiraling motion. For $\kappa>\kappa_{FV1}$ when the fixed points first appear
at $cDe_{crit}$ the one closer to the flow direction is a saddle node, and the
other is an unstable node. Thus a particle starting close to the FGP migrates
along the stable manifold of the saddle point up towards the flow direction
and then monotonically drifts along the flow-vorticity plane towards the
periodic orbit close to the vorticity axis. The region where this qualitative
behavior occurs is marked as $R_{vort}^{(2)}$ in figure 2 and the orientation
trajectories are similar to the experimental observations of Bartram et al.
(1975). Also, for $\kappa>\kappa_{FV1}$, upon increasing $c\cdot De$ beyond a
certain value indicated by the $\lambda_{2,\text{flow}}^{0-}=0$ boundary in
figure 2, the saddle point changes into a stable node and hence a particle
with an initial orientation close to FGP ultimately becomes stably flow
aligned. This is marked as $R_{flow-vort}^{(2)}$ in figure 2. In this region,
depending upon the initial orientation, a particle finally obtains either a
stable orientation close to the flow direction or undergoes periodic motion
close to the vorticity axis. For $c\cdot De$ beyond the
$\lambda_{2,\text{flow}}^{0--}=0$ boundary of figure 2 for
$\kappa>\kappa_{FV1}$, or for $c\cdot De>cDe_{\text{crit}}$ when
$\kappa<\kappa_{FV1}$, the fixed point in the FGP away from the flow direction
is a saddle node with its stable manifold perpendicular to the FGP. In the
part of this region marked as $R_{flow-vort}^{(3)}$ in figure 2, the final
particle behavior is similar to that in $R_{flow-vort}^{(2)}$ discussed above.
But, due to a saddle node instead of an unstable fixed point further from the
flow direction in the FGP, a greater proportion of the initial orientations in
$R_{flow-vort}^{(3)}$ than in $R_{flow-vort}^{(2)}$ lead the particle towards
a flow aligned state. Upon increasing $c\cdot De$ beyond $c\cdot De|_{cut-
off}^{Limit}$, in the regions marked as $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$,
the particle at all initial orientations apart from the log-rolling state at
the vorticity axis lead to a flow-aligned state. This is because the limit
cycle close to the vorticity axis does not exist and the vorticity axis is an
unstable spiral ($R_{flow}^{(1)}$) or node ($R_{flow}^{(2)}$). The flow-
alignment of the particle was observed in the experiments of Iso et al.
(1996a) at larger elasticity ($c$).

We are able to obtain qualitative agreement with the low shear rate
experiments of Gauthier et al. (1971) and Iso et al. (1996b) which correspond
to low $c\cdot De$ regime ($R_{vort}^{(1)}$ in figure 2). In these
experiments, the particle spirals towards the vorticity axis, and our theory
predicts spiraling towards the stable limit cycle near the vorticity axis.
While the authors (Gauthier et al., 1971; Iso et al., 1996b) state that the
particle approaches the vorticity axis, this behavior is not shown. Moreover,
in the time series plot of angle with the vorticity axis in the results of Iso
et al. (1996b) a deviation from zero and a small oscillatory behavior can be
observed, indicating the possibility of the existence of a stable limit cycle
near the vorticity axis in their experiments. The high shear rate experiments
of Bartram et al. (1975) were performed by releasing a $\kappa=9.1$ rod close
to the gradient axis. The particle initially travels nearly parallel to the
flow-gradient plane and approaches the flow direction. Upon perturbing, it
moves out of the flow-gradient plane along the flow-vorticity plane in a
monotonic fashion before it starts spiraling near the vorticity axis. This is
qualitatively similar to the orientation dynamics in the $R_{vort}^{(2)}$
region shown in figure 2 where a saddle point exists near the flow direction
(stable manifold in the flow-gradient plane) and a stable limit cycle is near
vorticity axis. Our theory is however unable to fully explain the high shear
rate experiments of Iso et al. (1996b) performed for $\kappa=34.4$, $De=3$ and
$c=0.39$. Iso et al. (1996b) observe the particle to remain at an angle
between 5 to 60∘ from the vorticity axis near the flow-vorticity plane after
spiraling or moving monotonically away from the FGP. The regions $R_{flow-
vort}^{(2)}$, $R_{flow-vort}^{(3)}$, $R_{flow}^{(1)}$, and $R_{flow}^{(2)}$
shown in figure 2) do have an attractor near the flow-vorticity plane, but it
is a saddle point. From equation (127) ${p}_{\text{3,saddle}}=0.7$ for $c\cdot
De=1.17$ and $\kappa=34.4$ i.e. the saddle point is at an angle of 25.2∘ from
vorticity axis. In our theory, the absence of a stable fixed point in the
flow-vorticity plane between the vorticity and flow directions may be due to
neglecting shear thinning, finite $c$ effects, finite polymer length, and
polymer entanglement. We use an Oldroyd-B equation to model the polymer
stress. Finite polymer length may be captured by the FENE-P model and polymer
entanglement likely at larger $c$ by the Giesekus model Bird et al. (1987).
Future numerical investigations may be used to test our theory and further
clarify the previous experimental findings.

Our findings have a major implication in using dilute (low $c$) polymeric
liquids to achieve desired properties such as strength and anisotropy of
products manufactured from dilute fiber-filled suspensions mentioned in the
introduction. At very low $c\cdot De$, all the fibers will eventually have
orientations close to the vorticity axis. Therefore, adding a small polymer
concentration to the fluid in roll-to-roll manufacturing can lead to low
resistance films with higher anisotropy and hence better quality. Since the
limit cycle in the low $c\cdot De$ regime becomes closer to the vorticity axis
as $c\cdot De$ increases, the anisotropy can be tuned by changing the shear
rate or polymer relaxation time (De is the product of shear rate and polymer
relaxation time). Even higher anisotropy can be obtained if very large $De$
can be achieved since at a large $c\cdot De$, the fibers with all initial
orientations ultimately align near the flow direction. For moderate values of
$c\cdot De$, the flow field that precedes a period of simple shear may pre-
orient the fibers somewhat closer to the flow-gradient plane or the vorticity
axis. This will determine whether the fibers lie in the basin of attraction of
attractor near the flow direction or the attractor near the vorticity axis.
The polymer stresses in the simple shear flow can then drive the particles to
a single final orientation leading to highly aligned fibers.

## Acknowledgement

This work was supported by NSF grants 1803156 and 2206851 and NASA grant
80NSSC23K0348.

Declaration of Interests: The authors report no conflict of interest.

## References

  * Abtahi & Elfring (2019) Abtahi, S Arman & Elfring, Gwynn J 2019 Jeffery orbits in shear-thinning fluids. Physics of Fluids 31 (10), 103106\.
  * Bartram et al. (1975) Bartram, E, Goldsmith, HL & Mason, SG 1975 Particle motions in non-newtonian media. Rheologica Acta 14 (9), 776–782.
  * Batchelor (1970) Batchelor, GK 1970 Slender-body theory for particles of arbitrary cross-section in stokes flow. Journal of Fluid Mechanics 44 (3), 419–440.
  * Bird et al. (1987) Bird, Robert Byron, Armstrong, Robert Calvin & Hassager, Ole 1987 Dynamics of polymeric liquids. vol. 1: Fluid mechanics .
  * Boger (1977) Boger, DV 1977 A highly elastic constant-viscosity fluid. Journal of Non-Newtonian Fluid Mechanics 3 (1), 87–91.
  * Breitenbach (2002) Breitenbach, Jörg 2002 Melt extrusion: from process to drug delivery technology. European journal of pharmaceutics and biopharmaceutics 54 (2), 107–117.
  * Brunn (1977) Brunn, P 1977 The slow motion of a rigid particle in a second-order fluid. Journal of Fluid Mechanics 82 (3), 529–547.
  * Chae & Kumar (2008) Chae, Han Gi & Kumar, Satish 2008 Making strong fibers. Science 319 (5865), 908–909.
  * Cox (1970) Cox, RG 1970 The motion of long slender bodies in a viscous fluid part 1. general theory. Journal of Fluid mechanics 44 (4), 791–810.
  * Cox (1971) Cox, RG 1971 The motion of long slender bodies in a viscous fluid. part 2. shear flow. Journal of Fluid Mechanics 45 (4), 625–657.
  * d’Avino et al. (2014) d’Avino, G, Hulsen, MA, Greco, F & Maffettone, PL 2014 Bistability and metabistability scenario in the dynamics of an ellipsoidal particle in a sheared viscoelastic fluid. Physical Review E 89 (4), 043006.
  * Férec et al. (2021) Férec, Julien, Bertevas, Erwan, Khoo, Boo Cheong, Ausias, Gilles & Phan-Thien, Nhan 2021 Rigid fiber motion in slightly non-newtonian viscoelastic fluids. Physics of Fluids 33 (10), 103320.
  * Gauthier et al. (1971) Gauthier, F, Goldsmith, HL & Mason, SG 1971 Particle motions in non-newtonian media. Rheologica Acta 10 (3), 344–364.
  * Guckenheimer & Holmes (2013) Guckenheimer, John & Holmes, Philip 2013 Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, , vol. 42. Springer Science & Business Media.
  * Gunes et al. (2008) Gunes, DZ, Scirocco, Rossella, Mewis, Joannes & Vermant, Jan 2008 Flow-induced orientation of non-spherical particles: Effect of aspect ratio and medium rheology. Journal of Non-Newtonian Fluid Mechanics 155 (1-2), 39–50.
  * Harlen & Koch (1993) Harlen, OG & Koch, Donald L 1993 Simple shear flow of a suspension of fibres in a dilute polymer solution at high deborah number. Journal of Fluid Mechanics 252, 187–207.
  * Huang et al. (2003) Huang, Zheng-Ming, Zhang, Y-Z, Kotaki, Masaya & Ramakrishna, Seeram 2003 A review on polymer nanofibers by electrospinning and their applications in nanocomposites. Composites science and technology 63 (15), 2223–2253.
  * Iso et al. (1996a) Iso, Yoichi, Cohen, Claude & Koch, Donald L 1996a Orientation in simple shear flow of semi-dilute fiber suspensions 2. highly elastic fluids. Journal of non-newtonian fluid mechanics 62 (2-3), 135–153.
  * Iso et al. (1996b) Iso, Yoichi, Koch, Donald L & Cohen, Claude 1996b Orientation in simple shear flow of semi-dilute fiber suspensions 1. weakly elastic fluids. Journal of non-newtonian fluid mechanics 62 (2-3), 115–134.
  * Jeffery (1922) Jeffery, George Barker 1922 The motion of ellipsoidal particles immersed in a viscous fluid. Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character 102 (715), 161–179.
  * Johnson et al. (1990) Johnson, Shirley J, Salem, Andrew J & Fuller, Gerald G 1990 Dynamics of colloidal particles in sheared, non-newtonian fluids. Journal of non-newtonian fluid mechanics 34 (1), 89–121.
  * Kim & Karrila (2013) Kim, Sangtae & Karrila, Seppo J 2013 Microhydrodynamics: principles and selected applications. Courier Corporation.
  * Leal (1975) Leal, LG 1975 The slow motion of slender rod-like particles in a second-order fluid. Journal of Fluid Mechanics 69 (2), 305–337.
  * Leal & Hinch (1971) Leal, LG & Hinch, EJ 1971 The effect of weak brownian rotations on particles in shear flow. Journal of Fluid Mechanics 46 (4), 685–703.
  * Magda et al. (1991) Magda, JJ, Lou, Jianzhong, Baek, SG & DeVries, KL 1991 Second normal stress difference of a boger fluid. Polymer 32 (11), 2000–2009.
  * Mutiso et al. (2013) Mutiso, Rose M, Sherrott, Michelle C, Rathmell, Aaron R, Wiley, Benjamin J & Winey, Karen I 2013 Integrating simulations and experiments to predict sheet resistance and optical transmittance in nanowire films for transparent conductors. ACS nano 7 (9), 7654–7663.
  * Nakajima et al. (1994) Nakajima, Toshinari, Kajiwara, Kanji & McIntyre, James Eric 1994 Advanced fiber spinning technology. Woodhead Publishing.
  * Phan-Thien & Fan (2002) Phan-Thien, Nhan & Fan, Xi-Jun 2002 Viscoelastic mobility problem using a boundary element method. Journal of non-newtonian fluid mechanics 105 (2-3), 131–152.
  * Stover & Cohen (1990) Stover, Carl A & Cohen, C 1990 The motion of rodlike particles in the pressure-driven flow between two flat plates. Rheologica Acta 29 (3), 192–203.
  * Strogatz (2018) Strogatz, Steven H 2018 Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. CRC press.
  * Wang et al. (2020) Wang, Shiyan, Tai, Cheng-Wei & Narsimhan, Vivek 2020 Dynamics of spheroids in an unbound quadratic flow of a general second-order fluid. Physics of Fluids 32 (11).
  * Wang et al. (2019) Wang, Yelong, Yu, Zhaosheng & Lin, Jianzhong 2019 Numerical simulations of the motion of ellipsoids in planar couette flow of giesekus viscoelastic fluids. Microfluidics and Nanofluidics 23, 1–16.
  * Yin et al. (2010) Yin, ZhouPing, Huang, YongAn, Bu, NingBin, Wang, XiaoMei & Xiong, YouLun 2010 Inkjet printing for flexible electronics: Materials, processes and equipments. Chinese Science Bulletin 55 (30), 3383–3407.