# Optimum dissipation by cruising in dense suspensions

Pappu Acharya<EMAIL_ADDRESS>Theoretical Chemistry, Lund
University, Lund SE-221 00, Sweden  Martin Trulsson
<EMAIL_ADDRESS>Theoretical Chemistry, Lund University, Lund
SE-221 00, Sweden

###### Abstract

Dense suspensions tend to shear jam at large packing fractions. However, it
has recently been shown that various oscillation protocols can unjam as well
as reduce viscosity and dissipation. In this numerical work, we ask ourselves
what is the optimum shear protocol in terms of dissipation. We show that many
cruising protocols’ dissipation are similar to shear protocols with a steady
primary shear and superimposed cross oscillations, even though the latter’s
viscosity reduction is more considerable. Furthermore, we find that
alternating between primary and perpendicular oscillations yields a much
higher dissipation than the two protocols mentioned above, yet has similar
viscosity as the cross-oscillatory one. While self-organization has been shown
to minimize viscosity, our findings challenge the idea that random
organization is the underlying mechanism for reducing dissipation. Instead,
shear “fragility” combined with geometry seems to be the key ingredients,
which explains the counter-intuitive decoupling of the minima of viscosity and
dissipation for the cruising protocol. This work paves the way for a new class
of highly-energy efficient flow protocols.

##  I Introduction

Dense suspensions are widespread in our society, both as industrial and
consumer products. Studying them is, hence, relevant for a better society. For
example, medical pills formation, vaccines, and manufacturing of ceramic and
paste all rely on processes that involve dense suspensions. Dense suspensions
are also widespread in nature, from flowing magma to beach sand and mud.
Mixing cornstarch and water, creating oobleck, is also an ideal and popular
toy experiment to increase interest and educate our kids in natural science.

Dense suspensions are composed of solid, Brownian, or non-Brownian, particles
suspended in a fluid with a high volume fraction of the solid phase. As the
fraction of solid particles increases, so does also the flow resistance,
_i.e.,_ the viscosity. Above a certain packing fraction $\phi_{J}$, a
suspension generally jams and stops to flow [1]. The viscosity usually shows
an inverse power law divergence as a function of the distance to this jamming
point, indicating a second-order non-equilibrium phase transition with an
associated diverging cooperative length [2, 3, 4]. Dense suspensions also
often show non-Newtonian behavior, where the flow resistance is nonlinear as a
function of the flow rate. A sub-linear relationship is denoted as shear-
thinning and super-linear shear-thickening. A typical example of a shear-
thinning suspension is blood and a shear-thickening is oobleck. The latter can
even show discontinuous shear-thickening (denoted DST) [5, 6, 7, 8] where
above a certain shear stress or flow rate a suspension does an abrupt
transition from a low-viscous fluid to a high viscous one or even become an
elastic solid. There are advantages and disadvantages of such a DST scenario.
An application making use of one of these advantages is Kevlar soaked with
such a DST suspension to generate a flexible yet impact-proof space suit,
called an STF-Armor [9]. The same concept has also been used to design bullet-
proof jackets [10]. Despite the attractive fringe benefit of shear thickening,
the main drawback comes when a suspension needs to flow.

For all dense suspensions, DST or not, it is beneficial if we can alter and
facilitate flowability. In the literature, many works have been focused on the
rheological properties of dense suspensions and some are able to reduce a
suspension’s viscosity by more than 90 percent. The procedure to reduce
viscosity involves many strategies, _e.g._ adding polydispersity [11] to the
particle sizes, changing their shapes [12, 13], adding some lubricating
particles, by adding another fluid [14] _etc._. One recent and successful
approach is to apply a cross shear to the suspension as it does not involve a
change of a suspension’s composition. However, there is an extra energy cost
involved in doing so. While many of the previous works [15, 16, 17] have been
focused on reducing the viscosity, only a few focused on the energy needed to
propagate the suspension, _i.e.,_ the energy dissipation, [16, 18].

In this work, we design and study a large set of numerical shear protocols to
find optimal ways of propagating a suspension in terms of energy dissipation
(per strain and time). We argue that oscillations are not necessary for the
dissipation reduction, rather, a change in the direction of the fluid flow is
sufficient, related to the fragility of dense suspensions [19, 20, 21]. We
illustrate this by applying a cruising protocol, which gives a comparable
dissipation reduction compared to previously best-reported protocols. This
cruising protocol uses only two directions, compared to three for oscillating
protocols, and does not involve a flow reversal [22, 23], which are known to
reduce the viscosity, but rather shear rotations. Furthermore, we show that
the optimal angle for our cruising protocol, _i.e.,_ shear rotations, depends
on the packing fraction, where in the vicinity of the steady-state jamming
point obtuse angles $(>90^{o})$ are superior. We give a clear explanation of
why some protocols are better than others with semi-theoretical expressions of
the dissipation for each protocol.

Figure 1: Schematic diagram of all four protocols (a) AO or alternating
oscillatory cross shear [18], (b) HO or half oscillatory cross shear, (c) SW
or square wave shear and (d) CR or cruising. In all four protocols, we used
$xz$ as the primary shear plane while $xy$ as the cross shear plane.

##  II the AO protocol’s dissipation

We start by re-analyzing one of the previously proposed protocols, denoted AO
(Alternating Oscillatory cross shear) [18]. In this protocol, the system is
sheared linearly in the primary direction while an oscillatory shear is
applied along the perpendicular direction. A full cycle is divided into two
parts where $\alpha$ is related to the time spent in the primary and
oscillatory directions. A strain rate of $\dot{\gamma}_{\rm
pri}=\dot{\gamma}/(1-\alpha)$ is applied in the primary direction for time
$(1-\alpha)T$, yielding a strain in the primary direction equal to
$\dot{\gamma}_{\rm pri}(1-\alpha)T=\dot{\gamma}T=\gamma_{\rm pri}$. Here,
$\dot{\gamma}$ is the corresponding shear rate for the steady-shear state,
yielding the same strain if sheared for a time $T$. During shearing in the
primary direction, there is no oscillatory shear along the perpendicular
direction. Then, for a time $\alpha T$, an oscillatory shear with a time
period of $\alpha T$ is applied as $\gamma_{\rm cross}=\gamma_{\rm
osc}\sin\omega t$, where $\omega=2\pi/\alpha T$ and $t$ is time. We find the
equation for dissipation per strain (and time) with $\alpha$ (see the Appendix
for details) as

$W=\frac{W_{\parallel}+W_{\perp}}{W_{\rm
ss}}=\frac{\eta_{\parallel}}{\eta_{\rm
ss}}\left(\frac{1}{1-\alpha}\right)+\frac{\eta_{\perp}}{\eta_{ss}}\frac{2\pi^{2}}{f^{2}}\left(\frac{1}{\alpha}\right),$
(1)

where $W_{\parallel}$, $W_{\perp}$, $W_{\rm SS}$, are the energy of
dissipation per strain for the primary shear, cross shear, and steady-state
(SS) respectively. Here, $\eta_{ss}$ is the steady state viscosity which only
depends upon the packing fraction of the system, $f=\gamma_{\rm
pri}/\gamma_{\rm osc}$ is the ratio of strain along the primary and
perpendicular directions, while, $\eta_{\parallel}$ and $\eta_{\perp}$ are the
viscosities along the primary and perpendicular directions and can be computed
directly from the simulations. We plot the dissipation from the simulations as
a function of $\alpha$ along with this expression in Fig. 6(a) in the
Appendix, which shows a good agreement for three different $f$. As we apply
primary shear to the suspension along the flow direction of the $xz$ plane, we
denote the primary direction of flow as $xz$ and the oscillatory direction as
$xy$ throughout this paper, therefore $\eta_{\parallel}$ and $\eta_{\perp}$
can be defined as

$\displaystyle\eta_{\parallel}$ $\displaystyle=$
$\displaystyle\frac{\sigma_{xz}}{\dot{\gamma}_{\rm pri}},$ (2)
$\displaystyle\eta_{\perp}$ $\displaystyle=$
$\displaystyle\frac{\int_{0}^{\frac{2\pi}{\omega}}\sigma_{xy}\cos(\omega
t)dt}{\dot{\gamma}_{0}\int_{0}^{\frac{2\pi}{\omega}}\cos^{2}(\omega t)dt},$

where $\sigma_{kl}$ denotes shear stresses. Note that the viscosities are
rate-independent in the over-damped limit and hence do not depend on $\alpha$.
The minimum point for each $f$ can be calculated by taking a derivative of Eq.
(1) with respect to $\alpha$ and equating it to zero. The real root for all
the values of $\eta_{\parallel}$ and $\eta_{\perp}$ can then be written as

$\alpha=\frac{1}{1+\frac{f}{\sqrt{2}\pi}\left(\frac{\eta_{\parallel}}{\eta_{\perp}}\right)}.$
(3)

At increasing $f$, $\eta_{\parallel}$ increases, while $\eta_{\perp}$ is
roughly constant. Hence, we can conclude that the minimum point in the
dissipation curve will shift toward $\alpha=0$ with increasing $f$ as
displayed in Fig. 6(a) in the Appendix. The minimum value of the dissipation
for different $f$ can then be written as

$W_{min}=\frac{(\sqrt{\eta_{\parallel}}+\frac{\sqrt{2}\pi}{f}\sqrt{\eta_{\perp}})\^{A}{{}^{2}}}{\eta_{ss}}.$
(4)

One of our first findings is that even though the AO protocol reduces the
viscosity in the primary direction, it comes at a cost in dissipation, where
the AO dissipate more than doing standard steady-state. This is at odds with
the original AO results of Ref. [18], which was due to a subtle error (see the
reference’s erratum).

A similar protocol but with a constant shear rate along the perpendicular
direction yields the dissipation as (see the Appendix for details)

$W_{\rm AO}=\frac{W_{\parallel}+W_{\perp}}{W_{\rm
ss}}=\frac{\eta_{\parallel}}{\eta_{\rm
ss}}\left(\frac{1}{1-\alpha}\right)+\frac{\eta_{\perp}}{\eta_{\rm
ss}}\frac{16}{f^{2}}\left(\frac{1}{\alpha}\right).$ (5)

Both protocols’ viscosities are very close with the latter slightly less
dissipative as $2\pi^{2}>16$. We again plot the dissipation from the
simulations with $\alpha$ along with this expression in Fig. 6(b) in the
Appendix, which shows a good agreement for three different $f$. Nevertheless,
both protocols dissipate a large amount of energy (more than the steady state)
and are not very fruitful in terms of energy reduction.

We continue by re-analyzing also the protocol denoted SO (Simple oscillatory
cross shear), where oscillations instead are superimposed on a primary shear
direction in a steady manner. In contrast to the AO protocol, we find the same
dissipation (data not shown) from our numerical simulations as reported in
Ref. [18]. For expressions of how dissipation depends on viscosities, see our
Appendix. Given that SO and AO give such different dissipation, this naturally
leads to the questions of (i) why AO and SO are so different? and (ii) what is
the optimum way of propagating a suspension? We focus on relative dissipation,
keeping the strain and time compared to the steady-state protocol, _i.e.,_ at
a fixed global shear rate $\dot{\gamma}=\gamma_{\rm pri}/T$. This as
dissipation per strain $W$ can, trivially, be reduced by allowing for a slower
global shear rate (equivalent to longer shear time). This can for example be
realized as $W_{SS}=\eta_{\rm ss}\dot{\gamma}$ and $\eta_{\rm ss}$ is rate-
independent in our model. This means we trade time for energy. In practice,
one wants to keep both at a minimum. Since $W_{X}=\eta_{X}\dot{\gamma}$
irrespective of the protocol, the relative dissipation (per strain and time)
_i.e._ $W=W_{X}/W_{SS}$ is actually measuring a relative “dissipative”
viscosity [16].

