# Perfect stimulated Raman adiabatic passage with imperfect finite-time pulses

Shruti Dogra<EMAIL_ADDRESS>Gheorghe Sorin Paraoanu
<EMAIL_ADDRESS>QTF Center of Excellence, Department of Applied
Physics, Aalto University School of Science, P.O. Box 15100, FI-00076 AALTO,
Finland

###### Abstract

We present a well-tailored sequence of two Gaussian-pulsed drives that
achieves perfect population transfer in STImulated Raman Adiabatic Passage
(STIRAP). We give a theoretical analysis of the optimal truncation and
relative placement of the Stokes and pump pulses. Further, we obtain the power
and the duration of the protocol for a given pulse width. Importantly, the
duration of the protocol required to attain a desired value of fidelity
depends only logarithmically on the infidelity. Subject to optimal truncation
of the drives and with reference to the point of fastest transfer, we obtain a
new adiabaticity criteria, which is remarkably simple and effective.

††preprint: APS/123-QED

##  I Introduction

The STIRAP (Stimulated Raman Adiabatic Passage) protocol got its first
validation in an experiment [1, 2] where partially overlapping Stokes and pump
laser beams were employed to transfer the population from a lower energy state
to a higher vibrational state without populating the intermediate level in a
three-level system consisting of molecular vibrational states. This was done
with a non-trivial pulse arrangement (usually referred to as counter-intuitive
sequence), where the Stokes pulse precedes the pump pulse. This selective and
precise adiabatic transfer of population has been a subject of much interest
from a theoretical as well as experimental perspective [3, 4]. Due to its
intrinsic robustness against practical imperfections, STIRAP has been widely
adopted in various different experimental systems [5, 6]. Most importantly,
the success of the protocol (even from the theoretical point of view) relies
on the fulfilment of the adiabaticity criteria [7, 8, 9, 10, 11, 12]. As per
quantum adiabatic theorem [13, 14], the system, which is initialised in an
eigenstate, follows the corresponding eigenstate of the instantaneous
Hamiltonian. However, a widely acceptable quantitative criteria for
adiabaticity is still lacking [15]. An interesting approach based on local
adiabaticity criteria is discussed in [16], where the Hamiltonian generating
adiabatic evolution is designed in a such a way that it fulfils the local
adiabaticity condition at infinitesimal time steps, which is further used to
obtain the adiabatic-evolution version of the Grover’s search algorithm.

Here we present an adiabaticity criteria which is demonstrably sufficient for
achieving perfect population transfer. Our criteria is markedly different from
the existing ones and is surprisingly effective despite its simplicity. The
key concept of our analysis is based on the most sensitive point of the
dynamics, which is in the middle of the sequence, where the rate of evolution
of the quantum state is the highest. Also, a high-fidelity STIRAP requires the
pulse sequence to be implemented in an optimal time, which involves optimal
truncation of the drives, as well as optimal width of the drives and relative
placement of the drives in the pulse sequence. Here we show how these issues
can be solved. Another important aspect is the power of the pulses, which we
obtain optimally with the help of our newly introduced adiabaticity criteria.
We analyse the situation in detail and arrive at analytical expressions that
lead to a perfectly tailored STIRAP.

The protocol studied here is experimentally implementable to any three-level
system. The set of parameters presented here can be directly used in a circuit
QED based experimental setup with a multi-level Josephson-junction artificial
atom [17]. There are also ways to suppress the non-adiabatic excitations by
employing shortcuts to adiabaticity [6]. In superadiabatic(sa)-STIRAP in
three-level systems, an additional counterdiabatic pulse is needed, realizing
direct coupling between the initial and the target states. A circuit QED based
setup implementing saSTIRAP protocol in a three-level system has been
demonstrated in Ref. [18] and its robustness against various experimental
imperfections has been analysed in Ref. [19].

The Hamiltonian governing the STIRAP for a three-level system in the
computational basis $\\{|0\rangle$, $|1\rangle$, $|2\rangle\\}$, in the
dispersive regime and under the rotating wave approximation, is given by [17]

$H_{0}=\frac{\hbar}{2}\left(\begin{array}[]{ccc}0&\Omega_{01}(t)&0\\\
\Omega_{01}(t)&2\delta_{01}&\Omega_{12}(t)\\\
0&\Omega_{12}(t)&2(\delta_{01}+\delta_{12})\end{array}\right),$ (1)

where the time-varying amplitudes of the driving fields are chosen as
Gaussians with equal standard deviation $\sigma$. These Gaussians are
separated in time by an amount $t_{s}$ given by

$\Omega_{01}(t)=\Omega_{01}^{0}e^{-t^{2}/2\sigma^{2}}\quad\textrm{and}\quad\Omega_{12}(t)=\Omega_{12}^{0}e^{-(t-t_{s})^{2}/2\sigma^{2}}.$

An adiabatic evolution is ideally infinitely slow and would require the system
to be in an eigenstate of the instantaneous Hamiltonian at all times. At the
two-photon resonance condition (i.e. $\delta_{01}=-\delta_{12}$), a convenient
choice of the eigenvector is the dark state $|{\rm
D}\rangle=\cos\Theta|0\rangle-\sin\Theta|2\rangle$, which does not have any
dependence on the intermediate level $|1\rangle$. Here the mixing angle
$\Theta$ is defined by $\Theta=\tan^{-1}\Omega_{01}(t)/\Omega_{12}(t)$.

##  II Optimal pulse duration

Adiabatic drive in principle demands infinitely long operation time for a
complete transfer of population. Ideally, as required by STIRAP, Gaussians
pulses are of infinite extent. However, to cope with the experimental
limitations on pulse generation and to minimize the losses due to decoherence,
one would have to truncate the Gaussians $\Omega_{01}(t)$ and $\Omega_{12}(t)$
optimally. Therefore there is a tradeoff between the loss in the transfer
fidelity that can be afforded and the total pulse time. Revisiting the mixing
angle, while assuming $\Omega_{01}^{0}=\Omega_{12}^{0},$ we write

$\tan\Theta(t)=\frac{\Omega_{01}(t)}{\Omega_{12}(t)}=\frac{e^{-rt/\sigma}}{e^{-r^{2}/2}},$
(2)

where we introduce the parameter $r=t_{s}/\sigma$.

We truncate this STIRAP pulse sequence (consisting of drives $\Omega_{01}(t)$
and $\Omega_{12}(t)$) from left at time
$t=t_{i}=-n_{t}\sigma+t_{s}=-(n_{t}-r)\sigma$, which we call initial time
point, and from right at $t=t_{f}=n_{t}\sigma$, which we call final time
point, where $n_{t}$ is a real number ($n_{t}\in\mathbb{R}$). The total pulse
duration is therefore $T=t_{f}-t_{i}=(2n_{t}-r)\sigma$. We fix the values of
$r$ and $\sigma$ and present the corresponding dynamics of $\Theta(t)$ versus
total pulse duration as shown in Fig. 1, where different curves correspond to
different values of $n_{t}$. The width of the Gaussian may be fixed to any
arbitrary value (here $\sigma=30$ ns) as this does not effect the variation of
$\Theta(t)$ in a given total time $T$. Ideally, during the STIRAP drive, the
mixing angle $\Theta(t)$ is expected to vary from $0$ to $\pi/2$, while in
reality, a finite-time sequence effectively varies $\Theta(t)$ from
$\Theta_{i}\rightarrow 0$ to $\Theta_{f}=\pi/2-\Theta_{i}$. A closer look at
Fig. 1 immediately concludes that a choice of small enough $n_{t}$ might
result in a large $\Theta_{i}$ _e.g._ in Fig. 1(b) the blue curve marked with
circles has $n_{t}=0$ and $\Theta_{i}>\pi/6$ while the orange curve marked
with squares corresponds to $n_{t}=4$ and $\Theta_{i}\approx 0$. Thus, for a
given $r$, lower values of $n_{t}$ result in a poor transfer fidelity. To make
it worse, real situations by default have $\Theta_{i}\approx 0$, therefore a
too small value of $n_{t}$ is susceptible to create errors which can be
difficult to trace. An elaborated picture of the ideal situation is presented
in Fig. 1(b), where curves corresponding to $n_{t}=1,4$ have
$\Theta_{i}=12.6^{0}$ and $0.64^{0}$ respectively, with total pulse duration
being $105$ ns and $315$ ns respectively. Clearly in this case $n_{t}=4$ has a
much more desirable outcome than that of $n_{t}=1$, despite the high time
cost. Another important factor that plays a role in the time management of the
STIRAP implementation is the relative separation between the two pulses
($r=t_{s}/\sigma$). Comparing curves corresponding to $n_{t}=1$ in Fig. 1(a)
and (b), it is found that for $r=-2$, $\Theta_{i}\approx 1^{0}$ and the total
pulse duration is $140$ ns. Thus an optimal combination of $n_{t}$ and $r$
provides an efficient STIRAP without compromising much with respect to the
time cost. A quite thorough picture can be obtained from the contour plot in
Fig. 2, wherein mixing angle corresponding to the final state
($\Theta_{f}=\pi/2-\Theta_{i}$) is plotted for different combinations of
$n_{t}$ and $r$. The final value of the mixing angle ranges from
$\Theta_{f}=45^{0}$ (corresponding to $r=0,n_{t}=0$ and thus no time
evolution) to $\Theta_{f}=90^{0}$, which corresponds to a complete transfer of
population from $|0\rangle\rightarrow|2\rangle$.

Figure 1: Time variation of $\Theta(t)$ versus total pulse duration $T$ for
different values of $n_{t}$ is shown at two different values of $r$: (a)
$r=-2$ and (b) $r=-1$. In each of these figures, blue curve with circular
markers, red curve with triangles, black curve with diamonds, and orange curve
with squares correspond to $n_{t}=0,1,2,4$ respectively (also specified at the
right end of each curve in part (b)).  Figure 2: Contour plot showing the
mixing angle at time $t=t_{f}$ corresponding to the normalized relative
separation $r$ between the two driving fields and the parameter $n_{t}$
characterizing the transfer time.

Next, we quantify the threshold for a desired fidelity of the STIRAP protocol.
From Eq. (2),

$\displaystyle-\frac{rt}{\sigma}$ $\displaystyle=$
$\displaystyle-\frac{r^{2}}{2}+\ln\left[\tan\Theta(t)\right].$ (3)

As stated earlier, due to the STIRAP driven dynamics, our three-level system
is in the dark state $|{\rm D}\rangle$ at all times. We parametrize the
initial dark state with close proximity to the ground state $|0\rangle$ by

$|\psi_{i}\rangle=\left(\begin{array}[]{c}\cos\Theta_{i}\\\ 0\\\
-\sin\Theta_{i}\end{array}\right)=\left(\begin{array}[]{c}\sqrt{1-\epsilon^{2}}\\\
0\\\ -\epsilon\end{array}\right),$ (4)

such that the ideal case is recovered when $\epsilon\longrightarrow 0$. Since
we know that during the evolution under the STIRAP Hamiltonian our system is
in the dark state at all times, the final state reads

$|\psi_{f}\rangle=\left(\begin{array}[]{c}\cos\Theta_{f}\\\ 0\\\
-\sin\Theta_{f}\end{array}\right)=\left(\begin{array}[]{c}\epsilon\\\ 0\\\
-\sqrt{1-\epsilon^{2}}\end{array}\right).$ (5)

For both the initial and final states, the parameter $\epsilon$ is a measure
of infidelity. Indeed, the fidelity is $F_{i}=F_{f}=|\langle
0|\psi_{i}\rangle|^{2}=|\langle 2|\psi_{f}\rangle|^{2}=1-\epsilon^{2}$ and
therefore the infidelity is $\epsilon^{2}$. At the final time point,
$t=t_{f}=n_{t}\sigma$, and from Eqs. (3), (4), and (5), one may easily arrive
at

$n_{t}=\frac{1}{r}\ln\frac{\epsilon}{\sqrt{1-\epsilon^{2}}}+\frac{r}{2}.$ (6)

Thus the total pulse duration is

$T=\frac{2\sigma}{r}\ln\frac{\epsilon}{\sqrt{1-\epsilon^{2}}},$ (7)

which, as expected, is directly proportional to the widths of the Gaussians. A
larger value of $|r|$ corresponds to faster truncation (smaller $n_{t}$) and
is overall advantageous in terms of the total pulse duration. For small enough
$\epsilon$ ($\epsilon^{2}<<1$), Eq. 7 leads to $\epsilon\propto
e^{rT/(2\sigma)}$. Thus, the infidelity decreases exponentially with total
time.

Fig. 3(a) contains plots of $n_{t}$ vs $-r$, where different curves correspond
to different values of $\epsilon$. It is interesting to note that there may
exist the STIRAP sequence even for negative values of $n_{t}$ when the
relative separation between the two Gaussians is large. However, the total
transfer time is positive as expected, which is clearly seen in Fig. 3(b)
showing the variation of $T/\sigma$ with $-r$ at corresponding values of
$\epsilon$. Consider the vertical green line at $r=-1.5$ in Fig. 3(b), and the
values of $T/\sigma$ while it intersects different curves plotted at different
values of $\epsilon$. The smaller $\epsilon$ is, the higher the fidelity,
which requires larger values of $T/\sigma$ for a fixed value of $r$. For
instance, assuming $\sigma=40$ ns, a total transfer time of $184.2$ ns,
$122.8$ ns, $79.9$ ns, $61.3$ ns, and $14.6$ ns are required to obtain final
values of the mixing angle $\Theta_{f}$ to be $89.9^{0}$, $89.2^{0}$,
$85.9^{0}$, $81.9^{0}$, and $48.6^{0}$ respectively.

Figure 3: Plots of (a) $n_{t}$ vs $-r$ and (b) $T/\sigma$ vs $-r$
corresponding to different values of $\epsilon$.

Another interesting situation arises when
$\Omega_{01}^{0}\neq\Omega_{12}^{0}$, which influences the left and right
truncation limits, however the total transfer time remains unchanged. This
situation is discussed in detail in Appendix B.

##  III Adiabaticity criteria

Next, we evaluate the optimal value of the pulse amplitude corresponding to
the optimal transfer time calculated in the last section. The total pulse area
then may be compared with the total energy required to achieve the selective
population transfer. The adiabatic criteria for a STIRAP implementation
implies that at any arbitrary time $t$ the effective area $A(t)$ is much
greater than the time rate of change of the mixing angle,

$A(t)=\sqrt{\Omega_{01}^{2}(t)+\Omega_{12}^{2}(t)}\quad>>\quad\dot{\Theta}(t),$
(8)

which upon integration, gives rise to the global adiabaticity condition, as
discussed in [17].

Figure 4: Plots of the time-varying amplitudes $\Omega_{01}(t)/(2\pi)$ and
$\Omega_{12}(t)/(2\pi)$, effective area $A(t)/(2\pi)$, and $\dot{\Theta}(t)$
versus time are shown for (a) $r=-1.5$ and (b) $r=-2$. Here $\sigma=35$ns,
$\Omega_{01}^{0}/(2\pi)=44$ MHz, and $\Omega_{12}^{0}/(2\pi)=37$ MHz.

For a deeper insight into the protocol, let us look at the time dependence of
these quantities. The plots of the time-varying amplitudes of the driving
fields ($\Omega_{01}(t)$ and $\Omega_{12}(t)$), effective area $A(t)$, and
rate of change of the mixing angle $\dot{\Theta}(t)$ are shown in Fig. 4,
where $\sigma=35$ ns, $\Omega_{01}^{0}/(2\pi)=44$ MHz, and
$\Omega_{12}^{0}/(2\pi)=37$ MHz. It is easy to notice that the rate of change
of mixing angle is maximum in the middle of the sequence where
$\Omega_{01}(t)=\Omega_{12}(t)$.

Ideally (for $\epsilon\rightarrow 0$), at $t=t_{I}$, the mixing angle is
$\Theta=\pi/4$ and the populations are $p_{0}=p_{2}=0.5$. Also the STIRAP
sequence is fastest and more prone to errors at this point, such that a non-
zero $p_{1}$ occurs close to this point. Thus $t=t_{I}$ plays an important
role. It is also clearly seen in Fig. 4 that the variation of the mixing angle
with time attains its maximum value at $t=t_{I}$, when
$\Omega_{01}(t)=\Omega_{12}(t)$.

Thus, in most of the cases, especially the ones corresponding to poor
performance of STIRAP, the time $t=t_{I}$ corresponds to close values of the
terms on the left and right hand sides of the inequality in Eq. (8). In Fig.
4, at $t=t_{I}$, a vertical blue line intersects various curves, such that
$A(t_{I})/(2\pi)$, $\Omega_{01}(t_{I})/(2\pi)=\Omega_{12}(t_{I})/(2\pi)$,
$\dot{\Theta}(t_{I})$ and the time axis are labelled by points A, O, Q, and P
respectively. An intuitive argument based on observation leads to a non-
trivial relation, that has to be obeyed for a better performance of the
STIRAP, given by,

$\Omega_{01}(t_{I})-2\pi\dot{\Theta}(t_{I})\quad\gtrsim\quad\sqrt{\Omega_{01}^{2}(t_{I})+\Omega_{12}^{2}(t_{I})}-\Omega_{01}(t_{I}),$
(9)

or alternatively,

$\frac{\Omega_{01}(t_{I})-2\pi\dot{\Theta}(t_{I})}{\Omega_{01}(t_{I})}\quad\gtrsim\quad(\sqrt{2}-1),$
(10)

as $\Omega_{01}(t)=\Omega_{12}(t)$ at $t=t_{I}$. Time ($t_{I}$) is given by,

$t_{I}=\frac{t_{s}}{2}+\frac{\sigma^{2}}{t_{s}}\ln\alpha,$ (11)

where $\alpha=\Omega_{01}^{0}/\Omega_{12}^{0}$. Further,

$\displaystyle\dot{\Theta}(t)$ $\displaystyle=$
$\displaystyle\frac{\dot{\Omega}_{01}(t)\Omega_{12}(t)-\Omega_{01}(t)\dot{\Omega}_{12}(t)}{\Omega_{01}^{2}(t)+\Omega_{12}^{2}(t)}$

and at $t=t_{I}$,

$\displaystyle\quad\dot{\Theta}(t_{I})$ $\displaystyle=$
$\displaystyle-\frac{t_{s}}{2\sigma^{2}}.$ (13)

Note that even for $\alpha\neq 1$, $\dot{\Theta(t)}$, attains its maximum
value at $t=t_{I}$, where $\Omega_{01}(t)=\Omega_{12}(t)$. Irrespective of the
asymmetry introduced by different pulse amplitudes ($\Omega_{01}^{0}$ and
$\Omega_{12}^{0}$), the rate of population transfer reaches its maximum at
$t=t_{I}$. From Eqs. 9– 13, the condition for a better STIRAP result is given
by,

$\sigma\frac{\Omega_{01}^{0}}{2\pi}\gtrsim\frac{-re^{t_{I}^{2}/2\sigma^{2}}}{2(2-\sqrt{2})}.$
(14)

For $\Omega_{01}^{0}/(2\pi)=\Omega_{12}^{0}/(2\pi)=\Omega$, $\alpha=1$, and
the above inequality results into

$\sigma\Omega\gtrsim\frac{-re^{r^{2}/8}}{2(2-\sqrt{2})}.$ (15)

This is the adiabaticity criteria for the STIRAP population transfer, which is
obtained by assuming the system to be in the dark state at all times (see Eq.
4). The STIRAP population transfer calculated using Eqs. 6 and 15 (for
$\alpha=1$) will be labelled in the following as parameter ‘Set 1’. The
dependence of the right side of Eq. 15 with respect to $-r$ is plotted as
shown in Fig. 5(a) with continuous black curve with markers. The corresponding
population transfer obtained from Set 1 (for $\epsilon=0.05$) is shown in Fig.
5(b) with black markers. It is noteworthy that the plot of $p_{2}$ versus $-r$
(see Fig.5(b)) is independent of the values of $\sigma$ and is in fact
dependent on $\epsilon$. On the other hand, the total pulse duration T is
directly proportional to $\sigma$ and the drive amplitude $\Omega$ is
inversely proportional to $\sigma$. Furthermore, a larger value of $\epsilon$
leads to less efficient population transfer in significantly shorter time.

Figure 5: Plots for (a) STIRAP adiabaticity criteria ($\sigma\Omega=R.H.S.$)
and (b) corresponding final population in the second excited state ($p_{2}$)
for Set 1 (Eqs. (6) and (15)) with continuous black curve with markers, for
Set 2 (Eqs. 16 and 17) diamond markers and the global adiabaticity condition
with continuous blue line are shown. The inset shows an elaborated view of the
respective quantities for $|r|<1$. (c) Total pulse duration in blue with
square markers and drive amplitude in red with circular markers are shown for
Set 2, where $\epsilon=0.05$, $\alpha=1$, and $\sigma=30$ ns.

The results of Eq. (15) are compared with the global adiabaticity criteria,
$\sigma\Omega>>\sqrt{\pi}/4$ [17]. In Fig 5(a), the continuous blue line
corresponds to $\sigma\Omega=\sqrt{\pi}/4$, which means that for the global
adiabaticity criteria to be satisfied, the product $\sigma\Omega$ must lie
significantly above the blue line. On comparison between the continuous blue
curve and black curve with markers in Fig 5(a), we find that it is possible to
have an effective population transfer even when the global adiabaticity
criteria is clearly violated in the region $|r|<0.5$. For instance,
considering a STIRAP evolution with $r=-0.4$, $\sigma=5$ ns, $\Omega=70$ MHz,
and optimally tailored time of $115$ ns (where $n_{i}=n_{f}=11.3$ with
$\epsilon=0.01$), starting from the ground state yields a final state
population, $p_{2}\approx 0.99$ that goes beyond the global adiabaticity
condition (stated earlier), and $T\Omega=8.02$, which is $<10$ and hence
violates the adiabaticity condition reported in [20, 21]. Also, for larger
values of $|r|$, Fig 5(a) presents much disparity between the two adiabaticity
criteria.

Further, on close observation of the $p_{2}$ population profile in Fig 5(b),
we find that the population transfer is not very efficient for $|r|<1$, and
that the above example of a perfect transfer at $r=-0.4$ is a mere
coincidence. These imperfections originate in the assumption of dark state
dynamics, which is not valid for small values of $|r|$ due to spurious
excitations. They can be compensated by higher power of the drive. Based on
these phenomenological considerations, we can design an optimal set of
parameters for $|r|<1$. We call this $r-$conditioned set of parameters and
label it as Set 2, given by

$\displaystyle n_{t}$ $\displaystyle=$
$\displaystyle-\ln\frac{\epsilon}{\sqrt{1-\epsilon^{2}}}+\frac{r}{2},\quad|r|<1$
(16) $\displaystyle=$
$\displaystyle\frac{1}{r}\ln\frac{\epsilon}{\sqrt{1-\epsilon^{2}}}+\frac{r}{2},\quad|r|\geq
1$ $\displaystyle\sigma\Omega$ $\displaystyle\gtrsim$
$\displaystyle\frac{e^{r^{2}/8}}{2(2-\sqrt{2})},\quad|r|<1$ (17)
$\displaystyle\gtrsim$
$\displaystyle\frac{-re^{r^{2}/8}}{2(2-\sqrt{2})},\quad|r|\geq 1,$

where $\alpha=1$. The plot for the right hand side of the inequality in Eq.
(17) versus $r$ is shown as diamond markers in Fig 5(a) and the corresponding
final population $p_{2}$, calculated from Set 2 is shown in 5(b). Clearly, Fig
5(b) reflects a perfect STIRAP transfer for parameter Set 2. Fig. 5(c) shows
the variation of total pulse duration (T) and the drive amplitude ($\Omega$)
as functions of $r$, where $\sigma=30$ ns.

We also simulated STIRAP where the product $\sigma\Omega$ is close to but
lesser than the respective right hand sides of Eq. (17), and find that we can
still arrive at a good enough population transfer in certain situations. This
is especially true for large values of $|r|$. Thus, we conclude that the
adiabaticity condition in Eq. (17) is sufficient but not necessary for a
perfect population transfer.

##  IV A perfect STIRAP protocol

A demonstration of the improvement achieved by employing the conditions in
Eqs. (16) and (17) is shown in Fig. 6, where population $p_{2}$ of the second
excited state ($|2\rangle$) at the final time $t=t_{f}$ is plotted as a
function of $\sigma$ and $-r$, with $\alpha=1$ and $\epsilon=0.05$. Fig. 6(a)
shows the $p_{2}$ with fixed $n_{t}=3$ and $\Omega=45$ MHz. Thus, any
arbitrary point on the $p_{2}$-map in Fig. 6(a) satisfies
$\sigma\Omega>\sqrt{\pi}/4$. The simulation of a perfectly tailored STIRAP
utilizing the conditions in Eqs. (16) and (17) is shown in Fig. 6(b). The
resultant population profile demonstrates a well-tailored STIRAP protocol for
the desired population transfer. Fig. 6 (b) shows clear improvement relative
to results shown in Fig. 6 (a).

For a practical implementation, the total time-cost and pulse-power evaluation
are also important. The corresponding maps of the total pulse duration and
maximum pulse-amplitude in the same ranges of $\sigma$ and $-r$ are shown in
the Appendix in Fig. 7. The wide range of resultant $T$ and $\Omega$ values
provide flexibility to the protocol. Larger values of $\epsilon$ further lead
to significant reduction in the time-cost. In turn, a choice of slightly
larger $\sigma$ can significantly reduce the amplitude $\Omega$.