Figure 2: Relative viscosity along the (a) parallel and (b) perpendicular
directions for different protocols with packing fraction. We also plot the
steady state value as a comparison. In all the simulations, we use $\mu=1$ and
$\alpha=0.8,0.7,0.5$ and $0.5$ for AO, HO, SW, and CR respectively. Here the
symbols are numerical data while the dashed lines are guides for the eye.

##  III Improving the AO protocol

The dissipation in the AO protocol is rather high despite the huge viscosity
reduction. To minimize energy dissipation, we systematically apply different
types of protocols to our system and ascertain an optimum protocol. Here we
give an introduction to various protocols. For all the protocols, we use
strains of $1\%$ in all directions and normalize the dissipation with the
propagating strain. We then normalize with the steady-state value at the same
global shear rate as discussed earlier.

###  III.1 Half oscillation

A full oscillation along the perpendicular direction costs an immense amount
of energy; therefore, we first try to reduce that by applying half oscillation
(see Fig. 1 (b)). The idea is to reduce the accumulated strain in the cross
direction from $4\gamma_{\rm osc}$ to $2\gamma_{\rm osc}$. In this half
oscillation protocol, a full cycle is composed of $1\%$ strain along the
primary direction (flow direction of the $xz$ plane), with a shear rate
$\dot{\gamma}_{\rm pri}=\dot{\gamma}/(1-\alpha)$ and a time $(1-\alpha)T$.
Along the perpendicular direction _i.e._ $xy$, we shear the system $1\%$ up
and $1\%$ down with the applied shear rate $\dot{\gamma}_{\rm
osc}=2\dot{\gamma}/\alpha$ for a time $\alpha T$. The total time for a cycle
is hence $T$ and the strain per cycle in the primary direction is $1\%$ just
like in the previously described AO protocol. For the half oscillation, the
relative dissipation relation can be written as (see the Appendix for details)

$W_{\rm HO}=\frac{W_{\parallel}+W_{\perp}}{W_{\rm
ss}}=\frac{\eta_{\parallel}}{\eta_{\rm
ss}}\left(\frac{1}{1-\alpha}\right)+\frac{\eta_{\perp}}{\eta_{\rm
ss}}\frac{4}{f^{2}}\left(\frac{1}{\alpha}\right).$ (6)

Comparing with Eq. (5), we see that we only change the prefactor of the second
term from $16$ to $4$. However, both $\eta_{\parallel}$ and $\eta_{\perp}$
depend on the protocol applied. Similar to the full oscillation, once the two
viscosities are known, we can easily minimize the dissipation with respect to
$\alpha$. All reported relative dissipation corresponds to the minimum values.
The dissipation is slightly less in this case compared to the AO, although the
viscosity reduction is lower in the former. We discuss these results later in
section V.

###  III.2 Square wave

In order to reduce the energy dissipation further, we design a new protocol
where the shear path looks like a square wave (SW), as shown in Fig. 1(c).
This could be considered as a quarter oscillation protocol but differs from
the previous two in that it does not include a shear reversal between two
primary shears. In the first cycle, we shear the system along the flow
direction of the $xz$ plane in the first half, and then along the flow
direction of the $xy$ plane in the second half. In the second cycle, however,
we shear the system along the flow direction of the $xz$ plane again in the
first half, but in the second half, we shear it along $-xy$. The key idea with
this is to further reduce the accumulated strain in the cross direction (per
cycle). We discuss this protocol’s viscosity and dissipation reductions later
in section V. The relative dissipation for a square wave follows a similar
trend as for the half oscillation protocol, where we once again affect
prefactor of the second term, and can be written as (see the Appendix for
details)

$W_{\rm SW}=\frac{\eta_{\parallel}}{\eta_{\rm
ss}}\left(\frac{1}{1-\alpha}\right)+\frac{\eta_{\perp}}{\eta_{\rm
ss}}\left(\frac{1}{\alpha}\right).$ (7)

As for the AO and HO protocols, we report the minimum dissipation values (with
respect to $\alpha$).

Figure 3: Relative dissipation for different protocols as a function of
packing fraction. We systematically decrease the dissipation by changing the
protocol, most importantly by decreasing the accumulated strain in the
perpendicular direction, and conclude that cruising is better than the other
protocols. Symbols correspond to numerical data using $\mu=1$, where
viscosities from Fig. 2 have been used to find the corresponding minimum
dissipation. Dashed lines are merely guides to the eye.

###  III.3 Cruising protocol

Although the square wave protocol reduces the viscosity by a large amount (see
section V), it will not propagate the suspensions when in cross mode, as the
shearing along $xy$ and $-xy$ cancels any such propagation. The simplest way
to overcome this issue is by continuously shear along one of these two
opposite directions. This leads to cruising protocol (CR), where, in the first
half of the cycle, we shear the system along the flow direction of the $xz$
plane and in the second half along the flow direction of the $xy$ plane as
shown in Fig. 1(d). In contrast to all the above protocols, one main advantage
of this protocol is that we propagate the suspension for each direction. In a
full cycle, the ratio between the accumulated strain and the strain in the
“propagating” direction equals $\sqrt{2}$ (compared to $\geq 2$ in the other
protocols). Although the viscosity reduction here is the lowest among all the
protocols, because of the highly efficient propagation, this is the most
beneficiary in terms of dissipation per unit strain (and at a fixed global
shear rate). We discuss the results later in the section V and make it clear
that this protocol is best among all the protocols we discussed and equivalent
to the best existing protocol called SO [18]. The dissipation is related to
the viscosities as

$W_{\rm CR}=\frac{1}{2}\left[\frac{\eta^{\prime}}{\eta_{\rm
ss}}\left(\frac{1}{(1-\alpha)}\right)+\frac{\eta^{\prime\prime}}{\eta_{\rm
ss}}\left(\frac{1}{\alpha}\right)\right].$ (8)

In the above equation, we have relabeled the viscosities as $\eta^{\prime}$
and $\eta^{\prime\prime}$ as none of them is now parallel or perpendicular to
the propagating axis. Although $W_{\rm CR}$ looks like $2$ times smaller than
$W_{\rm SW}$, it is important to understand that the viscosities depend on the
protocol. Since we have $\eta^{\prime}=\eta^{\prime\prime}=\eta$ in this
protocol, the minimum dissipation occurs at $\alpha=0.5$. The minimum
dissipation hence takes the simple form: $W_{\rm CR}=2\eta/\eta_{\rm ss}$.
This generally means that we need at least a $50\%$ reduction of the viscosity
in order to be a viable approach.

###  III.4 Relaxing the turn angle

For all protocols, we have so far exclusively used $\pm 90^{o}$ turns with
respect to the $xz$ plane when changing the flow direction. One can easily
relax this constraint and allow for any degree of turns, with $\pm(1/2)\theta$
degrees of turns around the propagating direction. The dissipation for the CR
protocol and using a particular angle $\theta$ can be written as (see the
Appendix for details)

$W_{\rm
CR}(\theta)=\frac{1}{(2+2\cos\theta)}\left[\frac{\eta^{\prime}}{\eta_{\rm
ss}}\left(\frac{1}{1-\alpha}\right)+\frac{\eta^{\prime\prime}}{\eta_{\rm
ss}}\left(\frac{1}{\alpha}\right)\right].$ (9)

Note that it will yield Eq. (8) at an angle $\theta=90^{0}$. The reduced form,
using $\alpha=0.5$ and where $\eta^{\prime}=\eta^{\prime\prime}=\eta$, is then

$W_{\rm CR}(\theta)=\frac{2}{(1+\cos\theta)}\frac{\eta}{\eta_{\rm ss}}$ (10)

which yields $W_{\rm CR}=1$ at $\theta=0^{o}$, as $\eta=\eta_{\rm ss}$ in that
case. This perfectly agrees with the idea that cruising has now become a
steady-shear protocol.

##  IV  Simulation details

The only unknowns in the above dissipation equations are the viscosities,
which can be obtained by numerical simulations. We simulate a dense athermal
suspension with various protocols and measure the energy dissipation per
strain and time. The latter is equivalent to having a fixed global shear rate
in the transport direction. In this work, we study a non-inertial suspension
composed of bi-disperse particles. The particles are iso-dense with the fluid
and confined in a three-dimensional periodic box using Lees-Edwards boundary
conditions at a packing fraction $\phi$. We use a bi-disperse $50:50$ mixture
of particles with a diameter ratio of $1:1.4$. A particle undergoes three
types of force and torque: Stokes drag, contact, and lubrication forces. The
Stokes drag force and torque are given by,

$\displaystyle\vec{F}_{i}^{s}$ $\displaystyle=$
$\displaystyle-6\pi\eta_{f}a_{i}(\vec{v}_{i}-\vec{v}_{f}),$
$\displaystyle\vec{\tau}^{s}_{i}$ $\displaystyle=$
$\displaystyle-8\pi\eta_{f}a_{i}^{3}(\vec{\omega}_{i}-\vec{\omega}_{f}).$ (11)

Here, $a_{i}$, $v_{i}$, and $\omega_{i}$ are the radius, linear velocity, and
angular velocity of particle $i$, while $v_{f}$ and $\omega_{f}$ are the
fluid’s linear and angular velocity (vorticity) and $\eta_{f}$ is the
viscosity of the fluid. We chose parameters such that the Stokes number
$\rho\dot{\gamma}\langle a\rangle^{2}/\eta_{f}\ll 1$ for $\alpha$ values in
the range [0.2,0.8], where $\langle a\rangle$ denotes the average radius and
$\rho$ the density of both solid particles and the fluid.

The contact force between particles is described by harmonic springs

$\vec{F}_{ij}=k_{n}\delta^{ij}_{n}\hat{n}^{ij}+k_{t}\delta^{ij}_{t}\hat{t}^{ij}.$
(12)

Here $\delta^{ij}_{n}$ is the normal overlap between two particles while
$\delta^{ij}_{t}$ is the relative tangential displacement. The springs’ normal
and tangential stiffness constants are denoted by $k_{n}$ and $k_{t}$. We
choose $k_{t}=(2/7)k_{n}$, and $k_{n}$ and $k_{t}$ values large enough such
that our particles full-fill near-hard sphere conditions (typically
$k_{n}/P\sim 10^{4}$). The unit vector along the normal and the tangential
direction of two particles are denoted by $\hat{n}^{ij}$ and $\hat{t}^{ij}$.
The tangential forces for each contact are restricted by the Coulomb friction
criteria $|\vec{F}^{t}_{ij}|<|\mu\vec{F}^{n}_{ij}|$, where, $\mu$ is the
coefficient of friction. When a lubricating film separates two particles, the
interaction force acting along normal and tangential directions can be written
as

$\displaystyle F_{ij}^{l,n}(h_{ij})$ $\displaystyle=$
$\displaystyle\frac{3}{2}\pi\eta_{f}a_{ij}\frac{(\vec{v}_{i}-\vec{v}_{j})\cdot\hat{n}_{ij}}{h_{ij}+\delta},$
$\displaystyle F_{ij}^{l,t}(h_{ij})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\pi\eta_{f}\ln\left(\frac{2a_{ij}}{h_{ij}+\delta}\right)(\vec{v}_{j}-\vec{v}_{i})\cdot\hat{t}_{ij},$
(13)