Figure 6: Population of the second excited state ($p_{2}$) is plotted at time
$t=t_{f}$ in a STIRAP protocol, with initial state $|0\rangle$. Results from
(a) STIRAP with $\Omega=45$ MHz and $n_{t}=3$ and (b) well-tailored STIRAP as
per Eqs. (16) and (17) are shown.

Alternatively, when evaluating the experimental feasibility various parameters
can be constrained and a perfect population transfer can be designed with the
help of these interweaved parametric equations and graphs.

For a desired value of $\epsilon$, and with the help of Fig. 3, Eqs. (6)or
(7), and Eq. (15), one can easily obtain an experimentally feasible set of
parameters $r$, $\sigma$, $n_{t}$, and $\Omega$ that leads to a perfect
STIRAP. It is noteworthy that the efficacy of this perfect STIRAP protocol
does not rigidly rely on the calculated parameters. In fact, parameters such
as $n_{t}$, $\Omega$ can be considered as the respective lower bounds to
achieve population transfer with infidelity $\epsilon^{2}$. Larger values of
these parameters will only make the transfer more efficient. This makes the
protocol robust against the experimental imperfections.

##  V Discussion and conclusions

We presented a well-tailored STIRAP protocol that leads to a perfect
population transfer $|0\rangle\rightarrow-|2\rangle$ alongside with
flexibility in the choice of parameters. For a given $\epsilon$, a combined
choice of parameters $n_{t}$ and $r$ already determines the final population
to be transferred. Furthermore, the choice of $\sigma$ determines the total
pulse duration ($T$) and the corresponding calculated value of the amplitude
$\Omega$ is responsible for the pulse power. A trade off between $\sigma$ and
$\Omega$ values can be settled by evaluating the experimental feasibility. We
also discussed the relatively general situation, where the Gaussian drives can
have unequal maximum amplitudes controlled by the parameter
$\alpha=\Omega_{01}^{0}/\Omega_{12}^{0}$. The results for $\alpha\neq 1$ are
discussed in the main text, while a detailed analysis is presented in the
appendices. The analysis presented here relies on a simple set of calculations
and observations, however the end results are non-trivial. In conclusion, our
calculations for the STIRAP drives lead to a perfect population transfer
within the reach of experimentally feasible scenario and without the help of
any additional shortcuts to the adiabaticity.

###### Acknowledgements.

We acknowledge financial support from the Academy of Finland under the Finnish
Center of Excellence in Quantum Technology QTF (projects 312296, 336810) and
from Grant No. FQXi-IAF19-06 (“Exploring the fundamental limits set by
thermodynamics in the quantum regime”) of the Foundational Questions Institute
Fund (FQXi), a donor advised fund of the Silicon Valley Community Foundation.

##  Appendix A Practical feasibility

Figure 7: Variation of (a) total pulse duration (T) and (b) maximum pulse-
amplitude ($\Omega$) as a function of $\sigma$ and $-r$.

The cost in terms of the total pulse duration and maximum pulse-amplitude is
shown in Fig. 7 as surface maps for a wide range of $\sigma$ and $-r$,
consistent with Fig. 6(b) of the main text.

##  Appendix B Gaussians with different amplitudes

Consider the pair of Gaussians together leading to a STIRAP pulse sequence for
$\Omega_{01}^{0}\neq\Omega_{12}^{0}$ and $\sigma_{1}=\sigma_{2}=\sigma$,

$\displaystyle\Omega_{01}(t)$ $\displaystyle=$
$\displaystyle\Omega_{01}^{0}e^{-t^{2}/2\sigma^{2}},$
$\displaystyle\Omega_{12}(t)$ $\displaystyle=$
$\displaystyle\Omega_{12}^{0}e^{-(t-t_{s})^{2}/2\sigma^{2}}.$ (18)

The mixing angle is given by

$\displaystyle\tan\Theta(t)$ $\displaystyle=$
$\displaystyle\frac{\Omega_{01}(t)}{\Omega_{12}(t)}=\frac{\Omega_{01}^{0}}{\Omega_{12}^{0}}e^{-rt/\sigma}e^{r^{2}/2},$
(19) $\displaystyle\textrm{or}\quad\frac{rt}{\sigma}$ $\displaystyle=$
$\displaystyle\frac{r^{2}}{2}-\ln\left[\tan\Theta(t)\right]+\ln\left[\Omega_{01}^{0}/\Omega_{12}^{0}\right].$
(20)

where, $r=t_{s}/\sigma$. For a finite time operation, this sequence of
Gaussians is truncated at the optimal time that provides the complete transfer
of population. The initial and final time points of the sequence are obtained
by truncating the sequence from left at $t=t_{i}=-n_{i}\sigma+t_{s}$ and from
right at $t=t_{f}=n_{f}\sigma$, such that the total operation time is
$T=\sigma(n_{i}+n_{f})-t_{s}$. Interestingly, we come across the same mixing
angle ($\Theta$) in the structure of the dark state $|\psi_{D}\rangle$, which
is an eigenstate of the instantaneous Hamiltonian with eigenvalue $0$. We
assume our three-level quantum system in the dark state with close proximity
to the ground state ($|0\rangle$),

$|\psi_{i}\rangle=\left(\begin{array}[]{c}\sqrt{1-\epsilon^{2}}\\\ 0\\\
-\epsilon\end{array}\right),$ (21)

such that $\epsilon\longrightarrow 0$. We know that during the evolution under
the STIRAP Hamiltonian, our system is ideally in the dark state at all times,
such that the final state is

$|\psi_{f}\rangle=\left(\begin{array}[]{c}\epsilon\\\ 0\\\
-\sqrt{1-\epsilon^{2}}\end{array}\right).$ (22)

At the final time point, $t=t_{f}=n_{f}\sigma$, and from Eqs. (4), (5), and
(20), one may easily arrive at

$n_{f}=\frac{1}{r}\ln\frac{\epsilon}{\sqrt{1-\epsilon^{2}}}+\frac{r}{2}+\frac{1}{r}\ln\frac{\Omega_{01}^{0}}{\Omega_{12}^{0}}.$
(23)

Similarly, at $t=t_{i}=-n_{i}\sigma+t_{s}$

$n_{i}=\frac{1}{r}\ln\frac{\epsilon}{\sqrt{1-\epsilon^{2}}}+\frac{r}{2}-\frac{1}{r}\ln\frac{\Omega_{01}^{0}}{\Omega_{12}^{0}},$
(24)

and the total time,

$T=(n_{i}+n_{f})\sigma-
t_{s}=(n_{i}+n_{f}-r)\sigma=\frac{2\sigma}{r}\ln\frac{\epsilon}{\sqrt{1-\epsilon^{2}}}.$
(25)

Consistent with parameter Set 2, the above equations are valid for $|r|\geq
1$. For $|r|<1$ we have,

$\displaystyle n_{f}$ $\displaystyle=$
$\displaystyle-\ln\frac{\epsilon}{\sqrt{1-\epsilon^{2}}}+\frac{r}{2}+\frac{1}{r}\ln\frac{\Omega_{01}^{0}}{\Omega_{12}^{0}}.$
(26) $\displaystyle n_{i}$ $\displaystyle=$
$\displaystyle-\ln\frac{\epsilon}{\sqrt{1-\epsilon^{2}}}+\frac{r}{2}-\frac{1}{r}\ln\frac{\Omega_{01}^{0}}{\Omega_{12}^{0}},$
(27) $\displaystyle T$ $\displaystyle=$ $\displaystyle(n_{i}+n_{f})\sigma-
t_{s}=(n_{i}+n_{f}-r)\sigma$ (28) $\displaystyle=$ $\displaystyle
2\sigma\ln\frac{\epsilon}{\sqrt{1-\epsilon^{2}}}.$

When $\Omega_{01}^{0}\neq\Omega_{12}^{0}$, truncation from the left and right
extremes correspond to slightly different values of $n_{i}$ and $n_{f}$ to
obtain the final state with infidelity $\epsilon^{2}$. This is due to the same
amount of fractional decrease in the amplitudes of the pulses expected at the
initial and final time points.

## References

  * Gaubatz  _et al._ [1988] U. Gaubatz, P. Rudecki, M. Becker, S. Schiemann, M. Külz, and K. Bergmann, Population switching between vibrational levels in molecular beams, Chemical Physics Letters 149, 463 (1988).
  * Gaubatz  _et al._ [1990] U. Gaubatz, P. Rudecki, S. Schiemann, and K. Bergmann, Population transfer between molecular vibrational levels by stimulated raman scattering with partially overlapping laser fields. a new concept and experimental results, The Journal of Chemical Physics 92, 5363 (1990), https://doi.org/10.1063/1.458514 .
  * Shore [1995] B. W. Shore, Examples of counter-intuitive physics, Contemporary Physics 36, 15 (1995), https://doi.org/10.1080/00107519508222135 .
  * Shore [2017] B. W. Shore, Picturing stimulated raman adiabatic passage: a stirap tutorial, Adv. Opt. Photon. 9, 563 (2017).
  * Vitanov  _et al._ [2017] N. V. Vitanov, A. A. Rangelov, B. W. Shore, and K. Bergmann, Stimulated raman adiabatic passage in physics, chemistry, and beyond, Rev. Mod. Phys. 89, 015006 (2017).
  * Bergmann  _et al._ [2019] K. Bergmann, H.-C. Nägerl, C. Panda, G. Gabrielse, E. Miloglyadov, M. Quack, G. Seyfang, G. Wichmann, S. Ospelkaus, A. Kuhn, S. Longhi, A. Szameit, P. Pirro, B. Hillebrands, X.-F. Zhu, J. Zhu, M. Drewsen, W. K. Hensinger, S. Weidt, T. Halfmann, H.-L. Wang, G. S. Paraoanu, N. V. Vitanov, J. Mompart, T. Busch, T. J. Barnum, D. D. Grimes, R. W. Field, M. G. Raizen, E. Narevicius, M. Auzinsh, D. Budker, A. Pálffy, and C. H. Keitel, Roadmap on STIRAP applications, J. Phys. B 52, 202001 (2019).
  * Kuklinski  _et al._ [1989] J. R. Kuklinski, U. Gaubatz, F. T. Hioe, and K. Bergmann, Adiabatic population transfer in a three-level system driven by delayed laser pulses, Phys. Rev. A 40, 6741 (1989).
  * Benseny and Mølmer [2021] A. Benseny and K. Mølmer, Adiabatic theorem revisited: The unexpectedly good performance of adiabatic passage, Phys. Rev. A 103, 062215 (2021).
  * Comparat [2009] D. Comparat, General conditions for quantum adiabatic evolution, Phys. Rev. A 80, 012106 (2009).
  * Tong  _et al._ [2007] D. M. Tong, K. Singh, L. C. Kwek, and C. H. Oh, Sufficiency criterion for the validity of the adiabatic approximation, Phys. Rev. Lett. 98, 150402 (2007).
  * Amin [2009] M. H. S. Amin, Consistency of the adiabatic theorem, Phys. Rev. Lett. 102, 220401 (2009).
  * Tong  _et al._ [2005] D. M. Tong, K. Singh, L. C. Kwek, and C. H. Oh, Quantitative conditions do not guarantee the validity of the adiabatic approximation, Phys. Rev. Lett. 95, 110407 (2005).
  * Born and Fock [1928] M. Born and V. Fock, Beweis des adiabatensatzes, Zeitschrift für Physik 51, 165 (1928).
  * Schiff [1949] L. I. Schiff,  _Quantum Mechanics_ (Mcgraw-Hill Book Co., Inc., New York, 1949).
  * Tong [2010] D. M. Tong, Quantitative condition is necessary in guaranteeing the validity of the adiabatic approximation, Phys. Rev. Lett. 104, 120401 (2010).
  * Roland and Cerf [2002] J. Roland and N. J. Cerf, Quantum search by local adiabatic evolution, Phys. Rev. A 65, 042308 (2002).
  * Kumar  _et al._ [2016] K. S. Kumar, A. Vepsäläinen, S. Danilin, and G. S. Paraoanu, Stimulated raman adiabatic passage in a three-level superconducting circuit, Nat. Commun. 7, 10628 (2016).
  * Vepsäläinen  _et al._ [2019] A. Vepsäläinen, S. Danilin, and G. S. Paraoanu, Superadiabatic population transfer in a three-level superconducting circuit, Sci. Adv. 5, eaau5999 (2019).
  * Dogra  _et al._ [2022] S. Dogra, A. Vepsäläinen, and G. S. Paraoanu, Experimental demonstration of robustness under scaling errors for superadiabatic population transfer in a superconducting circuit, arXiv:2203.12073 (2022).
  * Bergmann  _et al._ [1998] K. Bergmann, H. Theuer, and B. W. Shore, Coherent population transfer among quantum states of atoms and molecules, Rev. Mod. Phys. 70, 1003 (1998).
  * Falci  _et al._ [2017] G. Falci, P. G. Di Stefano, A. Ridolfo, A. D’Arrigo, G. S. Paraoanu, and E. Paladino, Advances in quantum control of three-level superconducting circuit architectures, Fortschritte der Physik 65, 1600077 (2017), https://onlinelibrary.wiley.com/doi/pdf/10.1002/prop.201600077 .