where $a_{ij}$ is the effective radius defined as
$a_{ij}=a_{i}a_{j}/(a_{i}+a_{j})$ while $h_{ij}$ is the gap between two
particles. We choose a lower cut-off $\delta=0.001a_{s}$ (where $a_{s}$ is the
radius of the smaller particle) to avoid divergence at $h_{ij}=0$. In our
simulations, first, we create a packing with no overlap between particles.
After applying various types of shear on such a system, we perform basic
molecular dynamic simulations and solve the trajectories using LAMMPS [24].
All quantities are for steady-state conditions, where we sheared the
suspension for at least $2$ unit of accumulated strain before we start to
measure. Data is then collected for a further $18$ unit of accumulated
strains. Dissipation in various directions are measured by
$W^{\prime}=\int\mathbf{\sigma}:\mathbf{\dot{\gamma}}\,dt$. Unless specified,
we strain each direction $1\%$ before we change direction (the exception being
the AO, where we have an intermediate straining of $2\%$ for the intermediate
part of the oscillation).

##  V Results

In this section, we discuss our numerical results. We focus on both the
viscosity reduction as well as the dissipation for tested protocols.

###  V.1 Viscosity Reduction

We begin by analyzing the viscosity reduction for all the protocols. Fig. 2(a)
and (b) show how the viscosities, in parallel and perpendicular directions,
vary with the packing fraction for our studied protocols respectively. In the
figure, we also plot the steady shear values of the viscosity for comparison.
The figure shows that all protocols successfully reduce the viscosity at a
fixed packing fraction. The largest reduction is found for the AO scheme,
while the least efficient turns out to be the case of the cruising. While the
AO and CR protocols have roughly the same viscosities in both directions, as
assumed in Ref. [18], this is not a universal result, as both HO and SW
protocols clearly show a lower perpendicular viscosity compared to their
parallel ones. In general, the more reversals (in the perpendicular direction)
we do between the primary strains, the larger the viscosity reduction. This
can be seen as AO corresponds to two reversals, HO to one, SW a half, and CR
to zero. This agrees with Ref. [18]’s findings.

Interestingly, the SW and CR show very similar viscosities in the “parallel”
direction but different in the “perpendicular” one 111The notation of parallel
and perpendicular in the cruising protocol is a bit misleading, as none of
them is along the propagating axis. Nevertheless, treating them on the same
footing is useful.. Both protocols have the same accumulated strain in the
perpendicular direction. The only difference is that the SW involves
alternating directions (a half-shear reversal per cycle) along the
perpendicular direction. Hence, a $1\%$ strain in the primary direction is not
sufficient to lose the micro-structural memory built up during cross shears,
explaining why we see a difference between these SW and CR viscosities.

Figure 4:  (a) Relative viscosity as a function of packing fraction for the
cruising protocol at different cruising angles. It is clear that the viscosity
reduction is larger for a wider angle, and it will hit a maximum at $180^{o}$.
(b) Relative dissipation as a function of packing fraction in the cruising
protocol for different cruising angles. The figure clearly shows that the
angles at which the minimum dissipation occurs are different for different
packing fractions. All the data are from direct simulation, using $\mu=1$ and
$\alpha=0.5$.  Figure 5: Comparison between the relative dissipation of SO and
CR as a function of re-scaled packing fraction. The collapse indicates that
both protocols are very similar to each other, and therefore ”random
organization” is not necessary to minimize the dissipation. (Inset) The angle
of minimum dissipation for different packing fractions. For the lowest
$\phi$’s the relative dissipation is $1$ (within the noise) and marked as red
square

###  V.2  Dissipation reduction

Suspensions dissipate energy when they flow, and the motivation of this work
is to reduce the dissipation per unit strain at a fixed global shear rate. We
here study the relative dissipation, _i.e._ we divide our dissipation with the
corresponding steady-shear value at the same strain and global shear rate. A
relative dissipation below one indicates an energy saving when transporting
the suspension. Fig. 3 shows the relative dissipation as a function of packing
fraction for our tested protocols. Non-intuitively the viscosity and
dissipation reduction are, in our cases, inversely correlated. The AO
(corresponding to a full oscillation), which gives more than $90\%$ reduction
in viscosity, turns out to be less energy efficient than steady simple shear.
This is at odds with the original results of Ref. [18], but it turned out that
there were some subtle mistakes in that paper when calculating the dissipation
for the AO protocol. The inefficiency of the AO comes from the fact that an
extra oscillation is applied, which costs a huge amount of energy without
propagating the suspension. Only when we apply the CR or SW protocols we get a
substantial energy saving for an extended range of packing fractions.
Furthermore, compared to the SO protocol of Ref. [18], another protocol that
gives energy savings, we see that both SO and CR are superior in terms of
energy saving, even if we for AO, HO, and SW have chosen the $\alpha$ giving
the minimum dissipation.

###  V.3 Effect of cruising angle

We focus now on the CR protocol and relax the constraint with $90^{o}$ angles
and test a few different angles with a scanning gap of $15^{o}$ _i.e._
$\theta=15^{o},30^{o},45^{o}....$ _etc._. For example, $\theta=15^{o}$ means
for the first half of a cycle, we shear the suspension along the flow
direction of the $xz$ plane, while for the second half, we shear the system at
an angle $15^{o}$ with respect to the $xz$ plane. Fig. 4(a) shows the effect
on viscosity, where the viscosity decreases monotonically as the angle widens,
reaching its minimum value at $\theta=180^{o}$. However, as seen before, this
information might sometimes be misleading when analyzing dissipation. Fig 4(b)
shows the relative dissipation as a function of packing fraction for five
different angles. The tendency between the angles is clearly non-monotonic.
For the lowest angle, the dissipation line is fairly flat and close to one,
even though it is in general a few percent below one at intermediate to high
packing fractions. At $\theta=0^{o}$, we have $W=1$, which is consistent with
this protocol’s physical and mathematical understanding as it is the same as
the steady shear. Close to the shear jamming point for steady state, the
dissipation seems to be lowest for angles close to $120-150^{o}$. The optimum
angle varies with packing fraction, where a higher packing fraction favors
wider angles and low packing fractions narrower angles. For very wide angles,
the propagation becomes inefficient due to geometrical reasons, and the
relative dissipation per strain diverges as $\theta\to 180^{o}$. This limit
corresponds to a pure oscillatory shear. In relative dissipation, the cruising
protocol with an optimized angle is very similar to SO, and most likely within
the numerical precision, for a wide range of packing fractions as shown in
Fig. 5. In the figure, we used $\phi/\phi_{c}$, for a better collapse, where
$\phi_{c}$ is the packing fraction at which shear jamming occurs for the
steady state. This since the $\phi_{c}$’s differ, due to microscopic details
[26, 27] 222We used two different $\phi_{c}$ values for our data and the data
of Ref. [18]. Both $\phi_{c}$ values were extracted by fitting the steady
state curves in with $\eta=\eta_{0}(\phi-\phi_{C})^{-2}$. Our data yielded
$\phi_{c}=0.5900$ while the data from Ref. [18] gave $\phi_{c}=0.5827$.

##  VI Discussion

In this work, we have focused on the flow of dense suspensions and their
dissipation. We found a new highly efficient shear cruising protocol that
keeps the energy dissipation at a minimum. We found that by adapting the
cruising protocol’s angle, we can get almost a magnitude in energy gain.

Obviously, we need to comment on the success of the SO protocol. Re-analyzing
that protocol, one sees that the SO protocol’s minimum dissipation well
corresponds to a CR motion with an optimum angle. The fact that the optimum
angle varies with packing fraction (see inset of Fig. 5) corresponds well with
the data of Ref. [18], where the minimum in units of oscillation ratio
$\omega\gamma/\dot{\gamma}$ (their units) varies with packing fraction. Our
re-interpretation is that this oscillation ratio can be translated into a
cruising angle which we also show for two available packing fractions from
Ref. [18] along with our minimum angles in the inset of Fig. 5.

Our results, hence, strongly suggest that the large dissipation reduction is
not per see due to the oscillations, but rather due to a sudden change in
direction. Since we are dealing with frictional particles, a $>90^{o}$ turn
will quite successfully deactivate frictional contacts. This picture
challenges the idea that a random organization [29, 30, 31] is the primary
cause of dissipation reduction. Instead, our data suggest that dissipation
reduction (or equivalent “dissipative” viscosity) is linked to dense
suspensions’ fragility and how it varies with both packing fraction, shear
angles, and strains [19, 20, 21]. Haven only considered a few protocols; this
clearly raises the question of whether even more efficient protocols exist. We
leave that as an open challenge for the dense suspension community.

Acknowledgments: We thank Christopher Ness, Romain Mari, Joakim Stenhammar,
and Zakiyeh Yousefian for the useful discussions. We also thank Christopher
Ness for the open-access LAMMPS script. This project was funded by the Wenner-
Gren Foundation (grant number UPD2021-0147). MT is financed by the Swedish
Research Council (grant number 2021-04997). The simulations were performed on
resources provided by (storage) the Swedish National Infrastructure for
Computing (SNIC) and the center for scientific and (computing) technical
computing at Lund University (LUNARC).

## APPENDIX A: THEORETICAL CALCULATION
OF THE DISSIPATION

Figure 6:  Relative dissipation for AO protocol with (a) sinusoidal
oscillation along the perpendicular direction (b) linear oscillation along the
perpendicular direction with various primary strains $\gamma_{\rm pri}$ and
$\gamma_{\rm osc}=1\%$ fixed. Symbols are numerical data using $\mu=1$ and
$\phi=0.56$. Lines are given by Eqs. (A10) and (A13). Note that none of the
tested strains’ relative dissipation goes below one at least up to
$\phi=0.56$.

###  1 The AO protocol

The dissipation per strain can be written as

$W_{\rm Tot}=W_{\parallel}+W_{\perp}$ (A1)

When comparing with the steady-state value, the relative dissipation is,

$W=\frac{W_{\rm Tot}}{W_{SS}}=\frac{W_{\parallel}+W_{\perp}}{W_{SS}}$ (A2)

The dissipation along the primary direction in a single cycle can be written
as

$W_{\parallel}(per\,cycle)=\sum_{0}^{(1-\alpha)T}\eta_{\parallel}\dot{\gamma}^{2}_{\rm
pri},$ (A3)

where $\dot{\gamma}_{\rm pri}=\dot{\gamma}/(1-\alpha)$. Putting this in Eq.
(A3), we arrive at

$W_{\parallel}(per\,cycle)=\sum_{0}^{(1-\alpha)T}\frac{\eta_{\parallel}\dot{\gamma}^{2}}{(1-\alpha)^{2}}=\frac{\eta_{\parallel}\dot{\gamma}^{2}}{(1-\alpha)}T.$
(A4)

Therefore, the dissipation per unit strain can be written as,

$W_{\parallel}(per\,strain)=\frac{\frac{\eta_{\parallel}\dot{\gamma}^{2}}{(1-\alpha)}T}{\dot{\gamma}T}=\frac{\eta_{\parallel}\dot{\gamma}}{(1-\alpha)}.$
(A5)

The viscosity along the perpendicular direction can be defined as