Defect-influenced particle advection in highly confined liquid crystal flows†

Magdalena Lesniewskaa, Nigel Mottramb, and Oliver Henricha‡

We study the morphology of the Saturn ring defect and director structure
around a colloidal particle with normal anchoring conditions and within the
flow of the nematic host phase through a rectangular duct of comparable size
to the particle. The changes in the defect structures and director profile
influence the advection behaviour of the particle, which we compare to that in
a simple Newtonian host phase. These effects lead to a non-monotonous
dependence of the differential velocity of particle and fluid, also known as
retardation ratio, on the Ericksen number.

##

††footnotetext: a Department of Physics, University of Strathclyde, Glasgow G4
0NG, UK.††footnotetext: b School of Mathematics & Statistics, University of
Glasgow, Glasgow G12 8QQ, UK. ††footnotetext: ‡ Email:
<EMAIL_ADDRESS>† Electronic Supplementary
Information (ESI) available: [details of any supplementary information
available should be included here]. See DOI: 10.1039/cXsm00000x/

##  1 Introduction

Microfluidics is concerned with the manipulation and control of fluid flow at
the microscale and sits as a multidisciplinary field at the intersection of
physics, engineering, and biology. During the last two decades it has seen a
tremendous rise in importance as it entered the main stream with many
practical applications ranging from medical diagnostics 1, 2 and drug delivery
3, 4 to chemical synthesis 5 and lab-on-a-chip technologies 6, 7, 8.

The essential physics of a microfluidic system is dictated by a competition
between various phenomena, which is captured by a series of dimensionless
numbers expressing their relative importance 9. The Reynolds number, for
instance, is often considered to be small in microfluidic applications.
However, relatively recently a focal point has been on inertial microfluidics
10, 11, 12, which gives rise to some interesting and counter-intuitive
phenomena.

When particle-laden flows are considered, aspects of confinement become
central to microfluidics. Characteristic dimensions of channels and chambers
are in the range of tens to hundreds of micrometres, and can be comparable to
the size of the particles being transported, and interactions between the
fluid and solid boundaries become increasingly important. In the simplest case
the fluid is Newtonian, without internal structure, and the geometry is a
uniform duct. The first theoretical results for migrating, rigid spheres in
unidirectional, two-dimensional flow were provided by Ho and Leal 13. These
were later extended to three dimensions and refined by Ganatos 14 and Staben
15, who also verified their theoretical results with an experimental study 16.
Owing to the rise of computer power and the advent of sophisticated simulation
methodologies, increasingly complex geometries can now be investigated 17.

For similar reasons, the study of particle migration in non-Newtonian fluids
with internal order structure has evolved only slowly. Suspensions of
colloidal particles in a nematic liquid crystalline host phase may serve here
as prototype of systems that cannot be described with standard continuum
theories, but require additional order parameters to capture the
microstructure and its change under flow conditions. After first theoretical
studies on the drag of colloidal particles in nematic hosts by Stark 18, the
main research focus has concerned the study of quiescent colloid-liquid
crystal suspensions 19, liquid crystalline droplets and shells 20, or flowing
pure liquid crystals 21, 22, 23. On the theoretical side, studies have been
extended to multiple particles and the full nemato-hydrodynamic problem that
solves for the tensor order parameter and velocity field 24. These approaches
are now complemented by various simulation methods 25, 26, 27. However,
practically nothing is known about particle migration in extreme confinement,
which forms the focus of this work.

Our paper is organised as follows: Section 2 introduces our theoretical
modelling framework, the Landau-de Gennes free energy and Beris-Edwards model,
while Section 3 gives details of our lattice Boltzmann simulation method.
Section 4 shows simulation results for scalar order parameter and director
field at various confinement ratios and flow velocities and them to those
obtained for simple Newtonian fluids. Section 5 summarises our results and
conclusions.

##  2 Theory

###  2.1 Landau-de Gennes free energy

The local order of the liquid crystal is described by a traceless and
symmetric second-order tensor $\bm{Q}(\bm{r},t)$ 28, 29. Its largest
eigenvalue $q<2/3$ is referred to as the scalar order parameter and provides a
measure of the liquid crystalline order at a certain position and time. The
eigenvector, $\bm{d}$, associated with $q$ is called the director and
describes the corresponding average orientation of the liquid crystal
molecules.

In equilibrium, the liquid crystal order is determined through minimisation of
its free energy, commonly described by the Landau-de Gennes free energy
functional

${\cal F}[\bm{Q}]=\int_{V}f(\bm{Q})\,dV+\int_{S}f_{s}(\bm{Q})\,dS,$ (1)

which includes the volume contribution $f=f_{b}+f_{g}$, that itself consists
of a bulk contribution $f_{b}$ and a gradient contribution $f_{g}$, and a
surface contribution $f_{s}$. The bulk free energy density is given by

$\displaystyle
f_{b}(\bm{Q})={\frac{A_{0}}{2}}\left(1-\dfrac{\gamma}{3}\right)Q_{\alpha\beta}^{2}-{\frac{A_{0}}{3}}\gamma\,Q_{\alpha\beta}Q_{\beta\pi}Q_{\pi\alpha}+{\frac{A_{0}}{4}}\gamma\,(Q_{\alpha\beta}^{2})^{2},$
(2)

where we use the Einstein summation convention, so that Greek indices that
appear twice are summed over. In Eq. 2, $A_{0}$ is a constant that sets the
overall energy scale and the parameter $\gamma$ controls the temperature
difference from the isotropic-nematic transition, and is related to a reduced
temperature $\tau$ by

$\tau=\frac{27}{\gamma}\left(1-\frac{\gamma}{3}\right).$ (3)

For $\gamma>3$ or $\tau<0$ the ordered, nematic state is the equilibrium
phase, whereas for $2.7\leq\gamma\leq 3$ or $0\leq\tau\leq 1$ the nematic
state is metastable. For $\gamma<2.7$ or $\tau>1$ the isotropic state is the
equilibrium phase.

The gradient free energy density $f_{g}$ contains the contributions of splay,
bend and twist deformations of the director field as well as order-elastic
effects due to gradients of the scalar order parameter,

$\displaystyle
f_{g}(\bm{Q})={\frac{1}{2}}\kappa_{0}(\partial_{\alpha}Q_{\alpha\beta})^{2}+{\frac{1}{2}}\kappa_{1}(\epsilon_{\alpha\sigma\nu}\partial_{\sigma}Q_{\nu\beta})^{2},$
(4)

where $\partial_{\alpha}=\partial/\partial x_{\alpha}$ and
$\epsilon_{\alpha\sigma\nu}$ is the Levi-Civita symbol in three dimensions. In
principle, the elastic constants $\kappa_{0}$ for splay and bend deformations,
and $\kappa_{1}$ for twist deformations can be different. However, in our
simulations we use the one elastic constant approximation
$\kappa_{0}=\kappa_{1}$.

The director is assumed to have a preferred normal orientation to the wall
surfaces and to the surface of the colloidal particle, known as a homeotropic
anchoring, and is described using a surface free energy term

$f_{s}(\bm{Q})={\textstyle\frac{1}{2}}w(Q_{\alpha\beta}-Q_{\alpha\beta}^{0})^{2},$
(5)

where $w$ is the surface anchoring strength with values $w_{wall}$ and
$w_{part}$ at the wall and particle surfaces, respectively. The preferred
orientation $Q^{0}_{\alpha\beta}$ is assumed to be uniaxial and is given by

$Q^{0}_{\alpha\beta}={\textstyle\frac{1}{2}}S_{0}(3\,n_{\alpha}n_{\beta}-\delta_{\alpha\beta}),$
(6)

where $\mathbf{n}$ is the surface unit normal, $\delta_{\alpha\beta}$ is the
Kronecker delta and $S_{0}$ is the preferred surface scalar order parameter
given by

$S_{0}=\frac{2}{3}\left(\frac{1}{4}+\frac{3}{4}\sqrt{1-\frac{8}{3\gamma}}\right).$
(7)

The anchoring strength at the surface of a colloidal particle is often
compared to the bulk fluid elastic constant by means of the dimensionless
parameter

$\omega=\frac{w\,R}{\kappa},$ (8)

where $R$ and $\kappa$ are the radius of the particle and the elastic
constant, respectively. For small values of this parameter, the presence of a
particle surface should have little impact on the local bulk liquid
crystalline ordering.

###  2.2 Beris-Edwards model

The time evolution of $Q_{\alpha\beta}$ is governed by the Beris-Edwards
equation 30

$\displaystyle\partial_{t}Q_{\alpha\beta}+\partial_{\gamma}(u_{\gamma}Q_{\alpha\beta})+S_{\alpha\beta}(\bm{W},\bm{Q})=\Gamma\,H_{\alpha\beta},$
(9)

where $\partial_{t}=\partial/\partial t$, $\bm{u}$ is the flow velocity,
$\bm{S}(\bm{W},\bm{Q})$ denotes the response to shear, $\bm{W}$ is the
velocity gradient tensor, $\bm{H}$ is the molecular field and $\Gamma$ is a
mobility parameter. The shear term is given by

$\displaystyle S_{\alpha\beta}$ $\displaystyle(\bm{W},\,\bm{Q})=(\xi
D_{\alpha\pi}+\Omega_{\alpha\pi})(Q_{\pi\beta}+{\textstyle\frac{1}{3}}\delta_{\pi\beta})$
$\displaystyle+(Q_{\alpha\pi}+{\textstyle\frac{1}{3}}\delta_{\alpha\pi})(\xi
D_{\pi\beta}-\Omega_{\pi\beta})-2\xi(Q_{\alpha\beta}+{\textstyle\frac{1}{3}}\delta_{\alpha\beta})Q_{\pi\sigma}W_{\sigma\pi}$
(10)

where $D_{\alpha\beta}=\frac{1}{2}(W_{\alpha\beta}+W_{\beta\alpha})$ and
$\Omega_{\alpha\beta}=\frac{1}{2}(W_{\alpha\beta}-W_{\beta\alpha})$ are the
symmetric and antisymmetric contributions to the velocity gradient tensor,
respectively, and $\xi$ is the so-called flow alignment parameter, a material
constant representing an effective molecular aspect ratio which determines
whether the liquid crystal exhibits a flow-aligned state at the Leslie angle
or tumbling state. The molecular field $\bm{H}$ is the functional derivative
of the free energy functional with respect to the order parameter,

$H_{\alpha\beta}=-\frac{\delta\cal F}{\delta
Q_{\alpha\beta}}+\frac{\delta_{\alpha\beta}}{3}{\rm Tr}\frac{\delta\cal
F}{\delta Q_{\alpha\beta}}.$ (11)

The second term in Eq. 11 involving the trace ensures tracelessness of the
tensor order parameter as it evolves through Eq.9. This leads to the following
molecular field:

$\displaystyle H_{\alpha\beta}=$ $\displaystyle-
A_{0}(1-\gamma/3)Q_{\alpha\beta}+A_{0}\gamma\,(Q_{\alpha\sigma}Q_{\sigma\beta}-{\textstyle\frac{1}{3}}Q_{\sigma\nu}^{2}\delta_{\alpha\beta})$
$\displaystyle-
A_{0}\gamma\,Q_{\sigma\nu}^{2}Q_{\alpha\beta}+\kappa_{0}\partial_{\alpha}\partial_{\sigma}Q_{\sigma\beta}+\kappa_{1}\partial_{\sigma}(\partial_{\sigma}Q_{\alpha\beta}-\partial_{\alpha}Q_{\sigma\beta})$

The governing equations of hydrodynamic motion are the equation of mass
conservation, also known as the continuity equation, and the Navier-Stokes
equation that describes the conservation of linear momentum. In tensor
notation they read

$\partial_{t}\rho+\partial_{\alpha}(\rho u_{\alpha})=0$ (13)

and

$\displaystyle\partial_{t}(\rho$ $\displaystyle
u_{\alpha})=\partial_{\beta}\Pi^{(LC)}_{\alpha\beta}+\partial_{\beta}\Pi^{(HD)}_{\alpha\beta},$
(14)