$\eta_{\perp}=\frac{\int_{0}^{2\pi/\omega}\sigma(t)\cos(\omega
t)dt}{{\gamma}_{\rm osc}\omega\int_{0}^{2\pi/\omega}\cos^{2}(\omega
t)dt}=\frac{W_{\rm OSC}(per\,cycle)}{(\gamma_{\rm
osc}\omega)^{2}\int_{0}^{2\pi/\omega}\cos^{2}(\omega t)dt}.$ (A6)

Therefore, the perpendicular dissipation in a full cycle takes the form,

$W_{\rm OSC}(per\,cycle)=\eta_{\perp}(\gamma_{\rm
osc}\omega)^{2}\frac{\pi}{\omega}=\eta_{\perp}\gamma^{2}_{\rm osc}\pi\omega.$
(A7)

Putting $\gamma_{\rm osc}=\gamma_{\rm pri}=\dot{\gamma}T$ (1% in our case) we
arrive at

$W_{\rm
OSC}(per\,cycle)=\eta_{\perp}\dot{\gamma}^{2}T\frac{2\pi^{2}}{\alpha}.$ (A8)

Therefore dissipation per strain can be written as,

$W_{\perp}=\frac{2\pi^{2}\eta_{\perp}\dot{\gamma}}{\alpha}.$ (A9)

The steady-state dissipation per strain is $W_{SS}=\eta_{ss}\dot{\gamma}$.
Rewriting the Eq. (A2) by putting $W_{\parallel}$, $W_{\perp}$, $W_{SS}$ we
arrive at the final expression,

$W=\frac{W_{\parallel}+W_{\perp}}{W_{\rm
ss}}=\frac{\eta_{\parallel}}{\eta_{\rm
ss}}\left(\frac{1}{1-\alpha}\right)+\frac{\eta_{\perp}}{\eta_{ss}}\frac{2\pi^{2}}{f^{2}}\left(\frac{1}{\alpha}\right),$
(A10)

which is Eq. (1) in the main text.

When a linear oscillation is applied along the perpendicular direction, the
strain rate along the perpendicular direction needs to maintain such that in
$\alpha T$ time, it covers a total of 4% strain (1% up, 2% down and 1% up
again). Therefore, the strain rate along the perpendicular direction is
related to the global strain rate as $\dot{\gamma}_{\rm
osc}=4\dot{\gamma}/\alpha$. The dissipation along the perpendicular direction
can then be written as

$W_{\rm OSC}(per\,cycle)=\sum_{0}^{\alpha T}\eta_{\perp}\dot{\gamma}_{\rm
osc}^{2}=T\frac{16\eta_{\perp}\dot{\gamma}^{2}}{\alpha}.$ (A11)

Therefore, the dissipation per strain along the perpendicular direction takes
the form

$W_{\perp}=\frac{16\eta_{\perp}\dot{\gamma}}{\alpha}.$ (A12)

Putting this back in Eq. (A2) yields

$W_{\rm AO}=\frac{W_{\parallel}+W_{\perp}}{W_{\rm
ss}}=\frac{\eta_{\parallel}}{\eta_{\rm
ss}}\left(\frac{1}{1-\alpha}\right)+\frac{\eta_{\perp}}{\eta_{\rm
ss}}\frac{16}{f^{2}}\left(\frac{1}{\alpha}\right),$ (A13)

which is Eq. (5) in the main text. We plot the relative dissipation as a
function of $\alpha$ for both sinusoidal oscillation and linear oscillation
along the perpendicular direction in Fig. 6. The $\gamma$ in the figure legend
is the strain in a single cycle along the primary direction, while we fix the
maximum $\gamma_{osc}=1\%$ for all the curves. Therefore, $\gamma$ basically
represents $f$ as discussed in the main text. A key thing to note from the
figure is the dissipation is a little bit lower in the case of linear
oscillation, as discussed in the main text and can be visualized from Eqs.
(A10) and (A13). It should be noted that although for different $f$, the
perpendicular viscosity is similar, we used the exact values from the
individual simulations.

###  2 The CR and SW protocols

In the cruising protocol, the suspension is sheared along the flow direction
of the $xz$ plane for a time period of $(1-\alpha)T$, which yields the same
expression for the primary shear as above. As described in the main text, we
relabeled the viscosities as $\eta^{\prime}$ and $\eta^{\prime\prime}$ as none
of them is now parallel or perpendicular to the propagating axis. Therefore,
dissipation for the shearing along the primary direction in a cycle can be
written as

$W_{\rm PRI}(per\,cycle)=\frac{\eta^{\prime}\dot{\gamma}^{2}T}{1-\alpha}.$
(A14)

In a similar way, the dissipation per cycle for shearing at an angle $\theta$
can be written as

$W_{\rm SEC}(per\,cycle)=\frac{\eta^{\prime\prime}\dot{\gamma}^{2}T}{\alpha}.$
(A15)

However, unlike AO, the strain covered during a full cycle is now a vector sum
of total strain in both directions and can be written as
$\gamma=\sqrt{(2+2\cos(\theta))}\dot{\gamma}T$. The total dissipation per
strain for both components can then be written as

$W=\frac{1}{\sqrt{(2+2\cos(\theta))}}\left(\eta^{\prime}\frac{\dot{\gamma}}{1-\alpha}+\eta^{\prime\prime}\frac{\dot{\gamma}}{\alpha}\right).$
(A16)

As the global shear rate is changed, we need to divide this by the steady
state value computed by the global shear rate. Therefore $W_{SS}$ per strain
can be written as $W_{SS}=\sqrt{(2+2\cos(\theta))}\dot{\gamma}\eta_{ss}$. The
final form of dissipation in the case of cruising is then

$W_{\rm
CR}(\theta)=\frac{1}{(2+2\cos\theta)}\left[\frac{\eta^{\prime}}{\eta_{\rm
ss}}\left(\frac{1}{1-\alpha}\right)+\frac{\eta^{\prime\prime}}{\eta_{\rm
ss}}\left(\frac{1}{\alpha}\right)\right],$ (A17)

which is Eq. (9) in the main text. In the square wave protocol, we need to
consider two cycles and take an average as the second half of the first cycle
is different from the second half of the second cycle. In the second half of
the first cycle, we shear the system along the flow direction of the $xy$
plane, while in the second cycle, we shear the system along the flow direction
of the $-xy$ plane in the second half. This cancels the progression for the
second half when we consider two complete cycles. Thus the strain it covers in
two cycles along the flow direction of $xz$ plane is $2\dot{\gamma}T$, _i.e._
in a single cycle it covers a $\dot{\gamma}T$ strain. Therefore, the global
strain rate is also $\dot{\gamma}$, and by a similar calculation, we end up
with the expression,

$W_{\rm SW}=\frac{\eta_{\parallel}}{\eta_{\rm
ss}}\left(\frac{1}{1-\alpha}\right)+\frac{\eta_{\perp}}{\eta_{\rm
ss}}\left(\frac{1}{\alpha}\right),$ (A18)

which is Eq. (7) in the main text.

###  3 The SO protocol

In the SO, unlike all the other protocols, the suspension is sheared in a
primary direction and simultaneously a sinusoidal shear is applied
perpendicular to this direction. In this case, the primary dissipation per
strain can be written as

$W_{PRI}(per\,strain)=\eta_{\parallel}\dot{\gamma}.$ (A19)

The perpendicular dissipation can be written as

$W_{\rm OSC}(per\,strain)=\frac{\eta_{\perp}\gamma_{\rm
osc}^{2}\omega^{2}}{2\dot{\gamma}}$ (A20)

Adding both these dissipations and dividing by the steady-state dissipation
per strain we finally arrive at

$W=\frac{\eta_{\parallel}}{\eta_{ss}}+\frac{\eta_{\perp}}{\eta_{ss}}\left(\frac{\gamma_{\rm
osc}\omega}{\dot{\gamma}}\right)^{2}.$ (A21)

## References

  * Guazzelli and Pouliquen [2018] E. Guazzelli and O. Pouliquen, Rheology of dense granular suspensions, Journal of Fluid Mechanics 852, P1 (2018).
  * Heussinger  _et al._ [2010] C. Heussinger, L. Berthier, and J.-L. Barrat, Superdiffusive, heterogeneous, and collective particle motion near the fluid-solid transition in athermal disordered materials, Europhysics Letters 90, 20005 (2010).
  * Trulsson  _et al._ [2017] M. Trulsson, E. DeGiuli, and M. Wyart, Effect of friction on dense suspension flows of hard particles, Phys. Rev. E 95, 012605 (2017).
  * Olsson and Teitel [2020] P. Olsson and S. Teitel, Dynamic length scales in athermal, shear-driven jamming of frictionless disks in two dimensions, Phys. Rev. E 102, 042906 (2020).
  * Seto  _et al._ [2013] R. Seto, R. Mari, J. F. Morris, and M. M. Denn, Discontinuous shear thickening of frictional hard-sphere suspensions, Phys. Rev. Lett. 111, 218301 (2013).
  * Jamali and Brady [2019] S. Jamali and J. F. Brady, Alternative frictional model for discontinuous shear thickening of dense suspensions: Hydrodynamics, Phys. Rev. Lett. 123, 138002 (2019).
  * Dong and Trulsson [2017] J. Dong and M. Trulsson, Analog of discontinuous shear thickening flows under confining pressure, Phys. Rev. Fluids 2, 081301 (2017).
  * Singh  _et al._ [2020] A. Singh, C. Ness, R. Seto, J. J. de Pablo, and H. M. Jaeger, Shear thickening and jamming of dense suspensions: The “roll” of friction, Phys. Rev. Lett. 124, 248005 (2020).
  * Cwalina  _et al._ [2015] C. D. Cwalina, R. D. Dombrowski, C. J. McCutcheon, E. L. Christiansen, and N. J. Wagner, Mmod puncture resistance of eva suits with shear thickening fluid (stf)–armortm absorber layers, Procedia Engineering 103, 97 (2015).
  * Wagner and Wetzel [2002] N. Wagner and E. Wetzel, Advanced body armor utilizing shear thickening fluids, in  _23rd army science conference, Orlando, FL_ , Vol. 3 (2002).
  * Pednekar  _et al._ [2018] S. Pednekar, J. Chun, and J. F. Morris, Bidisperse and polydisperse suspension rheology at large solid fraction, Journal of Rheology 62, 513 (2018).
  * Nagy  _et al._ [2017] D. B. Nagy, P. Claudin, T. Börzsönyi, and E. Somfai, Rheology of dense granular flows for elongated particles, Phys. Rev. E 96, 062903 (2017).
  * Trulsson [2018] M. Trulsson, Rheology and shear jamming of frictional ellipses, Journal of Fluid Mechanics 849, 718 (2018).
  * Xu  _et al._ [2013] M. Xu, H. Liu, H. Zhao, and W. Li, How to decrease the viscosity of suspension with the second fluid and nanoparticles?, Scientific reports 3, 3137 (2013).
  * Lin  _et al._ [2016] N. Y. Lin, C. Ness, M. E. Cates, J. Sun, and I. Cohen, Tunable shear thickening in suspensions, Proceedings of the National Academy of Sciences 113, 10774 (2016).
  * Dong and Trulsson [2020] J. Dong and M. Trulsson, Transition from steady shear to oscillatory shear rheology of dense suspensions, Phys. Rev. E 102, 052605 (2020).
  * Sehgal  _et al._ [2019] P. Sehgal, M. Ramaswamy, I. Cohen, and B. J. Kirby, Using acoustic perturbations to dynamically tune shear thickening in colloidal suspensions, Phys. Rev. Lett. 123, 128001 (2019).
  * Ness  _et al._ [2018] C. Ness, R. Mari, and M. E. Cates, Shaken and stirred: Random organization reduces viscosity and dissipation in granular suspensions, Science Advances 4, eaar3296 (2018).
  * Cates  _et al._ [1998] M. E. Cates, J. P. Wittmer, J.-P. Bouchaud, and P. Claudin, Jamming, force chains, and fragile matter, Phys. Rev. Lett. 81, 1841 (1998).
  * Seto  _et al._ [2019] R. Seto, A. Singh, B. Chakraborty, M. M. Denn, and J. F. Morris, Shear jamming and fragility in dense suspensions, Granular Matter 21 (2019).
  * Blanc  _et al._ [2023] F. Blanc, F. Peters, J. J. Gillissen, M. E. Cates, S. Bosio, C. Benarroche, and R. Mari, Fragile dense suspensions under shear rotation, arXiv 2212.02165 (2023).
  * Blanc  _et al._ [2011] F. Blanc, F. Peters, and E. Lemaire, Local transient rheological behavior of concentrated suspensions, Journal of Rheology 55, 835 (2011).
  * Peters  _et al._ [2016] F. Peters, G. Ghigliotti, S. Gallier, F. Blanc, E. Lemaire, and L. Lobry, Rheology of non-brownian suspensions of rough frictional particles under shear reversal: A numerical study, Journal of Rheology 60, 715 (2016).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of computational physics 117, 1 (1995).
  * Note [1] The notation of parallel and perpendicular in the cruising protocol is a bit misleading, as none of them is along the propagating axis. Nevertheless, treating them on the same footing is useful.
  * Torquato  _et al._ [2000] S. Torquato, T. M. Truskett, and P. G. Debenedetti, Is random close packing of spheres well defined?, Phys. Rev. Lett. 84, 2064 (2000).
  * Ness  _et al._ [2022] C. Ness, R. Seto, and R. Mari, The physics of dense suspensions, Annual Review of Condensed Matter Physics 13, 97 (2022).
  * Note [2] We used two different $\phi_{c}$ values for our data and the data of Ref. [18]. Both $\phi_{c}$ values were extracted by fitting the steady state curves in with $\eta=\eta_{0}(\phi-\phi_{C})^{-2}$. Our data yielded $\phi_{c}=0.5900$ while the data from Ref. [18] gave $\phi_{c}=0.5827$.
  * Pine  _et al._ [2005] D. J. Pine, J. P. Gollub, J. F. Brady, and A. M. Leshansky, Chaos and threshold for irreversibility in sheared suspensions, Nature 438, 997 (2005).
  * Corté  _et al._ [2008] L. Corté, P. M. Chaikin, J. P. Gollub, and D. J. Pine, Random organization in periodically driven systems, Nature Physics 4, 420 (2008).
  * Weijs  _et al._ [2015] J. H. Weijs, R. Jeanneret, R. Dreyfus, and D. Bartolo, Emergent hyperuniformity in periodically driven emulsions, Phys. Rev. Lett. 115, 108301 (2015).