respectively. Eq. 13 relates the local rate of change of the density $\rho$ to
the advection of mass by the fluid velocity $\bm{u}$. Eq. 14 is Newton’s
second law of momentum change for the fluid and involves the thermotropic
stress tensor $\Pi^{(LC)}_{\alpha\beta}$ and the hydrodynamic stress tensor
$\Pi_{\alpha\beta}^{(HD)}$. The thermotropic stress arises due to the liquid
crystal and is given by

$\Pi^{(LC)}_{\alpha\beta}=\sigma_{\alpha\beta}+\tau_{\alpha\beta}-\partial_{\alpha}Q_{\sigma\nu}\frac{\delta{\cal
F}}{\delta\partial_{\beta}Q_{\sigma\nu}}.$ (15)

In Eq. 15, $\sigma_{\alpha\beta}$ and $\tau_{\alpha\beta}$ are the symmetric
and antisymmetric stress contributions, respectively, defined as

$\displaystyle\sigma_{\alpha\beta}$ $\displaystyle=$ $\displaystyle-
p_{0}\,\delta_{\alpha\beta}-\xi
H_{\alpha\sigma}(Q_{\sigma\beta}+{\textstyle\frac{1}{3}}\delta_{\sigma\beta})-\xi(Q_{\alpha\sigma}+{\textstyle\frac{1}{3}}\delta_{\alpha\sigma})H_{\sigma\beta}$
(16)
$\displaystyle+2\xi(Q_{\alpha\beta}+{\textstyle\frac{1}{3}}\delta_{\alpha\beta})Q_{\sigma\nu}H_{\sigma\nu},$

where $p_{0}=-\left(\partial{\cal F}/\partial V\right)_{T}=-f$ is the
isotropic contribution from the nematic liquid crystal to the total pressure,
and

$\tau_{\alpha\beta}=Q_{\alpha\sigma}H_{\sigma\beta}-H_{\alpha\sigma}Q_{\sigma\beta}.$
(17)

The final term in Eq. 15 may be expanded as

$\begin{split}\partial_{\alpha}Q_{\sigma\nu}\frac{\delta{\cal
F}}{\delta\partial_{\beta}Q_{\sigma\nu}}=&-\kappa_{0}\partial_{\alpha}Q_{\sigma\beta}\partial_{\nu}Q_{\sigma\nu}\\\
&-\kappa_{1}\partial_{\alpha}Q_{\sigma\nu}\left(\partial_{\beta}Q_{\sigma\nu}-\partial_{\sigma}Q_{\nu\beta}\right).\end{split}$
(18)

The hydrodynamic stress tensor is defined as

$\Pi_{\alpha\beta}^{(HD)}=-p\,\delta_{\alpha\beta}-\rho
u_{\alpha}u_{\beta}+\mu(\partial_{\beta}u_{\alpha}+\partial_{\alpha}u_{\beta})+\zeta\partial_{\sigma}u_{\sigma}\delta_{\alpha\beta},$
(19)

where $\mu$ and $\zeta$ are the dynamic and bulk viscosity, respectively. The
hydrostatic pressure $p$ is related to the density via an ideal gas equation
of state as $p=c_{s}^{2}\rho$ with $c_{s}$ as lattice speed of sound as is
standard in the lattice Boltzmann method. The last term vanishes in
incompressible fluids as Eq. 13 becomes $\partial_{\alpha}u_{\alpha}=0$.

No-slip and no-penetration boundary conditions are applied on the walls and
particle surfaces, and the boundary conditions for $\bm{Q}$ are found from the
minimisation of the free energy 31

$n_{\gamma}\frac{\partial f}{\partial Q_{\alpha\beta,\gamma}}+\frac{\partial
f_{s}}{\partial Q_{\alpha\beta}}=0,$ (20)

where $Q_{\alpha\beta,\gamma}=\partial Q_{\alpha\beta}/\partial x_{\gamma}$.

##  3 Simulation Method

###  3.1 Simulation Setup

Fig 1 shows a diagram of the three-dimensional computational geometry, which
consists of a duct of $L_{x}=24$, $32$ or $48$ and $L_{y}\times
L_{z}=256\times 384$ lattice sites. Solid walls are positioned at $x=0$ and
$x=L_{x}$, $y=0$ and $y=L_{y}$. We define the measure of confinement as the
ratio of the particle diameter to the height of the duct, which leads in our
case to confinement ratios $2R/L_{x}=0.8$, $0.6$, and $0.4$. The value of
$L_{y}$ means that, with the particle at the centre of the duct, the system is
effectively unconfined in $y$-direction since $2R/L_{y}=0.075$. Periodic
boundary conditions are applied in the $z$-direction with the $z$-boundaries
acting as inlet and outlet of the duct. A pressure gradient $\Psi=\Delta
p/L_{z}$ is applied in $z$-direction, leading to a body force density acting
on all sites.

$L_{y}$ $L_{x}$ radius $R$ No slip, no-penetration,homeotropic anchoringflow

Fig. 1:  Overview of the computational geometry: We apply no-slip and no-
penetration boundary conditions and homeotropic anchoring conditions at the
walls which bound the region in $x$-direction and $y$-direction and at the
particle surface, with periodic boundary conditions at the $z$-boundaries. The
top part shows the top view, looking along the $x$-direction, and the bottom
part shows the side view, looking along the $y$-direction.

We use a hybrid lattice Boltzmann scheme 32 that applies a finite-difference
method for the dynamics of the $\bm{Q}$-tensor order parameter and solves the
hydrodynamic part of the problem by means of the lattice Boltzmann method. The
colloidal particle is discretised as a solid, mobile particle with a radius
$R=9.6$. On both the walls and the particle surface no-slip and no-penetration
boundary conditions are applied by using a bounce-back on links scheme 33, 34,
35. The surface free energy in Eq. 5 invokes a homeotropic anchoring condition
with a preferred orientation of the director normal to the surfaces.

There are technical limitations to our model that should be borne in mind.
While the centre of mass of the particle is integrated off-grid according to
Newton’s equation, the particle itself is discretised using a stair-case
geometry. This can result in some inaccuracies, especially for highly confined
regimes. For instance in our case $2R/L_{x}=0.8$ only about two lattice sites
are between the particle and the walls surfaces at its narrowest point. The
pressure obeys an ideal gas equation of state and is directly related to the
density via $p=c_{s}^{2}\rho$, while the effect of the constant pressure
gradient is modelled through an additional body force density on all sites.
Both treatments are common in the lattice Boltzmann methodology and allow for
an accurate modelling of a weakly compressible fluid, but the assumption of a
constant pressure gradient represents obviously a simplification over the real
situation. Thermal fluctuations have not been included since our simulations
were carried out using a temperature well away from the isotropic-nematic
transition line, and so elastic forces from the anchoring of the liquid
crystal dominate over thermal forces by orders of magnitude. Our approach uses
the one-constant approximation to simplify elasticity where the elastic
constants for splay, bend and twist deformations carry the same value. This
approximation is commonly used as first approach and does not compromise our
results qualitatively. However, relaxing this approximation could lead to
quantitative differences, and potentially also richer phenomenology. The
Beris-Edwards model uses a simplified approach to viscosities compared to the
Ericksen-Leslie theory, which has six viscosity coefficients
$\alpha_{1},\dots,\alpha_{6}$ 28 (only five are independent as the Parodi
relation applies 36). The viscosities in the Beris-Edwards model are
implicitly given through the isotropic dynamic shear viscosity $\mu$, the
rotational diffusion constant $\Gamma$, the flow alignment parameter $\xi$,
and the scalar order parameter $S_{0}$. They can be directly related to the
Ericksen-Leslie viscosities $\alpha_{1},\dots,\alpha_{6}$ 32, but parameterise
only a subset of possible values.

The simulations were run with our lattice Boltzmann code for complex fluids
Ludwig version 0.15.0 37. A typical simulation is first initialised with no
applied pressure gradient for $5\times 10^{4}$ iteration steps for each
anchoring strength. After this initial equilibration phase, the simulations
are restarted with various pressure gradients that are kept constant for
$4\times 10^{5}$ iteration steps. Typical runtimes are approximately 26 hours
using a hybrid MPI/OpenMP parallelisation with 2 MPI tasks each running on 20
OpenMP threads. The overview of the simulation parameters is included in Table
1. For further information about the exact implementation used in this work,
we guide the reader to the Ludwig code repository and related literature 38,
39.

Bulk energy scale | $A_{0}$ | $0.01$
---|---|---
Effective temperature | $\tau$ | $-0.29$
Elastic constants | $\kappa_{0},\kappa_{1}$ | $0.01$
Wall anchoring strength | $w_{wall}$ | $0.02$
Particle anchoring strength | $w_{part}$ | 0, 0.001, 0.01, 0.05
Anchoring parameter | $\omega$ | 0, 0.96, 9.6, 48
Flow alignment parameter | $\xi$ | $0.7$
Mobility parameter | $\Gamma$ | $0.5$
Density | $\rho$ | 1.0
Shear viscosity | $\mu$ | $5/6$
Bulk viscosity | $\zeta$ | $5/6$
Particle radius | $R$ | $9.6$
Table 1: Overview of simulation parameters

###  3.2 Parameter Mapping

Our simulation units can be mapped to physical units by calibrating the units
of pressure, time, and length. To achieve this, we assign the lattice spacing
$\Delta x$, the algorithmic time step $\Delta t$ and the reference pressure
$p^{*}$ from unity in lattice Boltzmann units (LBU) to their corresponding
values in SI units. The principle of this parameter mapping was also shown in
our previous work 27 using a different characteristic length scale.

The calibration of the length scale is straightforward as it is simply set by
considering the dimensions of the microfluidic duct. If we associate the
narrowest gap size $L_{x}=24$ in LBU corresponds to
$L_{x}\,\widehat{=}\,1.2\times 10^{-6}$ m in SI units, we obtain an LBU length
of $\Delta x\,\widehat{=}\,5\times 10^{-8}$ m $=50$ nm in SI units.

To obtain the pressure scale, we use the measurements of the Landau-de Gennes
parameters 29 (Appendix D therein) which give

$\displaystyle\frac{27}{2\,A_{0}\,\gamma}$ $\displaystyle\simeq$
$\displaystyle 5\times 10^{-6}\mbox{ J}^{-1}\mbox{m}^{3}=5\times 10^{-6}\mbox{
Pa}^{-1}.$

Using $A_{0}=0.01$ and $\gamma=3.1$ in our simulations results in a reference
pressure of $p^{*}=1$ LBU $\widehat{=}\,10^{8}$ Pa in SI units.

For the calibration of the timescale we use the following formula, which
relates the rotational viscosity $\gamma_{1}$ of the director to the
equilibrium scalar order parameter $q$ and the order parameter mobility
$\Gamma$:

$\gamma_{1}=\frac{2q^{2}}{\Gamma}$

We use $\Gamma=0.5$ in LBU and bulk energy density parameters that lead to
$q\approx 0.5$ since it is assumed that the system is well within the nematic
phase. Therefore, the rotational viscosity $\gamma_{1}=1$ in LBU. Typical
values for liquid crystals in SI units are $\gamma_{1}=0.1$ Pa s 28. Together
with $1$ Pa equating to a pressure of $10^{-8}$ in LBU, we obtain for the
algorithmic time step $\Delta t\,\widehat{=}\,10^{-9}$ s $=1$ ns.

The Ericksen number characterises the ratio of viscous to elastic forces and
is defined as

$Er=\frac{\eta\,u\,\Lambda}{\kappa},$

where $u$ is a characteristic flow velocity, in our case the velocity at the
centre of the duct $U_{c}$, $\eta$ is the dynamic viscosity, $\Lambda$ is a
characteristic length scale which is set by the narrowest gap size $L_{x}$
(see Table 2 for $\Lambda=2R$, which allows direct comparison with our
previous work 27), and $\kappa$ is the bulk elastic constant of the liquid
crystal.

Fig. 2: Normalised velocity profiles for $2R/L_{x}=0.6$, where $u$ is the
fluid velocity, and $u_{c}$ is the fluid velocity at the centre line. The
black line is the standard Poiseuille flow profile. The inset shows the
dependence of the centre line velocity $u_{c}$ on the Ericksen number $Er$.

The dynamic viscosity $\eta$ is calculated as an apparent viscosity, defined
as the ratio, $\eta=\mu\Phi_{0}/\Phi$ of the volumetric flux $\Phi_{0}$ of a
simple Newtonian fluid and the volumetric flux of the liquid crystalline
system $\Phi$, through a plane perpendicular to the flow in the $z$-direction,
namely

$\Phi=\int_{0}^{L_{x}}\int_{0}^{L_{y}}v_{z}(x)\;dx\,dy,\\\ $ (21)