# Numerical investigation of a segmented-blade ion trap with biasing rods

Jungsoo Hong  Myunghun Kim  Hyegoo Lee  Moonjoo Lee<EMAIL_ADDRESS>Department of Electrical Engineering, Pohang University of Science and
Technology (POSTECH), 37673 Pohang, Korea

###### Abstract

We report a numerical study of a linear ion trap that has segmented blades and
biasing rods. Our system consists of radio frequency (rf) blades, dc blades
with ten separate electrodes, and two biasing rods for compensating the ions’
micromotion. After calculating the optical access for the ions, we find rf and
dc voltages that result in a stable trapping configuration of 171Yb+ ions. We
also explore the micromotion compensation with the biasing rods, and calculate
the influence of blade misalignment to the trap potential. Our work offers
quantitative understanding of the trap architecture, assisting reliable
operation of an ion-trap quantum computer.

##  I Introduction

Trapped atomic ions are one of the best-controlled systems for quantum
information processing [1, 2, 3]. The ions are typically stored from days to
weeks in an ultrahigh vacuum environment. The record infidelity is as low as
10-6 for single-qubit gates [4] and 10-3 for two-qubit gates [5, 6], while the
state preparation and measurement errors are below 10-4 [4, 7]. The ion
qubit’s coherence time is longer than one hour in optimized conditions [8, 9].
Quantum entanglement of 20 ions has been achieved through addressing the
individual ions [10], and even entanglement of 24 ions was demonstrated in a
very compact setting [11]. Moreover, the “full connectivity”, mediated by
collective motion of the ions, makes it possible to execute quantum algorithms
in efficient ways [12].

Among a number of ion-trap architectures [13, 14], a macroscopic trap with
thin segmented blades [15] would be an excellent setting for quantum
experiments. This trap contains two diagonal rf blades in order to generate a
psuedopotential for confining the ions along the transverse directions. The
axial potential is formed by two dc blades aligned along the anti-diagonal
direction, via applying dc voltages to each of the segmented electrodes. Owing
to the comparatively smaller thickness of the blades than those of e.g., Refs.
[16, 17], the interblade distance could become short, which enables the
generation of deep trapping potentials. Seminal experiments with this trap
were reported, including the realization of a modular quantum network [18],
small programmable quantum computer [19], and global or parallel entangling
gates [20, 21].

We point out three major advantages of the segmented-blade trap. First, the
trap offers a great optical access for the qubit control and measurement [18,
22]. It is possible to address the ions from radial directions with a
numerical aperture (NA) as large as 0.6 [18]. Second, the trap provides the
capability of engineering various potential shapes along the trap axis.
Adjusting dc voltages to each segment, the trap can create harmonic, quartic,
and multi-well potentials [23]. Lastly, this trap exhibited great performance
at cryogenic temperatures: Recent experiment demonstrated the trapping of more
than 100 ions in a cryogenic environment [24], representing the scalability of
qubit numbers in this platform.

Here, we present a comprehensive numerical study on this segmented-blade ion
trap. In addition to the structure of Ref. [15], we include two biasing rods
[25, 26, 27, 28, 29] solely for compensating the micromotion. We investigate
the optical access, and calculate the potential energy for trapping 171Yb+
ions. Our calculation shows that the micromotion can be compensated by
applying dc voltages to the biasing rods and dc electrodes. We also describe
the influence of blade misalignment to the trap potential.

Figure 1:  (a) Overview of segmented blade trap. Length of the blade $L$ is
25.2 mm, and diameter of the biasing rod (gray) is 1 mm. (b) Top view of dc
blade. Zone in gold denotes metal-coated facet; white rectangles indicate
insulating ceramic surface. Segments from A to E are disconnected
electrically. Holes are used to insert screws for fixing the blade to mount
(not shown). Width of tapered end and groove, marked as $t$, is 50 $\mu$m. (c)
Central area of the trap in the $xy$ plane. Thickness of the blade, referred
to as $T$, is 300 $\mu$m; tapered length $l$ is 1 mm; angle between the blades
$\theta$ is 120∘. Blade-blade distance is given by $d$.  Figure 2:  (a)
Configuration when blade-blade distance $d<d_{0}$. Red circle corresponds to
point 1 and blue to point 2. Available NA is determined by
$\theta_{1}(<\theta_{2})$. (inset) Zoomed-in view of central area. (b) Actual
NA is limited by $\theta_{2}(<\theta_{1})$ at $d>d_{0}$. The angle
$\theta_{1}$ coincides with $\theta_{2}$ at $d=d_{0}$. (c) NAs are plotted as
a function of $d$, along $x$ direction in black and $y$ direction in red.
NAx/y,1/2/r means the NA limited along $x/y$ axis by the point $1,2$, and rod.
Actual NAs are indicated by solid lines. NAx,1/2 works for $x>0$, and NAx,r
for $x<0$. Inset shows $\theta_{\rm{r}}$.

##  II Result

###  II.1 Trap design

As presented in Fig. 1(a), we consider four segmented blades (gold) in a
linear trap configuration, with two biasing rods (gray). Two diagonal blades
are driven with a rf voltage, and the anti-diagonal ones with ten different dc
voltages. While all five segments in the rf blade are electrically connected,
the dc blades’ segments are disconnected by the grooves and nonmetal zones
(white area in Figs. 1(a) and (b)). We can apply different dc voltages to the
zones from A to E. The trap potential along the $xy$ directions is mostly
derived from the rf voltage with minor contribution from the dc segments; the
potential along the $z$ direction originates from the voltages to all the dc
electrodes. Fig. 1(c) shows a diagonal blade-blade distance of $d$ at an angle
$\theta=120^{\circ}$ between the rf and dc blades.

Figure 3:  Contour plot of calculated pseudopotential on $xy$ plane at $z=0$
(center of trap) in (a), and $xz$ plane at $y=0$ in (b), with
$V_{\rm{rf}}=600$ V. Scale bar indicates 1 mm. Slicecut of the trap potential
is shown with black lines along $x$ axis in (c), and $y$ axis in (d). Central
regions are focused in (e) and (f). Black solid line is result of the
numerical simulation; red dashed line is fit to a quadratic function; vertical
dashed lines indicate fitting range. (g), (h) Trap frequency and depth along
$x$ and $y$ directions as a function of blade-blade distance $d$.

###  II.2 Optical access

We discuss the optical access for measuring the ion fluorescence and
addressing individual ions. Figs. 2(a) and (b) show that the two tapering
corners 1 and 2 determine available NAs. When the blade-blade distance $d$ is
smaller than $d_{0}=2l\times t/(T-t)=400$ $\mu$m, we find that $\theta_{1}$ is
smaller than $\theta_{2}$, resulting in a situation where the NA is limited by
$\theta_{1}$. The size of the two angles is reversed at $d>d_{0}$ so that the
NA is determined by $\theta_{2}$ in this regime; the two angles become
identical for the case of $d=d_{0}$. Considering the trap geometry, we derive
analytical formulas for the NAs along the $x$ and $y$ directions. When $d<400$
$\mu$m, we obtain

$\displaystyle\textrm{NA}_{x,1}$
$\displaystyle=\sin{\left(\frac{\pi}{3}-\arctan{\left(\frac{t}{d}\right)}\right)}$
$\displaystyle\textrm{NA}_{y,1}$
$\displaystyle=\sin{\left(\frac{\pi}{6}-\arctan{\left(\frac{t}{d}\right)}\right)},$

and the expressions below describe the NAs at $d>400$ $\mu$m,

$\displaystyle\textrm{NA}_{x,2}$
$\displaystyle=\sin{\left(\frac{\pi}{3}-\arctan{\left(\frac{T}{2l+d}\right)}\right)}$
$\displaystyle\textrm{NA}_{y,2}$
$\displaystyle=\sin{\left(\frac{\pi}{6}-\arctan{\left(\frac{T}{2l+d}\right)}\right)},$

with the thickness of the tapered end $t$, tapered length $l$, and untapered
thickness of the blade $T$. The NA determined by the point $i$ along the $x/y$
axis is denoted by NAx/y,i $(i=1,2)$ plotted in Fig. 2(c). The actual NA
corresponds to the smaller value between the two NAs at a given distance. For
example, we attain NA${}_{x}=0.77$ and NA${}_{y}=0.35$ at $d=300$ $\mu$m that
is limited by the point 1, while at $d=500$ $\mu$m, the NAx increases to 0.80
and NAy to 0.39 determined by the point 2. It is notable that, for the
direction $x<0$, the optical access is limited by the biasing rods
$\theta_{\rm{r}}$ with NA$=0.46$.