with the flow being driven through the pressure gradient $\Psi=\Delta p/L_{z}$
with $\Delta p$ being the pressure difference between inlet and outlet. The
volumetric flow rate $\Phi_{0}$ of a Newtonian fluid with dynamic viscosity
$\mu$ through a gap $L_{x}$ driven by a pressure gradient $\Psi$ in plane
Poiseuille flow can be calculated as

$\displaystyle\Phi_{0}$ $\displaystyle=$
$\displaystyle\int_{0}^{L_{x}}\int_{0}^{L_{y}}\frac{L_{x}^{2}}{2\mu}\Psi\left(\frac{x}{L_{x}}-\left(\frac{x}{L_{x}}\right)^{2}\right)dx\,dy$
(22) $\displaystyle=$ $\displaystyle\frac{L_{x}^{3}\,L_{y}\,\Psi}{12\,\mu}.$
(23)

Fig. 2 shows velocity profiles for a confinement ratio $2R/L_{x}=0.6$,
normalised to the Poiseuille velocity profile that a simple Newtonian fluid
would have at the same pressure gradient. The apparent viscosity $\eta$ is the
ratio of the areas under the Poiseuille curve and curves at finite Ericksen
numbers.

$2R/L_{x}=0.4$ | $2R/L_{x}=0.6$ | $2R/L_{x}=0.8$
---|---|---
Er (Lx) | Er (2R) | Er (Lx) | Er (2R) | Er (Lx) | Er (2R)
1.65 | 0.658 | 4.38 | 2.63 | 6.15 | 4.92
8.30 | 3.32 | 9.84 | 5.90 | 10.37 | 8.30
18.10 | 7.24 | 21.25 | 12.75 | 17.95 | 14.36
35.55 | 14.22 | 39.18 | 23.51 | 22.32 | 17.86
52.72 | 21.09 | 51.86 | 31.11 | 64.85 | 52.16
69.77 | 27.91 | 102.12 | 61.27 | 86.06 | 69.10
Table 2: Conversion of Ericksen numbers for different confinement ratios using
different characteristic length scales, namely the size of the channel (odd
columns) or as in Ref. 27 the diameter of the particle (even columns).

##  4 Results and Discussion

We study the migration behaviour of a single particle moving in a nematic host
phase in highly confining ducts and investigate the effect that varying
pressure gradient $\Psi$, confinement ratio $2R/L_{x}$ and homeotropic
anchoring strength have. In a simple Newtonian fluid, or in a liquid crystal
at temperatures above the isotropic-nematic transition point, the motion of a
freely suspended spherical particle between two parallel plane walls has been
studied previously theoretically 14, 15, with simulations 17 and
experimentally 16. The main effect is that the relative retardation of the
particle motion to the fluid motion is primarily independent of the applied
pressure gradient, but greater for particles closer to either of the walls,
and for highly confined particles due to the proximity to the walls.

In a nematic liquid crystal with homeotropic anchoring conditions at the walls
the director orientation is forced to be parallel to the wall normals. The
degree of alignment depends on the strength of the anchoring, but also on the
velocity gradient, and therefore the pressure gradient. At low pressure
gradients, the nematic order will be enforced throughout the duct. But for
higher pressure gradients the director field flow-aligns at the Leslie angle.
Two conformations are persistent in flowing nematics, namely the so-called
bend state or H-state and the splay state or V-state. For both H- and V-state
the director flow-aligns to a positive (negative) Leslie angle in the lower
(upper) half of the channel. The difference between the two states is
determined by the way the director rotates between the positive and negative
Leslie angles at the centre: In the bend state, the director at the centre is
perpendicular to the walls, whereas in the splay state the director is almost
parallel to the walls at the centre. The bend state is generally adopted at
low flow velocities, whereas the nematic transitions to the splay state at
higher flow velocities.

Er=4.38  | Er=4.38
---|---
|
Er=51.86  | Er=51.86
|
|
Fig. 3: Director field for $2R/L_{x}=0.6$ confinement ratio before and after
the bend-to-splay transition. The top row shows the bend state (H-state) at
Ericksen number $Er=4.38$, while the bottom row shows the splay state
(V-state) at Ericksen number $Er=51.86$. The images in the left column
represent a slice through the middle of the channel in $xz$-plane (narrowest
duct dimension and flow direction) and have the view along the negative
$y$-dimension. Those in the right column show a slice in $yz$-plane (widest
and narrowest duct dimension) and have the view in positive $x$-direction. The
flow direction is from left to right in positive $z$-direction. The opacity of
the Saturn ring defect (green isosurface) has been reduced to show the
director structure. The colour code indicates the magnitude of the
$z-$-component of the director field in flow direction, where blue means
alignment in $xy$-plane and red means full alignment in $z$\- or flow
direction. The opacity of the defect rings (green isosurfaces) is slightly
reduced to enhance visibility of the local director field.

Fig. 3 shows the director field and defect structure at a medium confinement
ratio $2R/L_{x}=0.6$ and different Ericksen numbers. The top row displays the
bend state prior to the transition to the splay state, whereas the bottom row
shows the splay state after transitioning from the bend state. The left column
contains slices in the $xz$-plane (narrowest duct dimension and flow
direction) at $y=L_{y}/2$ with walls at the $x$-boundaries at the top and
bottom, whereas the right column shows slices in the $yz$-plane (widest duct
dimension and flow direction) at $x=L_{x}/2$ cropped to the vicinity of the
colloidal particle.

$2R/L_{x}=0.4$ |  $2R/L_{x}=0.6$ |  $2R/L_{x}=0.8$
---|---|---
$Er=1.65$ |  | $Er=4.38$ |  | $Er=6.15$ |
|  |  |  |  |
$Er=8.30$ |  | $Er=9.84$ |  | $Er=10.37$ |
|  |  |  |  |
$Er=18.10$ |  | $Er=21.25$ |  | $Er=22.32$ |
|  |  |  |  |
$Er=52.72$ |  | $Er=51.86$ |  | $Er=64.85$ |
|  |  |  |  |
$Er=69.77$ |  | $Er=102.12$ |  | $Er=86.06$ |
|  |  |  |  |
Fig. 4: Snapshots of the director field and defect structure in the steady
state at various Ericksen numbers and strongest particle anchoring parameter
$\omega=48$. The bright green region corresponds to the defect where liquid
crystalline order is reduced. In each cell the left images in each cell show
the side view looking in the negative $y$-direction with walls at the top and
bottom. The images on the right in each cell show the view from the top
looking in the positive $x$-direction. The confinement increases from left to
right, and Ericksen numbers increase from top to bottom. While in most cases
the particle stays at the centre of the duct throughout the simulation, there
are a small number of cases where they migrate away from it. Specifically,
there are two cases where the particle migrates fully to a wall (for the two
highest Ericksen number and confinement ratio $2R/L_{x}=0.8$), and three cases
where the particle migrates to a stable position between the wall and the
centre (for $Er=18.10$ and $2R/L_{x}=0.4$, for $Er=9.84$ and $2R/L_{x}=0.8$,
and for $Er=10.37$ and $2R/L_{x}=0.8$).

The director field is colour-coded with red indicating an orientation parallel
to the flow direction and blue indicating an orientation perpendicular to the
flow direction. The results presented in Fig. 3 give a representative example
of the director structures around a particle before and after the bend-to-
splay transition. The bend state shows the Saturn ring defect oriented
parallel to the walls with only very minor deformations, while the splay state
has the Saturn ring defect oriented approximately perpendicular to the walls.

As a quantitative overview of our findings, we include in Fig. 4, snapshots of
the particle and its defect in the steady state for varying confinement ratios
and Ericksen numbers. In each cell the left images show the side view looking
in the negative $y$-direction with walls at the top and bottom. The images on
the right show the view from the top looking in the positive $x$-direction.
The confinement increases from left to right from confinement ratios
$2R/L_{x}=0.4$ to $2R/L_{x}=0.6$ to $2R/L_{x}=0.8$, and Ericksen numbers
increase from top to bottom. The defect is shown as a green isosurface defined
by a local order parameter $q\leq 0.188$ and the particle anchoring strength
and dimensionless anchoring parameter are $w_{part}=0.05$ and $\omega=48$,
respectively, as lower anchoring strengths do not result in defects that could
be distinctively visualised.

At low Ericksen numbers, below the bend-to-splay transition, the particle has
a Saturn ring defect whose ring plane remains parallel to the walls. This is
the case for all confinement ratios and Ericksen numbers below $Er=10.37$, as
shown in the first and second row of Fig. 4. Two aspects are noteworthy:
Firstly, there is a slight increase of the defect isosurface radius downstream
of the particle, for which both Ericksen number (see images for $2R/L_{x}=0.4$
with $Er=1.65$ and $Er=8.30$) and confinement ratio (see image for
$2R/L_{x}=0.6,Er=4.38$ and image and Movie.1 for $2R/L_{x}=0.8,Er=6.15$) seem
responsible. However, confinement appears to play a more important role in
this context. Secondly, at slightly increased Ericksen numbers (see images for
$2R/L_{x}=0.6,Er=9.84$ and $2R/L_{x}=0.8,Er=10.37$), the Saturn ring becomes
angled such that the part downstream of the particle is closer to the bottom
wall, while the other part upstream of the particle remains virtually
unchanged. These two particular cases reached steady state positions that are
offset somewhere between the centre of the duct and the walls in the
$x$-direction, which contributes to this asymmetric appearance. This can be
again explained with the migration to the weak attractor region that we
observed in our previous work on controllable particle migration 27 in
practically unconfined conditions using a much wider duct and lower
confinement ratio $2R/L_{x}=0.15$. For direct comparison we provide in Table 2
an approximate conversion between particle Ericksen numbers, as used our
previous publication, and Ericksen numbers based on the smallest duct
dimension, as used in this work.

Fig. 5 shows a direct comparison of the defect rings around the particle for
the lowest and highest simulated Ericksen numbers below the bend-to-splay
transition, at confinement ratios (a) $2R/L_{x}=0.4$ (b) $0.6$ and (c) $0.8$.
The defects at the lowest Ericksen numbers, depicted in grey, are distinctive
Saturn rings that are oriented parallel to the walls in the $x-$direction. As
previously mentioned, increasing confinement leads to a defect ring that is
thicker at the downstream side of the particle, while it remains oriented
parallel to the wall at the $x$-boundary. Increasing the Ericksen number alone
does not change the orientation of the defect ring, but leads to a very slight
shift in position upstream (see grey and yellow defect rings in Fig. 5 a)).
However, increasing the Ericksen number and confinement ratio induces a
noticeable tilt of the defect ring, shown in Fig. 5 b) and c), as the particle
migrates into the off-centre steady state position somewhere between the
centre of the duct and one of the walls. Despite the difference in confinement
ratio and Ericksen number ($2R/L_{x}=0.6$ and $0.8$, $Er=9.84$ and $10.37$,
respectively) the shape of the defect rings is almost the same.

a) | b) | c)
---|---|---
|  |
|  |
Er=1.65 Er=8.30 | Er=4.38 Er=9.84 | Er=6.15 Er=10.37
Fig. 5: Disclination lines around the particle for the lowest  — and highest
— simulated Ericksen numbers below the bend-to-splay transition at confinement
ratios a) $2R/L_{x}=0.4$, b) $2R/L_{x}=0.6$ and c) $2R/L_{x}=0.8$,
respectively. The top row has the view in negative $y$-direction, the widest
duct dimension, with the walls in the narrowest duct dimension at the
$x$-boundaries situated closely above and below the particle. The bottom row
shows the view in positive $x$-direction. The flow is in the horizontal
positive $z$-direction from left to right.

Upon increasing the Ericksen number, a bend-to-splay transition takes place
somewhere between $8.30<Er<18.10$ (for $2R/L_{x}=0.4$), $9.84<Er<21.25$ (for
$2R/L_{x}=0.6$) and $10.37<Er<17.95$ (for $2R/L_{x}=0.8$). The defect ring is
now reoriented with its ring plane approximately perpendicular to the walls
and flow direction, as shown in Fig. 4, for instance in the third row, and
retains a similar shape at higher Ericksen numbers (see images for
$2R/L_{x}=0.4,Er=18.10$, $2R/L_{x}=0.6,Er=21.25$, and
$2R/L_{x}=0.8,Er=22.32$). The case for $2R/L_{x}=0.4,Er=18.10$ forms an
exception in that the particle moves very slightly away from the centre into a
stable off-centre position, while in the other cases the particle remains at
the centre of the duct, which can be also understood with the migration to the
previously observed weak attractor region at similar Ericksen numbers 27 (see
Table 2 for conversion of Ericksen numbers). A noticeable difference is that
with increasing confinement the defect ring appears compressed in the smallest
duct dimension due to the relative proximity of the walls (see image and
Movie.2 for $2R/L_{x}=0.8,Er=22.32$).