###  II.3 Simulation method

The total trap potential $\phi(\mathbb{r})_{\rm{tot}}$ is given by the sum of
the pseudopotential generated by a rf voltage applied to the two rf blades,
static potential $\phi(\mathbb{r})_{\rm{dc}}$ derived from dc voltages on each
segment of the dc blades, and another static potential
$\phi(\mathbb{r})_{\rm{rod}}$ from biasing rods:

$\displaystyle\phi(\mathbb{r})_{\rm{tot}}$
$\displaystyle=\phi(\mathbb{r})_{\rm{rf}}+\phi(\mathbb{r})_{\rm{dc}}+\phi(\mathbb{r})_{\rm{rod}}$
$\displaystyle=\frac{q^{2}}{4m\Omega^{2}}|\vec{E}(\mathbb{r})|^{2}+q\sum_{i}V(\mathbb{r})_{i}+qV(\textbf{r})_{\rm{rod}},$

where $q$ and $m$ represent the charge and mass of a single ion, respectively
and $\vec{E}(\mathbb{r})$ refers to the electric field stemmed from the rf
voltage, $V(\mathbb{r})_{i}$ is the voltage derived from the dc electrode $i$,
and $V(\textbf{r})_{\rm{rod}}$ is the dc voltage from biasing rods. We choose
the driving frequency $\Omega=2\pi\times 22.5$ MHz in this paper. All
calculations are done for 171Yb+ ions using the software ANSYS Electronics
Desktop 2020 R1.

###  II.4 Pseudopotential

Our simulation results of the rf pseudopotential are presented in Fig. 3. We
begin with the case of $d=500$ $\mu$m. Both rf blades are driven with an
amplitude of $V_{\rm{rf}}=600$ V, and all dc segments and biasing rods are
ground. In this configuration, as shown in Figs. 3(a) and (b), a
pseudopotential “tube” is formed along the $z$ axis. We then obtain the trap
frequency $\omega_{r}/2\pi$ via fitting the potential to a quadratic function
$\phi(r)=m\omega_{r}^{2}r^{2}/2$, for $r=x$ and $y$. The fitting range is
limited to the central region (from $-20$ to $+20~{}\mu$m with respect to the
potential minimum, marked in Figs. 3(c)-(f)) where the contribution of higher
order terms is negligible. The yielded trap frequencies are 3.2 MHz along both
the $x$ and $y$ axes. Concerning the trap depth $U_{0}$, we define the depth
as the maximum value of the potential energy in a given direction: We obtain
6.0 eV along the $x$ axis and 36 eV to the $y$ direction.

Figure 4:  Static potential $\phi(\textbf{r})_{\rm{dc}}$ derived from dc
voltages of $V_{\rm{A}}=V_{\rm{E}}=15$ V, $V_{\rm{B}}=V_{\rm{D}}=3$ V,
$V_{\rm{C}}=-10$ V, along $x$ axis in (a) and $y$ axis in (b) at $z=0$. Black
solid line is the simulation result; red dashed line is fit to a quadratic
function; the fitting range is defined by vertical dashed lines.  Figure 5:
(a) Contour plot of total potential $\phi(x,y,0)_{\rm{tot}}$ derived from rf
and dc voltages given in main text. Principal axis is rotated by $\zeta$ with
respect to $x$ axis. Scale bar indicates 100 $\mu$m. Equipotential line along
1 eV is shown, and several lines in inset. Color scale is from 0 to 40 eV in
(a) and 0 to 0.1 eV in inset. Total potential along $x$ axis in (b), $y$ axis
in (c), and $z$ axis in (d). Black solid line is simulation result; red dashed
line is fit to quadratic function; fitting range is defined by vertical dashed
lines. (e) Contour plot of total trap potential $\phi(x,y,\pm
1~{}\textrm{mm})_{\rm{tot}}$. Principal axis is rotated by
$180^{\circ}-\zeta^{\prime}$ with respect to $x$ axis.

We remark two features of this simulation result. First, the calculation shows
that the ion confinement is stronger along the $y$ axis. It is because the
distance between the rf and dc blades is smaller than that of the $x$
direction, leading to a higher electric field magnitude along the $y$ axis.
Second, as indicated in Figs. 3(e) and (f), it is notable that the minimum
point of the pseudopotential is shifted from the origin by (0.0, 1.4) $\mu$m.
It is due to the fact that the potential energy, created between the biasing
rod and rf blade at $x<0$ and $y>0$, shifts the rf-null point to $y>0$
slightly (see Appendix). Near the origin, this potential is approximated as
$\phi(0,y)_{\rm{rf}}\simeq-c_{y}\cdot y$ ($c_{y}$ is positive). The total
pseudopotential is then given by
$\phi(0,y)_{\rm{rf}}=\frac{1}{2}m\omega_{y}^{2}(y-\frac{c_{y}}{m\omega_{y}^{2}})^{2}+const.$
along the $y$ axis, explaining the shifting direction of the potential
minimum.

We continue such calculations for several $d$ (Figs. 3(g) and (h)) at
$V_{\rm{rf}}=600$ V. Both the trap frequency and depth decrease as $d$
increases. The trap depth is higher than several eV for all the distances, and
the trap frequency decreases from 5.2 MHz to 0.7 MHz along both the $x$ and
$y$ axes, as $d$ changes from 400 $\mu$m to 1 mm.

###  II.5 Confinement in three dimensions

The axial confinement of the ions is made by the static potential
$\phi(\mathbb{r})_{\rm{dc}}$. Here, we concentrate on the generation of a
single-well, harmonic potential along the $z$ direction. A positive dc voltage
set of $V_{\rm{A}}=V_{\rm{E}}$ and $V_{\rm{B}}=V_{\rm{D}}$ gives rise to a
stable symmetric potential along the $z$ axis. However, due to Earnshaw’s
theorem [30], the dc voltage set brings about a deconfinement effect along the
radial directions (Fig. 4) — the depth of pseudopotential is typically much
larger than that of the deconfinement potential, which ensures the trapping of
ions. Overall, the sum of the static potential $\phi(\mathbb{r})_{\rm{dc}}$
and pseudopotential $\phi(\mathbb{r})_{\rm{rf}}$ provides a stable total
trapping potential $\phi(\mathbb{r})_{\rm{tot}}$, which can be fitted with a
quadratic function along all three directions.

In Fig. 5, the total trap potential $\phi(\mathbf{r})_{\rm{tot}}$ is shown in
our chosen voltage set of $V_{\rm{rf}}=600$ V, $V_{\rm{A,E}}=15$ V,
$V_{\rm{B,D}}=3$ V, and $V_{\rm{C}}=-10$ V at $d=500$ $\mu$m. The biasing rods
are ground. The potential is well fitted with a quadratic function along the
three axes, and the fittings give the trap frequencies of
$(\omega_{x},\omega_{y},\omega_{z})/2\pi=(3.2,3.1,1.1)$ MHz for 171Yb+ ion,
with the depths of 6.1, 34, 9.1 eV at each direction, constituting a stable
trapping configuration. The potential minimum is located at $(0.2,1.6,0.0)$
$\mu$m. Despite the deconfinement effect of the dc voltages in Fig. 4, the
pseudopotential dominates the trapping potential along the $x$ and $y$
directions. Note that the impression of a negative voltage to $V_{\rm{C}}$
offers both deeper axial potential and higher trap frequency. When
$V_{\rm{C}}=0$ V, the axial trap frequency $\omega_{z}/2\pi$ decreases to
$0.63$ MHz and the trap depth along the $z$ axis to 6.5 eV.

It is noteworthy to mention that the principal axis is determined by the
applied dc voltages. In Fig. 5(a), we plot a slicecut of
$\phi(\textbf{r})_{\rm{tot}}$ at $z=0$. In this region, the effect of dc
voltage $V_{\rm{C}}=-10$ V governs the symmetry of the potential: The negative
dc voltages reduce confinement along the both dc electrodes, which causes that
the long principal axis aligns along the dc blades. The angle between this
axis and $x$ axis is $\zeta=35^{\circ}$ obtained by fitting the equipotential
line with an ellipse. This phenomenon is reversed, as presented in Fig. 5(e),
at $z=\pm 1~{}$mm where the effect of $V_{\rm{A/E}}=15$ V is maximized. The
positive dc voltage increases the confinement along the dc blades, and thus
the long principal axis of the equipotential line aligns along the rf blade.
The principal axis rotates from $\zeta$ to
$180^{\circ}-\zeta^{\prime}=139^{\circ}$ continuously as $z$ coordinate
changes from $z=0$ to $z=\pm 1$ mm.

Figure 6:  (a) Position of potential minimum as a function of $V_{\rm{rod}}$
with $\Delta V_{\rm{dc}}=0$. Solid lines are fit to linear functions. (b)
Potential minimum point as $\Delta V_{\rm{dc}}$ changes at fixed
$V_{\rm{rod}}=0$. (c) Position of potential minimum as $\Delta V_{\rm{dc}}$
changes like (b), however the position shift along the $x$ direction is
canceled by varying $V_{\rm{rod}}$. Position error is dominated by spatial
resolution of simulation, 100 nm (error bars not shown, similar with symbol
size).

###  II.6 Micromotion compensation

Next, we describe the simulation of micromotion compensation in our system. As
discussed in Refs. [31, 32], the micromotion could be compensated by adding
independent and different dc offset voltages to each rf blade. We approach
this from a different perspective: Via including the biasing rods (gray in
Fig. 1(a)), we have two independent control parameters, those rods and dc
electrodes, so as to move the position of potential minimum arbitrarily in the
transverse plane.

We first quantify the impact of biasing rods to the location of the potential
minimum. Applying $V_{\rm{rf}}=600$ V to rf blades with dc blades ground, we
change the voltage on both the biasing rods $V_{\rm{rod}}$ from $-100$ to
$100$ V, and search for the potential minimum of
$\phi(\textbf{r})_{\rm{tot}}$. The positions of potential minimum are plotted
in Fig. 6(a): The location is only shifted along the $x$ axis over $\sim$10
$\mu$m. The linear fitting in Fig. 6(a) gives a displacement per
$V_{\rm{rod}}$ of $45.3(4)$ nm$/V_{\rm{rod}}$ along the $x$ direction. In
terms of electric field, given $V_{\rm{rod}}=1$ V to both rods, an electric
field of (33.3(1), 0.0) V/m is generated at the origin, and at the rf-null
position of $(0.0,1.4)$ $\mu$m, slight electric field along the $y$ direction
is obtained as $(33.3(1),0.3(1))$ V/m. The shift of the potential minimum
induced by the small $y$ component is negligible along the $y$ axis.

Figure 7:  (a) Rf blade at $x<0,y>0$ is tilted by 1.5∘ with respect to the $z$
axis. Semi-transparent planes correspond to the cross sections shown in
(b)-(f). Contour plot of pseudopotential at $z=-10$ mm in (b), $z=-5$ mm in
(c), $z=0$ mm in (d), $z=5$ mm in (e), and $z=10$ mm in (f). Black solid line
is the slicecut along the $x$ and $y$ axis, fitting range is defined by black
dashed line, red dashed line is quadratic fitting. (g) Coordinates of rf null
are plotted along the $z$ axis.  Figure 8:  (a), (b) Cross sections of the
pseudopotential in the $yz$ plane (black lines). Contour plot of the
pseudopotential at $x=-50$ $\mu$m in (c), $x=-25$ $\mu$m in (d), $x=0$ $\mu$m
in (e), $x=25$ $\mu$m in (f), and $x=50$ $\mu$m in (g). White dot indicates
the minimum of pseudopotential. Black solid line is the slicecut along the $z$
axis at $y=0$ of the given plane, black dashed line defines the fitting range,
and red dashed line is quadratic fitting.  Figure 9:  Contour plot of
pseudopotential in the $yz$ plane at $x=0.2$ $\mu$m and slicecut at $y=1.6$
$\mu$m in (a), dc potential in (b), and total potential in (c). Scale bar
denotes $100$ $\mu$m. Black solid line is the slicecut, red dashed line is
quadratic fit, and black dashed line shows the fitting range.

We proceed to investigate the effect of the voltage change in dc electrodes to
the shift of the potential minimum. We simultaneously add a common offset
voltage $\Delta V_{\rm{dc}}$ to all the dc electrodes at $x>0,y>0$ and
$-\Delta V_{\rm{dc}}$ to every dc electrode in $x<0,y<0$. The results are
plotted in Fig. 6(b), showing that the voltage impression of $\Delta
V_{\rm{dc}}$ moves the potential minimum along the direction of
$\pm(\hat{x}+\hat{y})/\sqrt{2}$, in the central region of the trap. The
positions are moved over $\sim$10 $\mu$m as $\Delta V_{\rm{dc}}$ changes by 4
V. From linear fittings in Fig. 6(b), we obtain the slopes of 2.64(1)
$\mu$m$/\Delta V_{\rm{dc}}$ and 2.66(1) $\mu$m$/\Delta V_{\rm{dc}}$ along the
$x$ and $y$ axis, respectively.

Armed with such understanding of the biasing rods and dc electrodes, we seek
voltage sets that move the potential minimum solely along the $y$ axis. To
this end, with trials and errors, we obtain a relation between $V_{\rm{rod}}$
and $\Delta V_{\rm{dc}}$ as $V_{\rm{rod}}=-59.1\times\Delta V_{\rm{dc}}$,
which cancels the $x$-displacement by $\Delta V_{\rm{dc}}$ with
$V_{\rm{rod}}$. We remark that this value is slightly larger than the scaling
of the above values as $V_{\rm{rod}}=-58.3(4)\times\Delta V_{\rm{dc}}$. It is
because the $x$ component of an electric field by $V_{\rm{rod}}$ gently
decreases as $|y|$ $\mu$m increases. Note that Fig. 6(a), from which we find
the above relation of 45.3(4) nm/$V_{\rm{rod}}$, is obtained at $y=1.4$ $\mu$m
where the $y$ component of the field is close to zero (and the $x$ component
of an electric field per $V_{\rm{rod}}$ is almost maximum). However, we look
for $V_{\rm{rod}}$ that removes the $x$ displacement at different $y$
coordinates, where the small $y$ component of the field appears and $x$
component decreases accordingly: as $|y|$ $\mu$m becomes larger, given the
same $V_{\rm{rod}}$, the displacement along the $x$ axis (by $V_{\rm{rod}}$)
decreases, which requires higher $V_{\rm{rod}}$ to cancel out the $x$
displacement by identical $\Delta V_{\rm{dc}}$. This small discrepancy was
solved by our “trials and errors”—we simulate many times to find the voltage
ratio of $-59.1$ that suppresses the net $x$ displacement as much as possible.
Fig. 6(c) presents simulation results that the $y$ position is moved over
$\sim$10 $\mu$m when $\Delta V_{\rm{dc}}$ varies by $4$ V; the displacement
along the $x$ axis is canceled by impressing $V_{\rm{rod}}$ simultaneously.
The linear fitting of Fig. 6(c) yields the $y$-displacement per $\Delta
V_{\rm{dc}}$ of 2.65(1) $\mu$m/$\Delta V_{\rm{dc}}$ along the $y$ axis, which
is very close to the value of Fig. 6(b). Combining the results along both the
axes, it would be possible to locate the potential minimum at arbitrary points
in the transverse plane. We note that all error bars in this subsection are
originated from the precision of our numerical simulation.

###  II.7 Blade misalignment

####  II.7.1 Potential minimum

Our last investigation focuses on the influence of blade misalignment to the
trapping potential. Although all considerations above are presumed on perfect
alignment of the blades, in practice, it is not possible to avoid angular or
lateral misalignments of the blades. Such misalignment distorts the shape of
trap potentials, which would cause an instability of the ions. Here we
describe one of such misalignments, tilt between the blades. Fig. 7(a)
presents a situation in which one rf blade ($x<0,y>0$) is tilted at an angle
of $1.5^{\circ}$ with respect to the $z$ axis. The simulation results are
shown in Figs. 7(b)–(f) when $V_{\rm{rf}}=600$ V, and other dc electrodes and
biasing rods are ground. In case of $z<0$, the rf blade is located more
distantly from the origin with respect to other blades. Therefore, the effect
of rf voltages in this blade is reduced, moving the rf null position toward
the region of $x<0,y>0$ (Figs. 7(b), (c)). By contrast, for the regime of
$z>0$, this tendency appears in an opposite way so that the rf null is located
at $x>0,y<0$ (Figs. 7(e), (f)). Fig. 7(g) shows the calculation results of rf
nulls along the trap axis. We also present both the cross section in the $yz$
plane and slicecut along the $z$ direction as the $x$ coordinate changes
(Figs. 8(c)–(g)). The potential is very asymmetric along the $z$ axis, and the
trap depth is overall larger in the region of $z>0$ than that of $z<0$.

####  II.7.2 Axial micromotion

We point out that this kind of blade rotation causes an axial micromotion [33,
34], which would be very difficult to be compensated. In our consideration,
whereas the axis of $\phi(\textbf{r})_{\rm{dc}}$ is in parallel with the $z$
axis, the axis of the rf tube is not: In consequence, the axis of
$\phi(\textbf{r})_{\rm{tot}}$, along which the ions aligned, is not parallel
with the axis of $\phi(\textbf{r})_{\rm{rf}}$ — this induces an inevitable
axial micromotion. Therefore, it is crucial to align the blades as parallel as
possible over the course of actual trap development.

We quantify the axial micomotion in this tilted configuration. With the same
rf voltage and dc voltage set given in Sec. II.5, our calculation of
$\phi(\textbf{r})_{\rm{tot}}$ finds the potential minimum at
$\mathbf{\bar{r}}=(0.2,1.6,0.5)$ $\mu$m (Fig. 9), at which the ion is located.
We note that, under perfect alignment, the potential minimum is found at (0.2,
1.6, 0.0) $\mu$m in Sec. II.5: this $z$-position shift is induced by the tilt
misalignment, and the rf field at $\mathbf{\bar{r}}$ along the $z$ direction
causes the axial micromotion. Our numerical simulation provides an amplitude
of the rf field $E_{\rm{rf,z}}(\mathbf{\bar{r}})=59.8$ V/m along the $z$ axis
with an equivalent energy of $\phi(\mathbf{\bar{r}})_{\rm{rf}}=25.2$ neV. This
energy corresponds to an “excess” temperature of $T_{\rm{axial}}=585$ $\mu$K.

##  III Discussion

We remark four points regarding the present study. First, we describe how we
choose the location of the rods. The coordinates of the two rods’ centers are
$(x,y)=(-2.37,\pm 1.81)$ mm, determined by two practical reasons. One point is
that, in a given position, it should be possible to move the rf-null over tens
of microns in hundreds of volts, in which we can apply using a commercially
available power supply. Another reason consists of the optical access: Along
this direction, our numerical aperture (NA) is 0.40 limited by a recessed
viewport of our vacuum chamber. Therefore we aim to attain a NA governed by
the rods larger than 0.40, which does not influence imaging of the ions.

Next, during the investigation of micromotion compensation in Sec. II.6, one
might think that the $y$-position shift would be possible by installing two
biasing rods at the area $y>0$. We separately calculate the shift of the
potential minimum in such configuration, however, the displacement is very
small that the position only moves by $\sim 0.2$ $\mu$m in 1000 V. It is
because, the electric fields from biasing rods at $y>0$ are more shielded by
the blades than the fields from the rods at $x<0$, due to the asymmetric
angles between the blades.

Third, from the perspective of quantum networks and distributed quantum
computing, optical interface with the segmented-blade trap was mostly done
with high-NA lenses so far [18, 35, 36]. Alternatively, we point out that it
would also be interesting to combine the trap with optical cavities [37, 38,
39]. It would be informative if future theoretical studies could include a
strategy for coupling fiber cavities [40, 41, 42, 43], the influence of
dielectric surface charges to the trap potential [44], and the impact of high
voltage on piezoelectric transducer (for scanning the cavity length) to the
potential energy.

Lastly, while we herein concentrate on the calculation of the trap potential,
the vacuum pressure, particularly the local pressure near the ions, is also an
important factor for stable ion trapping. It is notable that the software
MolFlow+ [45] could be used for estimating the local pressure at the location
of the ions. By doing so, we expect that the trap mount and other vacuum
components could be designed and installed in an optimized geometry: the
system should be designed such that the conductance from the ion position to
the vacuum pumps is maximized.

##  IV Conclusion

In conclusion, we have presented a numerical investigation of a linear ion
trap, including four segmented blades and two biasing rods. We find that it is
possible to achieve an excellent optical access for the radial directions.
Also, the trap depth and frequency are calculated in a given voltage set,
which offers a stable trapping configuration of 171Yb+ ions. We then describe
an approach to compensate the micromotion, followed by a calculation regarding
the blade misalignment. Our study provides thorough understanding of this ion-
trapping system, forming the basis of stable ion-based quantum computer and
simulator.

##  V acknowledgments

We thank Kihwan Kim for helpful discussions. This work has been supported by
National Research Foundation (Grant No. 2019R1A5A1027055), Institute for
Information and Communication Technology Planning and Evaluation (IITP, Grant
No. 2022-0-01040), Samsung Science and Technology Foundation (SSTF-BA2101-07,
SRFC-TC2103-01), Samsung Electronics Co., Ltd. (IO201211-08121-01), and BK21
FOUR program.