With increasing Ericksen numbers the shape of the vertically oriented defect
ring remains largely unchanged for low and medium confinement, as shown in the
first and second column, forth and fifth row, of Fig. 4 ($2R/L_{x}=0.4$ and
$0.6$) for Ericksen numbers $Er=52.72,\,69.77$ and $Er=51.86,\,102,12$,
respectively. At $Er=102.12$ a slight change occurs such that the defect close
to the midplane of the duct in the $x$-direction are distorted and pulled in
the upstream direction, i.e. against the flow. This effect is a precursor of
the more dramatic elongation of the Saturn ring that will become even more
evident as the confinement ratio increases.

At even higher Ericksen numbers $Er=64.85$ and $Er=86.06$ and confinement
$2R/L_{x}=0.8$, shown in the third column forth and fifth row of Fig. 4, we
observe defects that differ substantially from those discussed before. In
these cases the particle migrates fully to one of the walls. This has also
been previously observed for similar Ericksen numbers in much lower
confinement 27. The migration to one of the walls results in a different
defect shape such that there is a pronounced elongation of the Saturn ring
defect in the upstream direction. There is also the indication of a small
satellite region of low order, upstream of the particle that never merges up
with the rest of the defect (see image and Movie.3 for $2R/L_{x}=0.8$ and
$Er=86.06$).

a) | b) | c)
---|---|---
|  |
|  |
Er=18.10 Er=35.55 | Er=21.25 Er=39.18 | Er=17.95 Er=22.32
Er=52.72 Er=69.77 | Er=51.86 Er=102.12 | Er=64.85 Er=86.06
Fig. 6: Saturn ring disclination lines around the particle for various
Ericksen numbers after the bend-to-spay transition has taken place. The
confinement ratios are a) $2R/L_{x}=0.4$, b) $2R/L_{x}=0.6$ and c)
$2R/L_{x}=0.8$. In the top row the view is along the negative $y$-direction,
the widest duct dimension, with the walls in the narrowest duct dimension at
the $x$-boundaries situated above and below the particle. The bottom row is
the view in the positive $x$-direction. The flow is in the horizontal
$z$-direction from left to right.

Before focusing on the director structure at high Ericksen numbers and large
confinement in more detail (see Fig. 7), we present briefly a synopsis of the
defect rings at different confinement ratios and Ericksen numbers. Fig. 6
shows superimposed, vertically oriented defect rings as they occur after the
bend-to-splay transition has taken place. At the lowest confinement ratio
$2R/L_{x}=0.4$, shown in Fig. 6 a), the defect ring remains relatively
undistorted across a range of medium to high Ericksen numbers. However,
comparing the images at the top with the view along the $y$-direction across
the narrowest gap to those at the bottom with the view along the $x$-direction
across the widest gap gives evidence that the shape of the Saturn ring defects
is sensitive to confinement. When confined, the defect rings are located
slightly downstream from the particle’s equator, whereas they remain situated
along the equator in the dimension of no or very small confinement
($2R/L_{y}=0.075$). This feature becomes more pronounced as the confinement
increases, discernible through the green defect rings at ratios $2R/L_{x}=0.6$
in Fig. 6 b), and more so at $2R/L_{x}=0.8$ in Fig. 6 c) where it results in
the compressed appearance (Fig. LABEL:6 c) top row). This applies to lower
(green isosurfaces) and medium (orange isosurfaces) Ericksen numbers.
Increasing both confinement and Ericksen numbers leads to the aforementioned
different appearance of the defect rings (purple and magenta isosurfaces).

It is worth mentioning that the confinement ratios we studied are larger than
those in similar studies 25, 26 ($2R/L_{x}=0.25$ and $2R/L_{x}=0.19$,
respectively), where the lower confinement has been chosen to eliminate
possible effects on the results. However, our case of $2R/L_{x}=0.4$ is
obviously already low enough to feature defect rings that appear undeformed
and occur at unaltered relative positions to the particle.

$Er=22.32$ | $Er=64.85$ | $Er=86.06$
---|---|---
|  |
|  |
Fig. 7:  Director field for $2R/L_{x}=0.8$ confinement ratio after the bend-
to-splay transition for increasing Ericksen numbers $Er=22.32,64.85$ and
$86.06$, respectively. The images in the top row represent a slice through the
middle of the channel in $xz$-plane (narrowest duct dimension and flow
direction) and have the view along the negative $y$-dimension. Those in the
bottom row show a slice in $yz$-plane (widest and narrowest duct dimension)
and have the view in positive $x$-direction. The flow direction is from left
to right in positive $z$-direction. The opacity of the Saturn ring defect
(green isosurface) has been reduced to show the director structure. The colour
code indicates the magnitude of the $z$-component of the director field in
flow direction, where blue means alignment in $xy$-plane and red means full
alignment in $z$\- or flow direction. The opacity of the defect rings (green
isosurfaces) is slightly reduced to enhance visibility of the local director
field.

The director structure at high Ericksen numbers and the largest confinement
ratio $2R/L_{x}=0.8$ is shown in Fig. 7. At this confinement ratio the walls
at the $x$-boundaries are close to the colloidal particle. The left column
shows the situation at moderately high Ericksen numbers $Er=22.32$. The defect
ring is vertically oriented, noticeably displaced upstream in the confined
$x$-dimension (Fig. 7 top left image), and situated at the equatorial region
of the particle in the non-confined $y$-dimension (Fig. 7 bottom left image).
Consideration of the director field reveals that in the confined dimension the
homeotropic anchoring conditions at the wall and particle surfaces prevent any
kind of flow alignment in the narrow gap between the particle and the walls.
Flow alignment occurs only further away upstream and downstream of the
particle. But while flow alignment and anchoring are compatible in the wake of
the particle, they are not at its bow, leading to the displacement of the
defect ring further upstream. This is different in the non-confined dimension
where the flow alignment occurs in a short distance from the particle and
entails a defect in the equatorial region.

At higher Ericksen numbers $Er=64.85$ the particle migrates readily to one of
the walls 27 and the shape of the defect changes markedly. The asymmetric
positioning of the particle in the duct is only partly responsible for this.
In fact, we observe large differential velocity between the particle and the
fluid, which means the particle acts now increasingly as obstacle and the
fluid has to flow around it. This results in a region with redirected flow and
flow alignment immediately upstream of the particle as visible in the cut
through the director field in $xz$-plane (Fig. 7 top centre image), but also
on either side of the particle as the cut through the director field in
$yz$-plane shows (Fig. 7 bottom centre image). The latter causes the defect
ring to become extended further upstream, albeit it never completely engulfs
the particle. This feature becomes more pronounced as we increase the Ericksen
numbers further to $Er=86.06$ (Fig. 7 third column). Here, a satellite region
of very low order emerges where the fully flow-aligned director meets the
region of redirected flow in the vicinity and upstream of the particle.

Fig. 8: Comparison of retardation ratios $v/u_{c}$ of particle velocity $v$ to
fluid velocity $u_{c}$ at the centre of the rectangular duct for confinement
ratios $2R/L_{x}=0.4$ (blue), $0.6$ (green) and $0.8$ (red) and different
particle anchoring strengths. Horizontal lines show results in a Newtonian
fluid from Staben et al. 15 (dashed-dotted lines) and with our approach (solid
lines). Open symbols indicate cases where the colloidal particle has been
fixed in $x$-direction for comparison as it would normally migrate away from
the centre of the duct to either an off-centre position or to the walls. The
vertical lines indicate the approximate position of the bend-to-splay
transition.

We conclude our study with an analysis of the migration behaviour of the
colloidal particle at different Ericksen numbers and confinement ratios and
compare it to that in a simple Newtonian fluid. For this purpose we draw on
the theoretical results obtained by Staben et al. 15, which have been
reproduced in a number of studies.

A suitable measure to characterise the migration behaviour is the retardation
ratio $v/u_{c}$ of particle velocity $v$ to fluid velocity $u_{c}$ at the
centre of the duct. In an isotropic Newtonian fluid under Poiseuille flow this
ratio is constant and depends only on the distance of the particle from the
walls of the duct and the confinement ratio. In particular, $v/u_{c}$ is
independent of the Reynolds number. Without confinement the retardation ratio
$v/u$ is unity as the particle acts as a tracer and is simply advected with
the fluid. At finite confinement ratios below $2R/L_{x}=1$ the movement of the
particle is slowed down in the parabolic Poiseuille flow due to the no-slip
boundary conditions on the walls of the duct.

Fig. 8 shows the retardation ratio $v/u_{c}$ for different confinement ratios
$2R/L_{x}$, particle anchoring strengths and Ericksen numbers $Er$. Using the
Ericksen number as abscissa has the advantage that the bend-to-splay
transition occurs at similar values aiding the comparison across different
confinement ratios. The straight, horizontal lines represent the results for a
particle in a simple Newtonian fluid. Dashed-dotted lines give the results
from Staben et al. 16 for particles at the centre of the duct. We measure
retardation ratios of $v/u_{c}=0.946,0.876$ and $0.759$ for confinement ratios
$2R/L_{x}=0.4,0.6$ and $0.8$, respectively, shown in Fig. 8 with solid lines.
These results compare well with those of Staben et al., which are
$v/u_{c}=0.945,0.871$ and $0.746$ for the same confinement ratios and
particles positioned at the centre of the duct. It is worth emphasising that
in our setup the largest confinement ratio $2R/L_{x}=0.8$ has less than 3
lattice sites between the particle surface and the walls on either side.
Nevertheless, the relative deviation between ours and Staben’s results for
Newtonian host phases is less than $1.8\%$ in the worst case, which means our
method is remarkably accurate given the relatively sparse discretisation.
However, it should be borne in mind that when modelling a liquid crystalline
host phases the sparse discretisation affects also the tensor order parameter
$\bm{Q}$ in addition to the fluid-solid interaction in a Newtonian host phase.
While these limitations affect the results in Fig. 8 to a certain extent,
there are nevertheless clear and robust trends that we will now discuss.

At low Ericksen numbers we observe retardation ratios $v/u_{c}$ that are close
or identical to their corresponding values in Newtonian fluids. Interestingly,
and primarily for no or low particle anchoring strength and low confinement
ratios $2R/L_{x}=0.4$ and $0.6$, the retardation ratio can be slightly larger
in the nematic host phase than in the Newtonian host phase (see light and
medium grey and blue data points in Fig. 8). This can be understood from the
fact that at these low Ericksen numbers the effective viscosity depends on the
local liquid crystalline order and leads to a velocity profile that deviates
from the parabolic profile of a Newtonian fluid (see Fig. 2).

As the Ericksen number increases, the nematic host phase undergoes a
transition from the bend to the splay phase. This occurs at Ericksen numbers
$8.30<Er<18.10\;(2R/Lx=0.4),9.84<Er<21.25\;(2R/Lx=0.6)$ and
$10.37<Er<17.95\;(2R/Lx=0.8)$, respectively and is indicated by the vertical
green dashed lines in Fig. 8. The transition is accompanied by a noticeable
drop in the retardation ratio, which reaches a minimum around Ericksen numbers
$Er\simeq 20$, so just beyond the bend-to-splay transition. The minimum is
smaller the larger the particle anchoring strength is, but only for medium and
large confinement (see medium grey and green curves as well as black and red
curves in Fig. 8) and not so for small confinement (see light grey and blue
curves in Fig. 8).

Beyond Ericksen numbers in the range of $Er\simeq 20$ the retardation ratio
$v/u_{c}$ begins to increase again, giving rise to an overall non-monotonous
dependence on the Ericksen number. This is the case across all confinement
ratios, and the retardation ratios begin to flatten out towards higher
Ericksen numbers, approaching or reaching the values of Newtonian fluids
again. This non-monotonous behaviour is therefore a consequence of the
decreasing importance of liquid crystalline elasticity and consistent with the
idea that with higher the Ericksen numbers the liquid crystal behaves
rheologically more like a simple fluid.

Regarding how the retardation ratio $v/u_{c}$ depends on the particle
anchoring strength the same tends as for the minima prevail. Higher anchoring
strengths entail smaller retardation ratios unless the confinement is small.
For our largest confinement ratio $2R/L_{x}=0.8$ and strongest particle
anchoring strength $w_{part}=0.05$ we observe a very strong decrease. This,
however, originates also from the migration of the particles to the walls. The
two empty circles in Fig. 8 (and similar empty symbols at the two other
confinement ratios) permit us to estimate how the trend would continue if the
particles had been prevented from leaving the region of maximum flow velocity
at the centre of the duct.

In order to explain these findings, we have to look at several separate
mechanisms: First of all, there is the transition from the bend to the splay
state, which all particles regardless of their anchoring conditions are
subject to. The data points for vanishing particle anchoring strength
$w_{part}=0$ (light, medium and dark grey in Fig. 8) are indicative of this.
The transition causes the general reduction of the retardation ratios from
their initially approximately Newtonian values at low Ericksen numbers to
their minima around $Er\simeq 20$. The reason for this decrease is the drop in
apparent viscosity and increase in flow velocity $u_{c}$ around the centre of
the duct, whereas the regions of the particle closer to the walls act as
anchor and do not allow the particle to pick up velocity $v$ at the same
proportion.

The second mechanism at work is the reorientation of the defect ring at the
bend-to-splay transition, provided the particle anchoring strength is large
enough for a defect to emerge. The vertical orientation of the defect ring
with its ring plane perpendicular to the flow direction and walls increases
the effective particle radius in the narrowest duct dimension and therefore
the effective confinement ratio. This leads to lower retardation ratio
$v/u_{c}$ the larger the anchoring strength is. However, our results suggest
this is only the case provided the confinement is not too large. For instance,
at $2R/L_{x}=0.4$ there is very little difference between vanishing and very
strong particle anchoring up to Ericksen numbers $Er\simeq 60$, while at
$2R/L_{x}=0.6$ and $0.8$ differences are clearly visible at all Ericksen
numbers. This subtlety can be understood by realising that at the different
confinement ratios and flow velocities both velocity and order parameter
gradients differ across the particle diameter. At a given flow velocity the
gradients are largest in large confinement and vice versa. At a given
confinement ratio the velocity gradient is largest at large flow velocities
and Ericksen numbers. It is precisely this nonlinear order-flow coupling and
the interactions between flow and order structure in the vicinity of the
particle that cause the observed minor variations in the retardation ratio.

Finally, there is also the trivial possibility of a direct interaction with
the wall anchoring when Ericksen numbers and confinement ratios are large.

##  5 Conclusions

In microfluidic setups of particle suspensions, confinement is often necessary
as it allows a certain degree of lateral control over the particle positions,
for instance when techniques like confocal or polarised microscopy are used.
Our present study has primarily the goal to address some knowledge gaps as to
how defects influence and alter the migration behaviour of colloidal particles
in moderate and large confinement.

In homeotropic anchoring conditions at the walls and surface of the particle
the director field is in the H- or bend state at low Ericksen numbers and has
a Saturn ring defect which is oriented parallel to the walls. Increasing the
confinement changes the appearance of the defect ring downstream. It can
either thicken the defect ring, or invoke a migration to off-centre off-wall
positions which we identify with the weak attractor region in our previous
work 27. The latter entails to a slight distortion and tilt away from the
centre plane.

At moderately high Ericksen numbers around $Er\simeq 20$ we observe the
transition from the bend or H-state to the splay or V-state. This leads
generally to a reorientation of the defect ring with a ring plane
perpendicular to the flow direction and walls. The defect ring is slightly
peeled off downstream in the confined dimension, but sits at the particle’s
equatorial region in the unconfined dimension, giving it the appearance of an
open mouth when viewed from the flow direction. These features are retained at
higher Ericksen numbers and in lower confinement. Highly confined particles
show a strong tendency to migrate to the walls, a behaviour we observed also
in our previous work 27. This leads to a highly asymmetric defect and induces
a satellite region of low order upstream of the particle, which acts partly as
an obstacle and forces the flow to divert around it.

The interaction between nematic order and flow on one hand, and the fluid-
solid interaction on the other hand results in a non-monotonous dependence of
the retardation ratio, the ratio of particle migration velocity to the maximum
velocity at the centre of the duct, on the Ericksen number. When the Ericksen
number is low, the retardation ratio is close to values observed in a
Newtonian host phase in all confinement ratios and particle anchoring
conditions. This is also the case for vanishing or low anchoring strength and
at high Ericksen numbers, where the nematic liquid crystal behaves
increasingly like a simple Newtonian fluid as the relative importance of
elastic effects decreases. Intermediate Ericksen numbers, however, are
characterised by a pronounced minimum in the retardation ratio. We attribute
this to a combination of two effects: Firstly there is the bend-to-splay
transition, to which particles in all anchoring conditions are subject.
Secondly, the defect ring undergoes a reorientation from horizontal alignment
with the ring plane parallel to the walls to a vertical orientation, which has
the ring plane perpendicular to the flow direction and the walls. This
increases the effective particle radius and therefore the confinement. The
second effect is only present when the defect ring is properly formed, i.e.
for stronger particle anchoring strengths, and when the confinement is lower.
This is because the increased retardation that the particle experiences is a
consequence of the interaction of the defect with the gradients of the flow
velocity and liquid crystalline order.

The present study leaves some questions untouched, for example how planar
degenerate or hybrid anchoring conditions affects the defect morphology and
migration behaviour of the particles. Planar degenerate wall anchoring is
fully compatible with the flow alignment that takes place at higher Ericksen
numbers. Hence, there is no bend-to-splay transition, rather a more gradual
transition to a state where the director field is flow-aligned at the Leslie
angle. Furthermore, planar degenerate anchoring conditions on the particle
surface invoke topologically different boojum defects that occur at low
Ericksen numbers symmetrically upstream and downstream of the particle on an
axis that goes through the centre of the particle. Similarly, particles with
homeotropic anchoring conditions as used in the present study, but larger
diameters, have also topologically different defects. Instead of the half-
integer bulk defect loops with topological charge $-1/2$ they have dipolar
full integer satellite defects with topological charge $-1$ that sit in a
distance from the particle surface. It is not clear how these topological
differences would affect the advection and migration behaviour.

## Conflicts of Interest

There are no conflicts to declare.

## Acknowledgements

We would like to thank Timm Krüger and Uroš Tkalec for very helpful
discussions. We acknowledge support from EPSRC under Grant No. EP/R513349/1.
M.L. acknowledges funding from the Mac Robertson Postgraduate Travel
Scholarship. For the purpose of complying with UKRI’s open access policy, the
authors have applied a Creative Commons Attribution (CC BY) license to any
Author Accepted Manuscript version arising from this submission. This work
used the ARCHIE-WeSt High Performance Computer (www.archie-west.ac.uk) based
at the University of Strathclyde.

## Notes and references

  * Chin _et al._ 2011 C. D. Chin, T. Laksanasopin, Y. K. Cheung, D. Steinmiller, V. Linder, H. Parsa, J. Wang, H. Moore, R. Rouse, G. Umviligihozo, E. Karita, L. Mwambarangwe, S. L. Braunstein, J. van de Wijgert, R. Sahabo, J. E. Justman, W. El-Sadr and S. K. Sia, _Nat. Med._ , 2011, 17, 1015–U138.
  * Hu _et al._ 2014 J. Hu, S. Wang, L. Wang, F. Li, B. Pingguan-Murphy, T. J. Lu and F. Xu, _Biosens. Bioelectron._ , 2014, 54, 585–597.
  * Zhang _et al._ 2013 Y. Zhang, H. F. Chan and K. W. Leong, _Adv. Drug Deliv. Rev._ , 2013, 65, 104–120.
  * Golombek _et al._ 2018 S. K. Golombek, J.-N. May, B. Theek, L. Appold, N. Drude, F. Kiessling and T. Lammers, _Adv. Drug Deliv. Rev._ , 2018, 130, 17–38.
  * Seemann _et al._ 2012 R. Seemann, M. Brinkmann, T. Pfohl and S. Herminghaus, _Rep. Prog. Phys._ , 2012, 75, 016601.
  * Mark _et al._ 2010 D. Mark, S. Haeberle, G. Roth, F. von Stetten and R. Zengerle, _Chem. Soc. Rev._ , 2010, 39, 1153–1182.
  * Qin _et al._ 2010 D. Qin, Y. Xia and G. M. Whitesides, _Nat. Protoc._ , 2010, 5, 491–502.
  * Sackmann _et al._ 2014 E. K. Sackmann, A. L. Fulton and D. J. Beebe, _Nature_ , 2014, 507, 181–189.
  * Squires and Quake 2005 T. M. Squires and S. R. Quake, _Rev. Mod. Phys._ , 2005, 77, 977–1026.
  * Di Carlo _et al._ 2007 D. Di Carlo, D. Irimia, R. G. Tompkins and M. Toner, _Proc. Natl. Acad. Sci. U. S. A._ , 2007, 104, 18892–18897.
  * Amini _et al._ 2014 H. Amini, W. Lee and D. Di Carlo, _Lab Chip_ , 2014, 14, 2739\.
  * Zhang _et al._ 2016 J. Zhang, S. Yan, D. Yuan, G. Alici, N.-T. Nguyen, M. Ebrahimi Warkiani and W. Li, _Lab Chip_ , 2016, 16, 10–34.
  * Ho and Leal 1974 B. Ho and L. Leal, _J. Fluid. Mech._ , 1974, 65, 365–400.
  * Ganatos _et al._ 1980 P. Ganatos, R. Pfeffer and S. Weinbaum, _J. Fluid. Mech._ , 1980, 99, 755–783.
  * Staben _et al._ 2003 M. E. Staben, A. Z. Zinchenko and R. H. Davis, _Phys. Fluids_ , 2003, 15, 1711–1733.
  * Staben and Davis 2005 M. E. Staben and R. H. Davis, _Int. J. Multiph. Flow_ , 2005, 31, 529–547.
  * Nikoubashman _et al._ 2013 A. Nikoubashman, C. N. Likos and G. Kahl, _Soft Matter_ , 2013, 9, 2603–2613.
  * Stark and Ventzki 2002 H. Stark and D. Ventzki, _Europhys. Lett._ , 2002, 57, 60–66.
  * Čopar _et al._ 2021 S. Čopar, M. Ravnik and S. Žumer, _Crystals_ , 2021, 11, 956.
  * Urbanski _et al._ 2017 M. Urbanski, C. G. Reyes, J. Noh, A. Sharma, Y. Geng, V. S. R. Jampani and J. P. F. Lagerwall, _J. Phys. Condens. Matt._ , 2017, 29, 133003\.
  * Mondal _et al._ 2018 S. Mondal, I. M. Griffiths, F. Charlet and A. Majumdar, _Fluids_ , 2018, 3, 39.
  * Kos and Ravnik 2020 Ž. Kos and M. Ravnik, _Sci. Rep._ , 2020, 10, 1446.
  * Fedorowicz _et al._ 2023 K. Fedorowicz, R. Prosser and A. Sengupta, 2023, arXiv:2304.02759.
  * Mondal _et al._ 2018 S. Mondal, A. Majumdar and I. M. Griffiths, _J. Colloid Interface Sci._ , 2018, 528, 431–442.
  * Stieger _et al._ 2014 T. Stieger, M. Schoen and M. G. Mazza, _J. Chem. Phys._ , 2014, 140, 054905.
  * Stieger _et al._ 2016 T. Stieger, S. Püschel-Schlotthauer, M. Schoen and M. G. Mazza, _Molecular Physics_ , 2016, 114, 259–275.
  * Lesniewska _et al._ 2022 M. Lesniewska, N. Mottram and O. Henrich, _Soft Matter_ , 2022, 18, 6942–6953.
  * de Gennes and Prost 1993 P.-G. de Gennes and J. Prost, _The Physics of Liquid Crystals_ , Oxford University Press, New York, 1993.
  * Wright and Mermin 1989 D. C. Wright and N. D. Mermin, _Rev. Mod. Phys._ , 1989, 61, 385–432.
  * Beris and Edwards 1994 A. N. Beris and B. J. Edwards, _Thermodynamics of flowing systems_ , Oxford Univ. Press, 1994.
  * Škarabot _et al._ 2007 M. Škarabot, M. Ravnik, S. Žumer, U. Tkalec, I. Poberaj, D. Babič, N. Osterman and I. Musevič, _Phys. Rev. E_ , 2007, 76, 1–8.
  * Marenduzzo _et al._ 2007 D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, _Phys. Rev. E_ , 2007, 76, 1–18.
  * Ladd 1994 A. Ladd, _J. Fluid. Mech._ , 1994, 271, 285.
  * Ladd 1994 A. Ladd, _J. Fluid. Mech._ , 1994, 271, 311.
  * Krüger _et al._ 2017 T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva and E. M. Viggen, _The Lattice Boltzmann Method: Principles and Practice_ , Springer International Publishing, Switzerland, 2017.
  * Parodi 1970 O. Parodi, _J. Phys._ , 1970, 31, 581–584.
  * 37 Ludwig GitHub repository, https://github.com/ludwig-cf/ludwig.
  * Desplat _et al._ 2001 J. C. Desplat, I. Pagonabarraga and P. Bladon, _Comput. Phys. Commun._ , 2001, 134, 273–290.
  * Adhikari _et al._ 2005 R. Adhikari, K. Stratford, M. E. Cates and A. J. Wagner, _Europhys. Lett._ , 2005, 71, 473–479.