## References

  * Leibfried  _et al._ [2003] D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Rev. Mod. Phys. 75, 281 (2003).
  * Blatt and Wineland [2008] R. Blatt and D. Wineland, Nature 453, 1008 (2008).
  * Monroe and Kim [2013] C. Monroe and J. Kim, Science 339, 1164 (2013).
  * Harty  _et al._ [2014] T. P. Harty, D. T. C. Allcock, C. J. Ballance, L. Guidoni, H. A. Janacek, N. M. Linke, D. N. Stacey, and D. M. Lucas, Phys. Rev. Lett. 113, 220501 (2014).
  * Ballance  _et al._ [2016] C. J. Ballance, T. P. Harty, N. M. Linke, M. A. Sepiol, and D. M. Lucas, Phys. Rev. Lett. 117, 060504 (2016).
  * Gaebler  _et al._ [2016] J. P. Gaebler, T. R. Tan, Y. Lin, Y. Wan, R. Bowler, A. C. Keith, S. Glancy, K. Coakley, E. Knill, D. Leibfried, and D. J. Wineland, Phys. Rev. Lett. 117, 060505 (2016).
  * Crain  _et al._ [2019] S. Crain, C. Cahall, G. Vrijsen, E. E. Wollman, M. D. Shaw, V. B. Verma, S. W. Nam, and J. Kim, Commun. Phys. 2, 97 (2019).
  * Wang  _et al._ [2017] Y. Wang, M. Um, J. Zhang, S. An, M. Lyu, J.-N. Zhang, L.-M. Duan, D. Yum, and K. Kim, Nat. Photonics 11, 646 (2017).
  * Wang  _et al._ [2021] P. Wang, C.-Y. Luan, M. Qiao, M. Um, J. Zhang, Y. Wang, X. Yuan, M. Gu, J. Zhang, and K. Kim, Nat. Commun. 12, 233 (2021).
  * Friis  _et al._ [2018] N. Friis, O. Marty, C. Maier, C. Hempel, M. Holzäpfel, P. Jurcevic, M. B. Plenio, M. Huber, C. Roos, R. Blatt, and B. Lanyon, Phys. Rev. X 8, 021012 (2018).
  * Pogorelov  _et al._ [2021] I. Pogorelov, T. Feldker, C. D. Marciniak, L. Postler, G. Jacob, O. Krieglsteiner, V. Podlesnic, M. Meth, V. Negnevitsky, M. Stadler, B. Höfer, C. Wächter, K. Lakhmanskiy, R. Blatt, P. Schindler, and T. Monz, PRX Quantum 2, 020343 (2021).
  * Wright  _et al._ [2019] K. Wright, K. M. Beck, S. Debnath, J. M. Amini, Y. Nam, N. Grzesiak, J.-S. Chen, N. C. Pisenti, M. Chmielewski, C. Collins, K. M. Hudek, J. Mizrahi, J. D. Wong-Campos, S. Allen, J. Apisdorf, P. Solomon, M. Williams, A. M. Ducore, A. Blinov, S. M. Kreikemeier, V. Chaplin, M. Keesan, C. Monroe, and J. Kim, Nat. Commun. 10, 5464 (2019).
  * Siverns and Quraishi [2017] J. D. Siverns and Q. Quraishi, Quantum Inf. Process. 16, 314 (2017).
  * Romaszko  _et al._ [2020] Z. D. Romaszko, S. Hong, M. Siegele, R. K. Puddy, F. R. Lebrun-Gallagher, S. Weidt, and W. K. Hensinger, Nat. Rev. Phys. 2, 285 (2020).
  * Hucul [2015] D. A. Hucul, Ph.D. thesis, University of Maryland (2015).
  * Huber  _et al._ [2008] G. Huber, T. Deuschle, W. Schnitzler, R. Reichle, K. Singer, and F. Schmidt-Kaler, New J. Phys. 10, 013004 (2008).
  * Schnitzler  _et al._ [2009] W. Schnitzler, N. M. Linke, R. Fickler, J. Meijer, F. Schmidt-Kaler, and K. Singer, Phys. Rev. Lett. 102, 070501 (2009).
  * Hucul  _et al._ [2015] D. Hucul, I. V. Inlek, G. Vittorini, C. Crocker, S. Debnath, S. M. Clark, and C. Monroe, Nat. Phys. 11, 37 (2015).
  * Debnath  _et al._ [2016] S. Debnath, N. M. Linke, C. Figgatt, K. A. Landsman, K. Wright, and C. Monroe, Nature 536, 63 (2016).
  * Lu  _et al._ [2019] Y. Lu, S. Zhang, K. Zhang, W. Chen, Y. Shen, J. Zhang, J.-N. Zhang, and K. Kim, Nature 572, 363 (2019).
  * Figgatt  _et al._ [2019] C. Figgatt, A. Ostrander, N. M. Linke, K. A. Landsman, D. Zhu, D. Maslov, and C. Monroe, Nature 572, 368 (2019).
  * Siverns  _et al._ [2017] J. D. Siverns, X. Li, and Q. Quraishi, Appl. Opt. 56, B222 (2017).
  * Lin  _et al._ [2009] G.-D. Lin, S.-L. Zhu, R. Islam, K. Kim, M.-S. Chang, S. Korenblit, C. Monroe, and L.-M. Duan, Europhys. Lett. 86, 60004 (2009).
  * Pagano  _et al._ [2018] G. Pagano, P. W. Hess, H. B. Kaplan, W. L. Tan, P. Richerme, P. Becker, A. Kyprianidis, J. Zhang, E. Birckelbaw, M. R. Hernandez, Y. Wu, and C. Monroe, Quantum Sci. Technol. 4, 014004 (2018).
  * Berkeland  _et al._ [1998] D. Berkeland, J. Miller, J. Bergquist, W. Itano, and D. Wineland, J. Appl. Phys. 83, 5025 (1998).
  * McLoughlin  _et al._ [2011] J. J. McLoughlin, A. H. Nizamani, J. D. Siverns, R. C. Sterling, M. D. Hughes, B. Lekitsch, B. Stein, S. Weidt, and W. K. Hensinger, Phys. Rev. A 83, 013406 (2011).
  * Schindler  _et al._ [2013] P. Schindler, D. Nigg, T. Monz, J. T. Barreiro, E. Martinez, S. X. Wang, S. Quint, M. F. Brandl, V. Nebendahl, C. F. Roos, M. Chwalla, M. Hennrich, and R. Blatt, New J. Phys. 15, 123012 (2013).
  * Xie  _et al._ [2021] Y. Xie, J. Cui, M. D’Onofrio, A. J. Rasmusson, S. W. Howell, and P. Richerme, Quantum Sci. Technol. 6, 044009 (2021).
  * Kim  _et al._ [2022] M. Kim, K. Kim, J. Hong, H. Lee, Y. Moon, W. C. Lee, S. Kim, T. Ha, J.-Y. Sim, and M. Lee, AIP Adv. 12, 115006 (2022).
  * Steane [1997] A. Steane, Appl. Phys. B 64, 623 (1997).
  * Debnath [2016] S. Debnath, Ph.D. thesis, University of Maryland (2016).
  * An [2017] S. An, Ph.D. thesis, Tsinghua University (2017).
  * Herschbach  _et al._ [2012] N. Herschbach, K. Pyka, J. Keller, and T. E. Mehlstäubler, Appl. Phys. B 107, 891 (2012).
  * Low [2019] P. J. Low, Master’s thesis, University of Waterloo (2019).
  * Craddock  _et al._ [2019] A. N. Craddock, J. Hannegan, D. P. Ornelas-Huerta, J. D. Siverns, A. J. Hachtel, E. A. Goldschmidt, J. V. Porto, Q. Quraishi, and S. L. Rolston, Phys. Rev. Lett. 123, 213601 (2019).
  * Siverns  _et al._ [2019] J. D. Siverns, J. Hannegan, and Q. Quraishi, Sci. Adv. 5, eaav4651 (2019).
  * Lee  _et al._ [2019a] M. Lee, K. Friebe, D. A. Fioretto, K. Schüppert, F. R. Ong, D. Plankensteiner, V. Torggler, H. Ritsch, R. Blatt, and T. E. Northup, Phys. Rev. Lett. 122, 153603 (2019a).
  * Walker  _et al._ [2020] T. Walker, S. V. Kashanian, T. Ward, and M. Keller, Phys. Rev. A 102, 032616 (2020).
  * Meraner  _et al._ [2020] M. Meraner, A. Mazloom, V. Krutyanskiy, V. Krcmarsky, J. Schupp, D. A. Fioretto, P. Sekatski, T. E. Northup, N. Sangouard, and B. P. Lanyon, Phys. Rev. A 102, 052614 (2020).
  * Steiner  _et al._ [2013] M. Steiner, H. M. Meyer, C. Deutsch, J. Reichel, and M. Köhl, Phys. Rev. Lett. 110, 043003 (2013).
  * Brandstätter  _et al._ [2013] B. Brandstätter, A. McClung, K. Schüppert, B. Casabone, K. Friebe, A. Stute, P. O. Schmidt, C. Deutsch, J. Reichel, R. Blatt, and T. E. Northup, Rev. Sci. Instrum. 84, 123104 (2013).
  * Takahashi  _et al._ [2017] H. Takahashi, E. Kassa, C. Christoforou, and M. Keller, Phys. Rev. A 96, 023824 (2017).
  * Lee  _et al._ [2019b] M. Lee, M. Lee, S. Hong, K. Schüppert, Y.-D. Kwon, T. Kim, Y. Colombe, T. E. Northup, D.-I. D. Cho, and R. Blatt, Phys. Rev. Applied 12, 044052 (2019b).
  * Ong  _et al._ [2020] F. R. Ong, K. Schüppert, P. Jobez, M. Teller, B. Ames, D. A. Fioretto, K. Friebe, M. Lee, Y. Colombe, R. Blatt, and T. E. Northup, New J. Phys. 22, 063018 (2020).
  * Kersevan [1991] R. Kersevan, Technical Report ST/M-91/17, Sincrotrone Trieste Technical Report (1991).

##  Appendix A rf-null position shift

Figure 10:  (a) Contour plot of rf quadrupole potential with the parameters
used in Figs. 3(a)–(f). (b) Central area. Red dashed lines are $x$ and $y$
axes. It is notable that rf null is shifted to $+y$ direction slightly.

##  Appendix B Blade mismatch along the radial direction

We investigate the influence of a mismatch in which one of the rf blades is
transversally displaced from the original position. Fig. 11(a) shows that the
rf blade at $x<0$ and $y>0$ is translated by $r-d/2$. We numerically calculate
$\phi(\textbf{r})_{\rm{rf}}$ as $r-d/2$ changes from $200$ $\mu$m to $-200$
$\mu$m, show the results in Figs. 11(b)–(f), and find rf-null points in Fig.
11(g). When $r$ is larger than the original position $d/2$, the rf null moves
toward the position of $x<0$ and $y>0$, and the tendency becomes opposite at
$r<d/2$.

We remark the similarity/difference between this transversal misalignment and
the tilt misalignment in Sec. II.7. In Fig. 7, the cross section of
$\phi(\textbf{r})_{\rm{rf}}$ and the rf-null positions are plotted as the $z$
coordinate varies. The behavior in a plane of Fig. 7(a) is very close to one
result of Fig. 11. Given $z$ position where the distance between the rf blades
is larger than $d$, the potential shape and shift of the rf null is similar
with the results at $r-d/2>0$ of Fig. 11; when the distance between the rf
blades is smaller than $d$ (Fig. 7), the similarity is found with the result
at $r-d/2<0$ (Fig. 11). Regarding the difference, in this transversal
mismatch, the rf tube and trap axis (defined by $\phi(\textbf{r})_{\rm{tot}}$)
run parallel, which does not cause an axial micromotion.

Figure 11:  (a) Rf blade at $x<0,y>0$ is translated transversally by a
distance $r-d/2$. Blade-blade distance is $d=500$ $\mu$m. Semi-transparent
plane denotes the position of cross sections shown in (b)–(f). Contour plot of
pseudopotential at $r-d/2=200$ $\mu$m in (b), $r-d/2=100$ $\mu$m in (c),
$r-d/2=0$ $\mu$m in (d), $r-d/2=-100$ $\mu$m in (e), and $r-d/2=-200$ $\mu$m
in (f). (g) Coordinates of rf null are plotted as functions of $r-d/2$ and
$r$.