# Elevating zero dimensional global scaling predictions to self-consistent
theory-based simulations

T. Slendebroek<EMAIL_ADDRESS>Oak Ridge Associated Universities,
Oak Ridge, TN 37830, USA  J. McClenaghan  General Atomics, San Diego,
California 92121, USA  O.M. Meneghini  General Atomics, San Diego, California
92121, USA  B.C. Lyons  General Atomics, San Diego, California 92121, USA
S.P. Smith  General Atomics, San Diego, California 92121, USA  T.F. Neiser
General Atomics, San Diego, California 92121, USA  N. Shi  General Atomics,
San Diego, California 92121, USA  J. Candy  General Atomics, San Diego,
California 92121, USA

###### Abstract

We have developed an innovative workflow, STEP-0D, within the OMFIT integrated
modelling framework. Through systematic validation against the International
Tokamak Physics Activity (ITPA) global H-mode confinement database, we
demonstrated that STEP-0D, on average, predicts the energy confinement time
with a mean relative error (MRE) of less than 19%. Moreover, this workflow
showed promising potential in predicting plasmas for proposed fusion reactors
such as ARC, EU-DEMO, and CFETR, indicating moderate H-factors between 0.9 and
1.2. STEP-0D allows theory-based prediction of tokamak scenarios, beginning
with zero-dimensional (0D) quantities. The workflow initiates with the PRO-
create module, generating physically consistent plasma profiles and
equilibrium using the same 0D quantities as the IPB98(y,2) confinement
scaling. This sets the starting point for the STEP (Stability, Transport,
Equilibrium, and Pedestal) module, which further iterates between theory-based
physics models of equilibrium, core transport, and pedestal to yield a self-
consistent solution. Given these attributes, STEP-0D not only improves the
accuracy of predicting plasma performance but also provides a path towards a
novel fusion power plant (FPP) design workflow. When integrated with
engineering and costing models within an optimization, this new approach could
eliminate the iterative reconciliation between plasma models of varying
fidelity. This potential for a more efficient design process underpins
STEP-0D’s significant contribution to future fusion power plant development.

††preprint: AIP/123-QED

##  I Introduction

Viable fusion power plants require a burning or ignited plasma core to produce
a net positive electrical output. Experimental devices, up to this point, have
achieved production of 3 MW of fusion power over the span of 4 seconds [1].
Recent JET deuterium-tritium experiments have demonstrated enhanced fusion
performance, achieving 10 MW for 5 seconds [2], by optimizing the deuterium-
tritium ratio, which in turn increases the beam-target fusion rate. Despite
these noteworthy accomplishments, the transition from short-lived fusion power
bursts to sustainable operation over extended periods is far from
straightforward. Operating in a pulsed manner for approximately three hours to
steady-state, powered by thermalized fusion reactions with a net electric
output, remains an ongoing challenge. Various approaches have been proposed to
meet this challenge such as larger size like EU-DEMO [3], higher magnetic
field approaches such as ARC [4] or advanced plasma operation approaches, such
as CAT [5]. The space of possible tokamak configurations, in terms of
engineering and plasma parameters, is vast and expands as advances are made in
technology such as the improvements in magnetic field with the use of high
temperature superconductors and advances in understanding of the plasma core
through high confinement regimes. How do these advances and improvements
translate to the performance of the plasma core? 0D systems codes [6, 7] are
frequently used to scope the wide and multi-dimensional parameter space of
tokamaks. These 0D systems codes predict the performance of the core plasma by
relying on scaling laws, notably the IPB98(y,2) confinement scaling [8]. This
scaling law uses 8 physically meaningful parameters that describe the plasma
and geometry of the tokamak and is derived from a weighted least squares
regression of the International Tokamak Physics Activity (ITPA) global H-mode
confinement database as a power law. The use of this power law is robust and
accurate for the given data range. However, the parameter space of a fusion
power plant (FPP) lies outside the database range and thus extrapolations of
the scaling law are required which adds uncertainty which is difficult to
quantify.

In contrast to extrapolation with scaling laws, predictions based on
integrated modeling workflows remain valid provided that the underlying
physics-based models are used within their applicable bounds, which extend to
the FPP operational domain. For example, the flux-driven transport solver
TGYRO [9] and pedestal stability code EPED [10] have been successfully used in
integrated modeling workflows for both existing experiments and future burning
plasma devices [[11], [12],[13]]. Notably, initial performance predictions of
burning plasmas in SPARC [14] using a reduced model of gyrokinetic transport
(TGLF [15]) have recently been successfully validated using high-fidelity
gyrokinetic simulations with CGYRO [16] in a novel integrated modeling
workflow [[17], [18]]. Integrated modeling has been used to reach a self
consistent 1.5D solution by iterating over core sources, transport, pedestal
and equilibrium codes. Predictions of the self-consistent simulations on
present day devices have shown good agreement with experimental results [[19],
[20], [21]]. A key aspect of these integrated simulations is the initial
conditions from which the iteration process is started. For example, a widely
used starting point to predict the performance of existing devices is to begin
from a past experimental condition and then modify original actuator settings
[22]. Another approach is to start from databases of existing simulations
[23]. Both of these methods require significant effort to finding data either
experimental or simulated and are limited in their ability to rapidly explore
a wide range of parameters.

In this paper we present a method for bootstrapping 1.5D stationary
simulations of fusion plasmas, starting from the scalar (i.e. 0D) engineering
parameters that are commonly used for cross-machine studies and extrapolation
to fusion power plants using 0D systems codes (see section I-A). We used this
approach to efficiently initialize self-consistent 1.5D simulations (see
section I-B). We validate this approach against DIII-D shots (see section
II-A) and existing cross-machine databases such as the ITPA global H-mode
confinement database [24] from which the IPB98(y,2) confinement scaling [8] is
derived (see section II-B). Predictions of proposed fusion reactors were made
in section II-C. We also summarize this research and provide perspective on
future work (see Section III).

###  I.1 Bootstrapping a 1.5D simulation from 0D parameters with PRO-create

Figure 1: Zero dimensional engineering parameters are transformed by PRO-
create (PC1 & PC2) into an initial guess for the plasma configuration, which
are evolved to a self-consistent, stationary physics solution by STEP.

_Ab initio_ 1.5D modeling of tokamak plasmas requires a starting equilibrium
solution, kinetic profiles (species densities, temperatures, and rotation),
and sources profiles (heating, current drive, and torque). A new module named
PRO-create has been created to generate a physically plausible (though not
necessarily self-consistent) starting point for 1.5D modeling in the OMFIT
framework [25]. Figure 1 illustrates the workflow from zero dimensional
engineering parameters to the fully self-consistent stationary 1.5D solution.
PRO-create has two stages, the first equilibrium stage PC1 is initialized with
a simple equilibrium using a constant pressure gradient and a polynomial
current density profile, and a second kinetic pressure profile stage PC2, with
bootstrap and current drive contributions to the current density profile.

In the first stage PC1, PRO-create initializes with the fixed boundary
equilibrium solver CHEASE [26], assuming a constant pressure gradient
$\frac{dP}{d\psi_{N}}$, where $P$ is the pressure and $\psi_{N}$ is the
normalized poloidal flux. The plasma boundary can be given as any arbitrary
closed boundary. The graphical user interface (GUI) of PRO-create’s
parameterization of the fixed boundary is shown in Figure 2. The interactive
boundary shape tool implemented in PRO-create allows the user to extend the
shape to include X-points, up-down asymmetry and squareness. The current
density profile term labeled PC1 is a polynomial of the form:

$\begin{split}<J_{\phi}/R>(\psi_{N})=&c_{I}\bigg{(}\underbrace{I_{p}(1-\psi_{\rm
N})^{\alpha}}_{\text{PC1 stage}}+\\\
&\underbrace{<J_{BS}/R>(\psi_{N})+<J_{CD}/R>(\psi_{N})}_{\text{PC2
addition}}\bigg{)}\end{split}$ (1)

where $<J_{\phi}/R>$ is the flux surface average of the toroidal current
density over the major radius, $I_{p}$ is the total plasma current, $c_{I}$ is
a scalar to ensure the integral of averaged $<J_{\phi}/R>$ over the area is
equal to the plasma current, and $\alpha$ is a scalar that controls the
broadness of the current profile (defaulted to $\alpha=1.4$,  this is only an
initial guess and not consequential for the end result see reference [27] for
a discussion about the parameter). In the term labeled PC2, the $J_{BS}$ is
calculated with the NEO 2021 adaption of the Sauter formula [28] and the
$<J_{CD}/R>$ term denotes the contribution of current drive from the auxiliary
heating system, and is calculated using equation 10.

Figure 2: Graphical User Interface of PRO-create’s interactive Boundary
parameterization.

In the second stage PC2, H-mode kinetic profiles are then defined on the
normalized $\psi_{N}$ grid for the next equilibrium iteration. The initial
temperature profiles have the form of a core polynomial and a hyperbolic
tangent at the pedestal:

$T_{e}(\psi_{\rm N})=T_{i}(\psi_{\rm N})=T_{1}(\psi_{\rm N})+T_{2}(\psi_{\rm
N})$ (2)

With $T_{1}$ and $T_{2}$:

$T_{1}(\psi_{\rm N})=\frac{1}{2}c_{t}\left(\tanh\left(\frac{\psi_{\rm
N}-x_{\text{ped}}}{w_{\text{ped}}}\right)+\tanh(1)\right)$ (3)

$T_{2}(\psi_{\rm N})=H(1-\overline{\psi}_{\rm
N})(T_{\text{core}}-T_{1}(\psi_{\rm N}=0))\left(1-\overline{\psi}_{\rm
N}^{a_{\text{in}}}\right)^{a_{\text{out}}}.$ (4)

Here, $H$ the heavy-side step function,
$c_{t}=\frac{T_{\text{ped}}-T_{\text{edge}}}{\tanh(1)}$,
$T_{\text{core}},T_{\text{ped}},T_{\text{edge}}$ the core, pedestal, and edge
temperatures respectively; $\overline{\psi}_{\rm N}=\psi_{\rm
N}/(x_{\text{ped}}-w_{\text{ped}})$, with $x_{\text{ped}}=1-w_{\text{ped}}$,
where $w_{\text{ped}}$ is the pedestal width as defined by [10] (code that
determines the pedestal height that conforms to the peeling and ballooning
stability boundary). The core polynomial exponents for the temperature profile
are $a_{\text{in},T}=1.2,$ $a_{\text{out},T}=1.4$, the same profile shape is
used for the density but with $a_{\text{in},\,ne}=1.1$ and
$a_{\text{out},\,ne}=1.1$. These parameters are set to these values to be
consistent with the EPED1 parameterization see reference [27]. The initial
pedestal height is obtained from the EPED-NeuralNet [29] (EPED-NN, neural net
version of EPED1). The ion densities are set by enforcing quasi neutrality as
well as a given $Z_{\text{eff}}$ profile. For a two species system of one
fully ionized main impurity species and one fuel species the ion densities are
expressed as:

$n_{\text{impurity}}(\psi_{N})=n_{e}(\psi_{N})\left(\frac{Z_{\text{eff}}(\psi_{N})-1}{Z_{\text{impurity}}^{2}-Z_{\text{impurity}}}\right)$
(5)

$n_{\text{fuel}}(\psi_{N})=n_{e}(\psi_{N})-Z_{impurity}\,n_{\text{{impurity}}}(\psi_{N})$
(6)

where $n_{\text{impurity}}$ and $z_{\text{imp}}$ are the impurity density and
charge state respectively and $n_{\text{fuel}}$ is the fuel species density
(hydrogen / deuterium / tritium).

The broad initial guess for the angular frequency (rotation) profile is set
as:

$\omega_{0}(\psi_{N})=0.9\omega_{0\,\text{axis}}e^{-3\psi_{N}}+0.1\omega_{0\,\text{axis}}$
(7)

where $\omega_{0\,\text{axis}}$ is the angular frequency (rotation) in rad/s
at the magnetic axis.

With the kinetic pressure profile and the complete current density profile
described in equation 1 a second iteration of CHEASE is taken to find the
equilibrium consistent with the kinetic pressure profile and current density
profile.

Table 1: Heating and current drive parameter assumptions Heating system | $f_{cd}$ | $f_{e}=P_{\text{source}}/P_{e}$ |  $\text{width}_{e}$ / $\text{width}_{i}$ |  $n_{ge}$ / $n_{gi}$
---|---|---|---|---
Neutral Beam Injection | 0.18 |  $\text{Sivukhin}(\rho)$ [30] | 0.27/0.2 | 1.0/0.7
Electron Cyclotron Heating | 0.2 | 1.0 | 0.15 $\frac{0.606[m]}{\text{minor radius}[m]}$ | 1
Ion Cyclotron Heating | 0.125 | 0.25 | 0.15 $\frac{0.606[m]}{\text{minor radius}[m]}$ | 1

Idealized heating and current drive profiles are implemented for neutral beam
injection (NBI), electron cyclotron heating (ECH), ion cyclotron heating
(ICH), and lower-hybrid heating (LH) systems each with different current drive
efficiencies and ion-to-electron energy deposition ratio. The heating profiles
are defined by equation 8 as Gaussians of nth order. These heating profiles
are required during the flux-matching in the self-consistent evolution in
STEP.

$q_{e}(\rho)=c_{iq}f_{e}(\rho)P_{\text{source}}\,\,\exp^{\left(\frac{\rho-\rho_{0}}{w_{0}}\right)^{ng_{e}}}$
(8)

$q_{i}(\rho)=q_{e}(\rho)(1-f_{e}(\rho))$ (9)

Here $c_{iq}$ is a volume density scalar to ensure the volume integrated
electron energy density source $q_{e}$ is equal to the absorbed power of the
heating system $P_{\text{source}}$, $f_{e}$ is the electron fraction of the
total power, $w_{0}$ the Gaussian width parameter, $\rho_{0}$ the normalized
position of the center of the peak, and $n_{g}$ the order of the Gaussian. The
parameters for each heating system can be found in Table LABEL:table:hcd. The
width and the Gaussian order $n_{g}$ of the NBI parameterization were deduced
from running 50 instances of FREYA [31] (Neutral Beam deposition code) on
DIII-D experimental shots. The width of the ECH and ICH is found by simulation
using TORAY[32] (ECH ray tracing code) and GENRAY [33] (ray tracing code for
ICH and LH) scaled down by minor radius to capture the tokamak size aspect.
The current drive of each heating system assumes the same profile shape as the
energy density source but with a current drive efficiency $\eta_{CD}$ for each
system [34] (chapter 6 of reference [[34]]).

$<J_{CD}/R>(\rho)=c_{ij}\eta_{CD}P_{\text{source}}\,\,\exp^{\left(\frac{\rho-\rho_{0}}{w_{0}}\right)^{n_{ge}}},$
(10)

with $c_{ij}$ an area density scalar to ensure the area integrated current
density is equal to the total current drive of the system. The particle source
and momentum have the same profile as the NBI heating but with particle source
amplitude $S_{NBI}=P_{NBI}/(eE_{NBI})$ with $P_{NBI}$ the absorbed heating
power of the NBI and $E_{NBI}$ the energy of the beam particles and for the
total injected momentum
$T_{NBI}=S_{NBI}\sqrt{\frac{2eE_{NBI}}{m_{beam}}}m_{beam}\sin{\theta_{V}}$
with $\theta_{V}$ the vertical tilt angle of the NBI with respect to the
plasma cross section and $m_{beam}$ the mass of the injected NBI particles.

The output of PRO-create is stored in the IMAS data schema (IDS) [35] format
using the OMAS implementation [36], the equilibrium information is stored in
the _equilibrium_ IDS, the kinetic profiles in _core profiles_ and the heating
and current drive in the _core sources_ IDS.

###  I.2 Self-consistent theory based solution with STEP

The 1.5D solution provided by PRO-create is physically plausible but not yet
consistent so it is used as the starting point for the STEP [37] workflow
which iterates over core sources, transport, pedestal and equilibrium codes
(as illustrated in Figure 3) until it converges on a self-consistent solution.
Each physics step reads and writes data to and from the same centralized data
structure (as illustrated in Figure 4). The physics steps used in this paper
are summed up below:

  * •

For the equilibrium step the fixed boundary code CHEASE [26] is used, CHEASE
makes the pressure and current profiles self-consistent based on the solution
of the grad-shavranov equation using the kinetic $P(\rho)$ and
$<J_{\phi}/R>(\rho)$ as inputs in the _equilibrium_ IDS. The solution provides
the reference magnetic geometry used to map profiles. An anomalous resistivity
was added to keep to keep the safety factor above one, essentially emulating
the sawtoothing process; nevertheless, the influence of sawteeth on
temperature was not considered in this approach.

  * •

For the current evolution, ONETWO [38] is used which calculates the self-
consistent ohmic current profile and bootstrap current profile. ONETWO is
setup to evolve the current till steady state, this assumes the shot has
undergone a few current diffusion time steps from the start of the flattop.

  * •

In order to make the density, temperature and rotation profiles consistent
with the fluxes according to validated [39] reduced models of neoclassical and
turbulent (gyrokinetic) fluxes TGYRO [9] is used. TGYRO calls TGLF SAT2 [15]
with electromagnetic contributions (EM) and NEO [40] for the turbulent and
neoclassical fluxes respectively. The particle, energy and momentum source is
gathered from the _core sources_ IDS to which TGYRO adds the synchrotron and
line radiation as well as the bremsstrahlung. The transport grid is set with 8
points evenly spaced from $0.2\leq\rho\leq 0.8$ with the boundary condition
$\rho=0.8$, this boundary point is set as the solution of EPED-NN [29] which
calculates the width and the height of the pedestal. The pedestal height is
blended with the boundary condition at $\rho=0.8$ see reference [41] for more
details. The solutions of the fluxes are stored in the _core transport_ IDS.

In the TGYRO step the line-averaged density is matched to the experimental or
target value by shifting the pedestal up and down to match the synthetic
interferometer (this is done after each converged TGYRO solution). Each
tokamak has its own interferometer path, for consistency the synthetic
interferometer’s path was chosen as a horizontal path through the magnetic
axis.

In Table LABEL:table:evolution an overview can be seen of the most physically
relevant plasma parameters and profiles that are used and calculated in PRO-
create and STEP.

Figure 3: Specific STEP workflow cycle used in this paper with the following steps: ONETWO:current evolution, CHEASE: equilibrium, TGYRO: pedestal and transport. Figure 4: Data flow of the STEP workflow, initializing with PRO-create to form a starting point in the OMAS data structure implementation where each physics step reads and writes data to the same centralized data structure (adapted from reference [36]). Table 2: Overview of plasma parameters and profiles used in PRO-create and STEP Quantity | Evolved or Static | Description
---|---|---
$I_{p}$ | static : 0D Input | Total toroidal plasma current
$B_{t0}$ | static : 0D Input | Toroidal magnetic field on-axis
$\kappa$ | static: 0D Input | Elongation
$\delta$ | static: 0D Input | Miller up down symmetric triangularity
$Z_{\text{eff}}$ | static: 0D Input | Effective charge Z for all ion species
$n_{e,\text{line}}$ | static : 0D Input | Line averaged density
$P_{\text{source}},S_{\text{source}}$ | static: 0D Input | Total absorbed power and total particle source
$q_{e}(\rho)$ | static : Equation 8 | Electron heating energy density
$q_{i}(\rho)$ | static : Equation 9 | Ion heating energy density
$n_{e,\text{sep}}$ | static : $n_{e\,ped}/4$ | Separatrix density
$q_{e,ohmic}$ | evolved : ONETWO | Ohmic heating profile
$n_{e,ped}$ | evolved : Synth diagnostic | Pedestal density, EPED definition [10]
$n_{e}(\rho)$ | evolved : TGYRO | Electron density profile
$T_{e}(\rho),T_{i}(\rho)$ | evolved : TGYRO | Electron and Ion temperature
$\omega_{0}(\rho)$ | evolved : TGYRO | Angular frequency (rotation)
$<J_{tor}/R>(\rho)$ | evolved : ONETWO | Flux surface average of the toroidal current density over R
$\Psi(R,Z)$ | evolved : CHEASE | 2D flux surfaces map
$\Psi_{\text{bound}}$ | static : Miller parameterization | 2D fixed boundary

Predictions of plasma pulses using this integrated modeling workflow starting
with very limited information requires us to make a number of assumptions, yet
we will show that it yields good predictions for global confinement when
staying within the validated parameter space of underlying physical models and
assumptions.

##  II Results

###  II.1 Sample DIII-D case from experiment

Table 3: 0D input parameters for DIII-D discharge 81499 at 4000 ms $I_{P}[MA]$ | $B_{t}[T]$ | $Z_{\text{eff}}$ | $n_{e,\text{line}}[10^{20}\,m^{-3}]$ | $R_{\text{major}}[m]$ | $a_{\text{minor}}[m]$ | $\kappa$ | $\delta$ | $P_{NBI}[MW]$ | $E_{NBI}[keV]$ |
---|---|---|---|---|---|---|---|---|---|---
1.34 | 1.91 | 2.32 | 0.48 | 1.69 | 0.63 | 1.68 | 0.32 | 5.74 | 72.6 |

The PRO-create/STEP workflow has been validated on a number of DIII-D
discharges. Shot # 81499 is shown as an example of this validation set.
Information about this shot is publicly available [42]. This shot is a
conventional H-mode with global parameters as described in Table
LABEL:table:0dparams. Information about this shot includes plasma profiles
like the density, rotation and temperature profiles but also heating and
current drive profiles from interpretative analysis in addition to the EFIT01
[43] plasma boundary (EFIT01 is a magnetic only reconstruction of the
equilibrium using EFIT and is valid for the plasma boundary shape). Setting up
STEP from the experimental equilibrium including experimental heating and
current profiles, is referred to as STEP-experimental. This approach is
compared to setting up STEP based on 0D quantities, which is referred to as
STEP-0D. The ITPA global H-mode database records the plasma boundary in up-
down symmetric Miller geometry while the experiment is a lower single-null
plasma with up-down asymmetric geometry. To reconcile the difference in plasma
volume between STEP-0D and STEP-experiment the symmetric Miller elongation of
STEP-0D is changed to match the recorded plasma volume from the database using
a simple optimization algorithm. This is important when comparing volume based
quantities like the thermal stored energy. The pedestal model EPED and TGLF
both utilize the symmetric Miller parameterization, thereby neglecting to
account for the asymmetry caused by the x-point. The equilibrium shape of shot
# 81499 from STEP-0D with adjusted elongation compared to the EFIT01 and STEP-
experimental can be seen in Figure 5. The predicted plasma density and
temperatures of STEP-0D and STEP-experimental can be seen in Figure 6. The
electron density and temperature profiles are nearly identical with some
slight difference in the ion temperature gradient scale length around
$0.3<\rho<0.5$. This discrepancy can be explained by contrasting the
analytical heating and current drive applied in STEP-0D with the actual
heating profiles in the STEP-experimental run; specifically, the former
demonstrates a higher electron-to-ion heating ratio, while the latter exhibits
more off-axis heating. The matching pedestal density indicates that matching
the experimental line averaged density from the experiments using the
synthetic diagnostic in the TGYRO step works as intended.

Figure 5: Plasma equilibrium of shot 81499 at 4000 ms for the experimental
EFIT01 in black, STEP-experimental in orange and STEP-0D in blue. Figure 6:
Core profile comparison between experiment (black), STEP-experimental (orange)
and STEP-0D (blue) for DIII-D shot 81499 at 4000 ms, (top) electron density,
(middle) electron temperature and (bottom) ion temperature.

###  II.2 Global H-mode confinement database

To extend the validation of the STEP modeling beyond single experimental use
cases STEP-0D integrated modeling has been applied to a subset of the ITPA
global H-mode database. This database contains data from multiple machines
with confinement times ranging from milliseconds to 1200 milliseconds. The
database was filtered to exclude elongation smaller than 1.4 ($\kappa>1.4$ as
the EPED-NN training’s data does not contain small elongation), filtering
shots that do not have deuterium as main fuel species (to exclude Hydrogen and
Helium plasmas) and filtering plasmas that did not record triangularity.

This subset contains 600 discharges from 7 different tokamaks with carbon
walls and divertor, with an exception of CMOD which has a molybdenum wall
(fully ionized Neon was used to model the CMOD impurities). For all the cases
in the database the PRO-create and STEP settings were unchanged and the STEP
workflow steps through the cycle described in Figure 3 three times. Although
PRO-create and STEP support up-down asymmetric plasma boundaries the data
recorded in the ITPA H-mode database includes only the up-down symmetric
Miller quantities. The same approach as in Section II was used by matching the
total volume by changing the elongation slightly. The database includes the
line averaged density but no information about the pedestal or separatrix
density therefore PRO-create sets an initial guess for these quantities as
fractions of the line averaged density,
$n_{e\,\text{ped}}=n_{e\,\text{line}}/1.3$,
$n_{e\,\text{sep}}=n_{e\,\text{ped}}/4$. A flat profile is prescribed,
$Z_{\text{eff}}(\rho)=Z_{\text{eff}}$ where $Z_{\text{eff}}$ comes from the
database. If this value is missing $Z_{\text{eff}}=2$ is chosen. The initial
on axis rotation is set as
$\omega_{0\text{axis}}=10^{5}\,[\text{rad}/\text{s}]$ with an exception for
shots without any NBI. For those the rotation is not evolved and set to zero.
The heating and current drive profiles are updated throughout the run with the
analytic form as presented in section A heating and current drive codes
weren’t used as there is not enough information in the database to set-up each
heating system so core deposition at $\rho_{0}=0.0$ was assumed for all cases.
EPED predicts the peak pressure of the ELM cycle, it is unknown at what point
in the cycle the data is taken therefore the pedestal height is reduced to an
effective 80 % of the EPED(-NN) prediction. This reduction is supported by
experimental evidence as can be seen in Figure 7 where 76728 EPED-NN
predictions were made and compared to automated tanh function fitting to
Thomson scattering measurements in real space for 655 DIII-D shots where the
mean ratio of experimental to predicted sits at $p_{\text{e, ped,
exp}}/p_{\text{e, ped, EPED-NN}}=0.794$.

Figure 7: The histogram shows the ratio between the experimental pedestal
electron pressure and the predicted electron pedestal using EPED-NN. The mean
and median of the ratio are 0.794 and 0.758, respectively. The 76728
experimental data points were obtained by automated tanh function fitting to
Thomson scattering measurements in real space for 655 DIII-D shots.

To compare STEP-0D to the experiments, the experimental energy confinement
time is plotted against the STEP-0D simulated energy confinement time as a
regression plot in Figure 8. The mean relative error of the STEP-0D results to
the experiments for this subset is 19%, for context the IPB98(y,2) scaling law
has a mean relative error of 22% for the same subset. The discharges from AUG,
CMOD, DIII-D JET and JT60 follow the
$\tau_{E,\text{exp}}=\frac{W_{th,\text{exp}}}{P_{lth,\text{exp}}}=\tau_{E,\text{STEP}}=\frac{W_{th,\text{STEP}}}{P_{lth,\text{exp}}}$
line (with $P_{lth,\text{exp}}$ the experimentally defined loss of power of
the thermal species) with some scatter. The three COMPASS experimental points
with the shortest confinement times received a high amount of injected power.
Both the EPED and EPED-NN models predict significantly higher pedestal
pressures (6-12 kPa) in these plasmas compared to the predictions made by Komm
et al. [44]. One plausible explanation for the difference in $\tau_{e}$
between the experimental results and the STEP predictions is that the
pedestals in the experimental plasmas were constrained by smaller Type III
ELMs [45] driven by ballooning instabilities, which are not accounted for in
EPED. To support this claim the importance of the pedestal height on the
confinement time in the pedestal stored energy ($W_{ped}$) component is
compared to the core stored energy ($W_{core}$) component by splitting them
up, as is illustrated in Figure 9. Although the core gradients and pedestal
height are not entirely independent, the pedestal height plays a significant
role in determining both the total stored energy and the height of the core
profiles. The ratio of total stored thermal energy ($W_{th}$) to the pedestal
contribution can be seen in Figure 10, which concludes that the pedestal
contribution is significant for conventional type I ELMy H-mode shots.

Figure 8: Experimental confinement time compared to the STEP-0D result for the
subsection of the ITPA database, the black dotted lines denote a 50% error
from the experimental energy confinement time, the mean relative error of the
STEP-0D results to the experiments for this subset is 19%. Figure 9:
Visualization of the Pedestal (magenta) and core contribution (purple) to the
thermal stored energy. Figure 10: Histogram of the ratio of total thermal
stored energy to the pedestal stored energy for the STEP results of the
validation database subset. Figure 11: Radar plot showing the exponents of
weighted linear regressions for the scaling law tau IPB98(y,2) (blue), the
validation subset (black) and the STEP-0D result (orange), all exponents are a
positive contributions with the exception of the loss power.

Currently there is no equivalent pedestal model for type III ELM shots as well
as for ELM free operational regimes like negative triangularity, EDA-H, grassy
ELM regime. Exploration of ELM-free discharges [46] is important for fusion
power plant sized devices as ELM heat fluxes scale unfavorably in fusion power
plant parameter spaces. The formulation of such models is therefore critical
in performing any predictive simulation and scoping of parameter spaces for
fusion power plant operation. STEP would prove a useful tool for coupling
these new type of edge models to benchmark against experiment and improve the
models.

Weighted least squares regression (WLS) of the experimental database and STEP
results can be made to infer differences between the IPB98(y,2) scaling law
and modeling results. The entire database used to generate the IPB98(y,2)
scaling consists of 10 000 shots, its WLS exponents can be seen in blue in
Figure 11 where the weight is a function of the number of plasma discharges in
the database $w_{i}=1/(2+\text{Int}(\sqrt{N_{i}}/4))$(as described in McDonald
et al. [47]). A WLS on the subset used in this study gives the exponents in
black. The best fitting exponents for the WLS regression ($R^{2}=0.998$) of
the STEP results are shown in orange in Figure 11. The exponents of the WLS
regression of STEP indicate agreement with the data subset with the exception
to elongation $\kappa$ and the inverse aspect ratio $\epsilon$. Emphasis
should not be given to the exact exponents of the STEP power law, as for
example sensitivity analysis of the pedestal density (Figure 12) on the
pedestal pressure indicates that the underlying physical behavior is not
power-law like but depends on the operating regime. Peeling limited pedestal
pressures scale inversely with pedestal density in the ballooning limited
regime whereas peeling limited regimes scale with the pedestal density.
Additionally it is difficult to attribute what causes the difference in power
law exponents like elongation and inverse aspect ratio as they are heavily
correlated to the other parameters. Similar analysis but then for L-mode
plasmas is done see reference [48]. This successful validation of the PRO-
create STEP workflow for ELMy H-mode plasmas gives confidence in predicting
the fusion power output of fusion power plant sized devices. Utilization of
this workflow finds discrepancies between 0D-system code predicted fusion
power and the predicted fusion power (see reference [[37]] for an application
of this). Another ongoing effort of this STEP - PRO-create workflow is the
design of a plasma scenario for a new tokamak device called BEST (Burning
plasma Experimental Reactor).

Figure 12: EPED neural net sensitivity scan, showing the pedestal pressure as
a function of the pedestal density where the EPED-NN input parameters are kept
constant except for the pedestal density for a DIII-D nominal operation point
(dashed lines), the shaded region in red denotes the uncertainty of the model.

###  II.3 Confinement predictions for fusion reactors

With the validated STEP workflow starting from 0D parameters predictions can
be made for proposed fusion reactors. To demonstrate this capability
predictions were made for CFETR (Table 1 baseline case [49]), ARC (ARC RPL
mode Table 2 right most column [50]) and EU-DEMO (Table 1 EU-DEMO 2018 [3]).
The 0D parameters of each case are summarized in table 4. The plasma rotation
profile $\omega_{0}(\rho)$ for all cases was set to zero due to the large size
device in the case of CFETR and EU-DEMO and for ARC zero torque from H&CD
systems. A Helium specie was included in all cases as 0.03 fraction of
electron density. The same settings for each of the physics codes were used
(see section B), EPED was ran instead of EPED-NN as these cases lie outside
the NN trainings dataset. For TGLF the SAT 2 model was used including
electromagnetic contributions where TGYRO evolved the density and temperatures
while keeping the pedestal density fixed. The predicted plasma performance and
confinement factors can be seen in Figure 13.

Table 4: 0D input parameters for CFETR, ARC, EU-DEMO where –> denotes changes to the reference parameters. Device | $I_{P}[MA]$ | $B_{t}[T]$ | $Z_{\text{eff}}$ | $n_{e,\text{ped}}[10^{20}\,m^{-3}]$ | $R_{\text{major}}[m]$ | $a_{\text{minor}}[m]$ | $\kappa$ | $\delta$ | $P_{aux}[MW]$ | ion species
---|---|---|---|---|---|---|---|---|---|---
CFETR | 13.0 | 6.5 | 2.2 | 0.7 | 7.2 | 2.2 | 2.01 | 0.43 | 80 | D,T,He,Ar
ARC | 14 | 11.5 | 1.13 | 2.5 | 4.2 | 1.2 | 1.62 | 0.35 | 22.7 | D,T,He,Kr
EU-DEMO | 17.75 | 5.86 | 2.19 | 0.45 | 9.07 | 2.92 | 1.65 | 0.33 | 50 –> 80 | D,T,He,Kr

#### ARC

In the ARC reference, the simulation simulates a lower-hybrid heating of 22.7
MW, employing the same heating profiles utilized in the STEP analysis. This
was done by setting the transport solver’s boundary point at $\rho=0.9$ with
$T_{i}=T_{e}$ = 2.5 keV. which was derived through scaling L-mode plasmas from
CMOD. This method was in line with the accurate replication of the L-mode
scaling law, as reported by Angioni in 2023 [48]. However, this assumption may
be overly optimistic for an L-mode since analysis with EPED suggests that the
peeling-ballooning limit is at 147 kPa at the pedestal (220 kPa at
$\rho=0.9$). Furthermore, the resulting fusion power from the simulation is
particularly sensitive to the transport boundary condition at $\rho=0.9$. This
sensitivity is apparent in the results presented in Table 5. Here, the fusion
power in H-mode at the peeling-ballooning limit was comparable to that in the
reference scenario. Yet, when considering 80% or less of the peeling-
ballooning limit pressure, the profiles collapse due to the loss of fusion
heating and the highly radiative krypton species. The detailed transport
analysis of the L-mode edge is thus crucial to ascertain the validity of the
reference scenario, considering the high sensitivity of the fusion power to
the boundary condition. Moreover, the H-mode peeling-ballooning limit might be
overly optimistic, given that the plasma power falls far short of the
threshold required for the L to H-mode transition scaling.

Table 5: Sensitivity of the fusion power to the transport boundary condition at $\rho=0.9$ where the H-mode peeling ballooning (PB) limit sits at 224 kPa and reference at 200 kPa. cases | PB-limit | Reference | 80% of PB-limit
---|---|---|---
P($\rho=0.9$) [kPa] | 224 | 200 | 181
Fusion power [MW] | 420 | 290 | 2.3

#### CFETR

The CFETR case employs 30 MW of NBI with a beam particle energy of 1 MeV,
coupled with 50 MW of ECH. The NBI deposition is at $\rho_{0}=0.0$ and the EC
deposition is at $\rho_{0}=0.3$. Under these conditions, the workflow yields a
fusion power output of 583 MW. A higher fusion power output of 952 MW, as
suggested by the reference, is attainable under a hybrid operation scenario.
This scenario is predicated on optimizing the plasma current profile, which
enhances overall plasma performance. While it is feasible to conduct similar
optimizations to reach this power output, such measures fall beyond the scope
of the current study.

#### EU-DEMO

In the 2018 scenario for EU-DEMO, a specific heating and current drive system
was not selected, leading to the adoption of 30 MW of NBI with a beam particle
energy of 1 MeV. Concurrently, 20 MW of ECH was implemented with deposition
positions at $\rho_{0}=0.0$ and $\rho_{0}=0.3$ respectively. The pedestal
density was assumed to be 80% of the pedestal top as specified by the case.
The reference fusion power was projected by a 0D system code to be
approximately 2012 MW. However, when 50 MW of auxiliary heating was applied,
the workflow only predicted a fusion power of 550MW, indicating an under-
powered system. An increase in auxiliary power to 80 MW, with 50 MW from NBI
and 30 MW from ECH, resulted in a modest increase in fusion power to 802 MW.
This is notably less than the expectations set by the 0D design study cited as
a reference. While further optimization is beyond the scope of the current
study, the gap to achieve 2000 MW of fusion power appears substantial,
suggesting the necessity of substantial revisions or improvements.

Figure 13: Fusion reactor plasma performance results from the STEP-0D
simulations with in blue $Q_{\text{plasma}}=P_{\text{fusion}}/P_{\text{aux}}$
and $H_{98y2}=\tau_{E,\text{STEP}}/\tau_{E,\text{IPB98y2}}$ in red. The EU-
DEMO case is with the increased auxiliary heating to 80 MW and the ARC case is
at the result at the peeling ballooning limit.

##  III Conclusions

A new STEP workflow has been developed that predicts plasma profiles and
equilibrium starting from 0D tokamak engineering parameters. This workflow
consists of generating physically plausible plasma profiles consistent with
the plasma equilibrium using PRO-create followed by iteration between physics
models for the equilibrium, sources, core transport, and pedestal to obtain a
self-consistent solution using STEP. Systematic validation against the ITPA
global H-mode confinement database has shown that STEP-0D can predict the
energy confinement time with an average mean relative error (MRE) of less than
19%. The successful validation of the PRO-create STEP workflow for ELMy H-mode
plasmas gives confidence in predicting the fusion power output of fusion power
plants. Predictions of this workflow were made on CFETR, ARC and EU-DEMO where
Conventional H-modes with moderate H-factors between 0.9 and 1.2 were found.
The conventional methodology employed in the design of fusion power plants
(FPP) starts with the identification of a operating point that is found by
running a 0D system codes optimization. The chosen design is then further
refined using more sophisticated, high-fidelity models. However, this method
proves to be inefficient at times as the results of high-fidelity models do
not always coincide with the optimal solutions identified by the 0D system
codes. This discrepancy becomes particularly evident when the design deviates
from the assumptions underpinning the system code’s scalings and the reduced
set of equations, such as those valid for a given aspect ratio or operating
regime. This paper presents a viable alternative: the ability to obtain self-
consistent 1.5D physics-based stationary solutions that start from the same 0D
input parameters utilized by the system codes. We envision that this
capability, when integrated with engineering and costing models in an
optimization workflow, could form the foundation for a novel FPP design
workflow. This new workflow would eliminate the need for iterating between
plasma models of varying fidelity and thus avoid wasting time reconciling
inevitable discrepancies. Further development of the physics capabilities of
this work involves coupling different pedestal models to include regimes of
ELM-free operation and developing and validating a robust workflow for
advanced tokamak scenarios (that rely on large radius internal transport
barriers to achieve high confinement and steady-state operation).

### ACKNOWLEDGEMENTS

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Fusion Energy Sciences, using the DIII-D National
Fusion Facility (DE-FC02-04ER54698), a DOE Office of Science user facility,
under Awards DE-SC0017992 (AToM) and DE-FG02-95ER54309 (GA Theory Grant). Part
of the data analysis was performed using the OMFIT integrated modeling
framework [25]. The Authors thank Cihan Akcay for editorial suggestions.

### DISCLAIMER

This report was prepared as an account of work sponsored by an agency of the
United States Government. Neither the United States Government nor any agency
thereof, nor any of their employees, makes any warranty, express or implied,
or assumes any legal liability or responsibility for the accuracy,
completeness, or usefulness of any information, apparatus, product, or process
disclosed, or represents that its use would not infringe privately owned
rights. Reference herein to any specific commercial product, process, or
service by trade name, trademark, manufacturer, or otherwise, does not
necessarily constitute or imply its endorsement, recommendation, or favoring
by the United States Government or any agency thereof. The views and opinions
of authors expressed herein do not necessarily state or reflect those of the
United States Government or any agency thereof.

## References

  * Keilhacker, M and Watkins, ML and JET Team [1999] Keilhacker, M and Watkins, ML and JET Team, “D-T experiments in the JET tokamak,” Journal of Nuclear Materials 266, 1–13 (1999).
  * Maslov  _et al._ [2022] M. Maslov, E. Lerche, C. Challis, J. Hobirk, A. Kappatou, D. King, D. Keeling, F. Rimini, E. De La Luna, I. Monakhov, _et al._ , “T-rich scenario for the record fusion energy plasma in JET DT,” Bulletin of the American Physical Society (2022).
  * Siccinio  _et al._ [2022] M. Siccinio, J. Graves, R. Kembleton, H. Lux, F. Maviglia, A. Morris, J. Morris, and H. Zohm, “Development of the plasma scenario for EU-DEMO: Status and plans,” Fusion Engineering and Design 176, 113047 (2022).
  * Sorbom, BN and Ball, J and Palmer, TR and Mangiarotti, FJ and Sierchio, JM and Bonoli, P and Kasten, C and Sutherland, DA and Barnard, HS and Haakonsen, CB et al. [2015] Sorbom, BN and Ball, J and Palmer, TR and Mangiarotti, FJ and Sierchio, JM and Bonoli, P and Kasten, C and Sutherland, DA and Barnard, HS and Haakonsen, CB et al., “ARC: A compact, high-field, fusion nuclear science facility and demonstration power plant with demountable magnets,” Fusion Engineering and Design 100, 378–405 (2015).
  * Buttery, RJ and Park, JM and McClenaghan, JT and Weisberg, D and Canik, J and Ferron, J and Garofalo, A and Holcomb, CT and Leuer, J and Snyder, PB et al. [2021] Buttery, RJ and Park, JM and McClenaghan, JT and Weisberg, D and Canik, J and Ferron, J and Garofalo, A and Holcomb, CT and Leuer, J and Snyder, PB et al., “The advanced tokamak path to a compact net electric fusion pilot plant,” Nuclear Fusion 61, 046028 (2021).
  * Chan  _et al._ [2015] V. Chan, A. Costley, B. Wan, A. Garofalo, and J. Leuer, “Evaluation of CFETR as a Fusion Nuclear Science Facility using multiple system codes,” Nuclear Fusion 55, 023017 (2015).
  * Reux, Cédric and Di Gallo, Luc and Imbeaux, Frédéric and Artaud, J-F and Bernardi, P and Bucalossi, Jérôme and Ciraolo, Guido and Duchateau, J-L and Fausser, C and Galassi, Davide et al. [2015] Reux, Cédric and Di Gallo, Luc and Imbeaux, Frédéric and Artaud, J-F and Bernardi, P and Bucalossi, Jérôme and Ciraolo, Guido and Duchateau, J-L and Fausser, C and Galassi, Davide et al., “DEMO reactor design using the new modular system code SYCOMORE,” Nuclear Fusion 55, 073011 (2015).
  * Christiansen, JP and Cordey, JG and Thomsen, K and Tanga, A and DeBoo, JC and Schissel, DP and Taylor, TS and Kardaun, OJWF and Wagner, F and Ryter, F et al. [1992] Christiansen, JP and Cordey, JG and Thomsen, K and Tanga, A and DeBoo, JC and Schissel, DP and Taylor, TS and Kardaun, OJWF and Wagner, F and Ryter, F et al., “Global energy confinement H-mode database for ITER,” Nuclear Fusion 32, 291 (1992).
  * Candy  _et al._ [2009] J. Candy, C. Holland, R. Waltz, M. R. Fahey, and E. Belli, “Tokamak profile prediction using direct gyrokinetic and neoclassical simulation,” Physics of Plasmas 16, 060704 (2009).
  * Snyder  _et al._ [2009] P. Snyder, R. Groebner, A. Leonard, T. Osborne, and H. Wilson, “Development and validation of a predictive model for the pedestal height,” Physics of Plasmas 16, 056118 (2009).
  * Rodriguez-Fernandez  _et al._ [2020] P. Rodriguez-Fernandez, N. Howard, M. Greenwald, A. Creely, J. Hughes, J. Wright, C. Holland, Y. Lin, F. Sciortino, S. Team, _et al._ , “Predictions of core plasma performance for the SPARC tokamak,” Journal of Plasma Physics 86, 865860503 (2020).
  * Snyder  _et al._ [2019] P. Snyder, J. Hughes, T. Osborne, C. Paz-Soldan, W. Solomon, M. Knolker, D. Eldon, T. Evans, T. Golfinopoulos, B. Grierson, _et al._ , “High fusion performance in Super H-mode experiments on Alcator C-Mod and DIII-D,” Nuclear Fusion 59, 086017 (2019).
  * Knolker  _et al._ [2020] M. Knolker, P. Snyder, T. Evans, T. Wilks, D. Eldon, B. Grierson, A. Järvinen, X. Jian, F. Laggner, J. McClenaghan, _et al._ , “Optimizing the super H-mode pedestal to improve performance and facilitate divertor integration,” Physics of Plasmas 27, 102506 (2020).
  * Creely  _et al._ [2020] A. Creely, M. J. Greenwald, S. B. Ballinger, D. Brunner, J. Canik, J. Doody, T. Fülöp, D. Garnier, R. Granetz, T. Gray, _et al._ , “Overview of the SPARC tokamak,” Journal of Plasma Physics 86, 865860502 (2020).
  * Staebler  _et al._ [2020] G. M. Staebler, J. Candy, E. A. Belli, J. E. Kinsey, N. Bonanomi, and B. Patel, “Geometry dependence of the fluctuation intensity in gyrokinetic turbulence,” Plasma Physics and Controlled Fusion 63, 015013 (2020).
  * Candy, Belli, and Bravenec [2016] J. Candy, E. A. Belli, and R. Bravenec, “A high-accuracy Eulerian gyrokinetic solver for collisional plasmas,” Journal of Computational Physics 324, 73–93 (2016).
  * Rodriguez-Fernandez, Howard, and Candy [2022] P. Rodriguez-Fernandez, N. Howard, and J. Candy, “Nonlinear gyrokinetic predictions of SPARC burning plasma profiles enabled by surrogate modeling,” Nuclear Fusion 62, 076036 (2022).
  * Rodriguez-Fernandez  _et al._ [2022] P. Rodriguez-Fernandez, A. Creely, M. Greenwald, D. Brunner, S. Ballinger, C. Chrobak, D. Garnier, R. Granetz, Z. Hartwig, N. Howard, _et al._ , “Overview of the SPARC physics basis towards the exploration of burning-plasma regimes in high-field, compact tokamaks,” Nuclear Fusion 62, 042003 (2022).
  * Luda, T and Angioni, C and Dunne, MG and Fable, E and Kallenbach, A and Bonanomi, N and Schneider, PA and Siccinio, M and Tardini, G and EUROfusion MST1 Team et al. [2020] Luda, T and Angioni, C and Dunne, MG and Fable, E and Kallenbach, A and Bonanomi, N and Schneider, PA and Siccinio, M and Tardini, G and EUROfusion MST1 Team et al., “Integrated modeling of ASDEX Upgrade plasmas combining core, pedestal and scrape-off layer physics,” Nuclear Fusion 60, 036023 (2020).
  * Kim  _et al._ [2023] H.-T. Kim, F. Auriemma, J. Ferreira, S. Gabriellini, A. Ho, P. Huynh, K. Kirov, R. Lorenzini, M. Marin, M. Poradzinski, N. Shi, G. Staebler, Z. Stancar, G. Stankunas, V. K. Zotta, E. Belli, F. J. Casson, C. D. Challis, J. Citrin, D. van Eester, E. Fransson, D. Gallart, J. Garcia, L. Garzotti, R. Gatto, J. Hobirk, A. Kappatou, E. Lerche, A. L.-O. C. Maggi, M. Maslov, M. Nocente14, R. Sharma1, A. D. Siena13, P. Strand11, E. Tholerus, D. Yadykin11, and J. contributors, “Validation of D-T fusion power prediction capability against 2021 JET D-T experiments,” Nuclear Fusion Special Issues (2023).
  * Staebler  _et al._ [2022] G. M. Staebler, M. Knolker, P. Snyder, C. Angioni, E. Fable, T. Luda, C. Bourdelle, J. Garcia, J. Citrin, M. Marin, _et al._ , “Advances in prediction of tokamak experiments with theory-based models,” Nuclear Fusion 62, 042005 (2022).
  * Garofalo, Andrea M and Gong, Xianzu and Grierson, Brian A and Ren, Qilong and Solomon, Wayne M and Strait, Edward J and Van Zeeland, Michael A and Holcomb, Christopher T and Meneghini, Orso and Smith, Sterling P et al. [2015] Garofalo, Andrea M and Gong, Xianzu and Grierson, Brian A and Ren, Qilong and Solomon, Wayne M and Strait, Edward J and Van Zeeland, Michael A and Holcomb, Christopher T and Meneghini, Orso and Smith, Sterling P et al., “Compatibility of internal transport barrier with steady-state operation in the high bootstrap fraction regime on DIII-D,” Nuclear Fusion 55, 123025 (2015).
  * Murakami, Masanori and Park, Jin Myung and Giruzzi, G and Garcia, J and Bonoli, P and Budny, RV and Doyle, EJ and Fukuyama, A and Hayashi, N and Honda, M et al. [2011] Murakami, Masanori and Park, Jin Myung and Giruzzi, G and Garcia, J and Bonoli, P and Budny, RV and Doyle, EJ and Fukuyama, A and Hayashi, N and Honda, M et al., “Integrated modelling of steady-state scenarios and heating and current drive mixes for ITER,” Nuclear Fusion 51, 103006 (2011).
  * Verdoolaege, Geert and Kaye, Stanley M and Angioni, Clemente and Kardaun, Otto JWF and Maslov, Mikhail and Romanelli, Michele and Ryter, François and Thomsen, Knud and Contributors, JET and EUROfusion MST1 Team, the et al. [2021] Verdoolaege, Geert and Kaye, Stanley M and Angioni, Clemente and Kardaun, Otto JWF and Maslov, Mikhail and Romanelli, Michele and Ryter, François and Thomsen, Knud and Contributors, JET and EUROfusion MST1 Team, the et al., “The updated ITPA global H-mode confinement database: description and analysis,” Nuclear Fusion 61, 076006 (2021).
  * Meneghini, O and Smith, SP and Lao, LL and Izacard, O and Ren, Q and Park, JM and Candy, J and Wang, Z and Luna, CJ and Izzo, VA et al. [2015] Meneghini, O and Smith, SP and Lao, LL and Izacard, O and Ren, Q and Park, JM and Candy, J and Wang, Z and Luna, CJ and Izzo, VA et al., “Integrated modeling applications for tokamak experiments with OMFIT,” Nuclear Fusion 55, 083008 (2015).
  * Lütjens, Hinrich and Bondeson, Anders and Sauter, Olivier [1996] Lütjens, Hinrich and Bondeson, Anders and Sauter, Olivier, “The CHEASE code for toroidal MHD equilibria,” Computer physics communications 97, 219–260 (1996).
  * Snyder  _et al._ [2004] P. Snyder, H. Wilson, J. Ferron, L. Lao, A. Leonard, D. Mossessian, M. Murakami, T. Osborne, A. Turnbull, and X. Xu, “ELMs and constraints on the H-mode pedestal: peeling–ballooning stability calculation and comparison with experiment,” Nuclear fusion 44, 320 (2004).
  * Redl  _et al._ [2021] A. Redl, C. Angioni, E. Belli, O. Sauter, A. U. Team, and E. M. Team, “A new set of analytical formulae for the computation of the bootstrap current and the neoclassical conductivity in tokamaks,” Physics of Plasmas 28, 022502 (2021).
  * Meneghini, Orso and Smith, Sterling P and Snyder, Philip B and Staebler, Gary M and Candy, Jeffrey and Belli, E and Lao, L and Kostuk, Mark and Luce, T and Luda, Teobaldo et al. [2017] Meneghini, Orso and Smith, Sterling P and Snyder, Philip B and Staebler, Gary M and Candy, Jeffrey and Belli, E and Lao, L and Kostuk, Mark and Luce, T and Luda, Teobaldo et al., “Self-consistent core-pedestal transport simulations with neural network accelerated models,” Nuclear Fusion 57, 086034 (2017).
  * Stix [1972] T. H. Stix, “Heating of toroidal plasmas by neutral injection,” Plasma Physics 14, 367 (1972).
  * Goldston  _et al._ [1981] R. J. Goldston, D. McCune, H. Towner, S. Davis, R. Hawryluk, and G. Schmidt, “New techniques for calculating heat and particle source rates due to neutral beam injection in axisymmetric tokamaks,” Journal of Computational Physics 43, 61–78 (1981).
  * Kritz, AH and Hsuan, H and Goldfinger, RC and Batchelor, DB [1982] Kritz, AH and Hsuan, H and Goldfinger, RC and Batchelor, DB, “Ray tracing study of electron cyclotron heating in toroidal geometry,” in  _Heating in Toroidal Plasmas 1982_ (Elsevier, 1982) pp. 707–723.
  * Smirnov, AP and Harvey, RW [2001] Smirnov, AP and Harvey, RW, “The GENRAY ray tracing code,” CompX Report CompX-2000-01 (2001).
  * Kikuchi, Lackner, and Tran [2012] M. Kikuchi, K. Lackner, and M. Q. Tran,  _Fusion physics_ (International Atomic Energy Agency, 2012).
  * Imbeaux, Frédéric and Pinches, SD and Lister, JB and Buravand, Y and Casper, Thomas and Duval, Basil and Guillerminet, B and Hosokawa, Masanari and Houlberg, Wayne and Huynh, Philippe et al. [2015] Imbeaux, Frédéric and Pinches, SD and Lister, JB and Buravand, Y and Casper, Thomas and Duval, Basil and Guillerminet, B and Hosokawa, Masanari and Houlberg, Wayne and Huynh, Philippe et al., “Design and first applications of the ITER integrated modelling & analysis suite,” Nuclear Fusion 55, 123006 (2015).
  * Meneghini, O and Snoep, G and Lyons, BC and McClenaghan, J and Imai, CS and Grierson, B and Smith, SP and Staebler, GM and Snyder, PB and Candy, J et al. [2020] Meneghini, O and Snoep, G and Lyons, BC and McClenaghan, J and Imai, CS and Grierson, B and Smith, SP and Staebler, GM and Snyder, PB and Candy, J et al., “Neural-network accelerated coupled core-pedestal simulations with self-consistent transport of impurities and compatible with ITER IMAS,” Nuclear Fusion 61, 026006 (2020).
  * Lyons, B.C. and et al. [2022] Lyons, B.C. and et al., “Flexible, integrated modeling of tokamak stability, transport, equilibrium, and pedestal physics,” Phys. Plasmas to be submitted (2022).
  * Pfeiffer  _et al._ [1980] W. Pfeiffer, R. Davidson, R. Miller, and R. Waltz, “ONETWO: a computer code for modeling plasa transport in tokamaks,” Tech. Rep. (General Atomics, San Diego, CA (United States), 1980).
  * Neiser  _et al._ [2020] T. Neiser, O. Meneghini, S. Smith, M. Fasciana, G. Staebler, and J. Candy, “Big Data Validation of the TGLF Transport Model,” in  _APS Division of Plasma Physics Meeting Abstracts_ , APS Meeting Abstracts, Vol. 2020 (2020) p. ZO04.007.
  * Belli and Candy [2008] E. Belli and J. Candy, “Kinetic calculation of neoclassical transport including self-consistent electron and impurity dynamics,” Plasma Physics and Controlled Fusion 50, 095010 (2008).
  * Meneghini  _et al._ [2016] O. Meneghini, P. Snyder, S. Smith, J. Candy, G. Staebler, E. Belli, L. Lao, J. Park, D. Green, W. Elwasif, _et al._ , “Integrated fusion simulation with self-consistent core-pedestal coupling,” Physics of Plasmas 23, 042507 (2016).
  * Roach, CM and Walters, M and Budny, RV and Imbeaux, F and Fredian, TW and Greenwald, M and Stillerman, JA and Alexander, DA and Carlsson, J and Cary, JR et al. [2008] Roach, CM and Walters, M and Budny, RV and Imbeaux, F and Fredian, TW and Greenwald, M and Stillerman, JA and Alexander, DA and Carlsson, J and Cary, JR et al., “The 2008 public release of the international multi-tokamak confinement profile database,” Nuclear fusion 48, 125001 (2008).
  * Lao, LL and John, HE St and Peng, Q and Ferron, JR and Strait, EJ and Taylor, TS and Meyer, WH and Zhang, C and You, KI [2005] Lao, LL and John, HE St and Peng, Q and Ferron, JR and Strait, EJ and Taylor, TS and Meyer, WH and Zhang, C and You, KI, “MHD equilibrium reconstruction in the DIII-D tokamak,” Fusion science and technology 48, 968–977 (2005).
  * Komm  _et al._ [2017] M. Komm, P. Bílková, M. Aftanas, M. Berta, P. Böhm, O. Bogár, L. Frassinetti, O. Grover, P. Háček, J. Havlicek, _et al._ , “Contribution to the multi-machine pedestal scaling from the COMPASS tokamak,” Nuclear Fusion 57, 056041 (2017).
  * Sartori, R and Saibene, G and Horton, LD and Becoulet, M and Budny, R and Borba, D and Chankin, A and Conway, GD and Cordey, G and McDonald, D et al. [2004] Sartori, R and Saibene, G and Horton, LD and Becoulet, M and Budny, R and Borba, D and Chankin, A and Conway, GD and Cordey, G and McDonald, D et al., “Study of type III ELMs in JET,” Plasma Physics and Controlled Fusion 46, 723 (2004).
  * Paz-Soldan, Carlos and the DIII-D Team. [2021] Paz-Soldan, Carlos and the DIII-D Team., “Plasma performance and operational space without ELMs in DIII-D,” Plasma Physics and Controlled Fusion 63, 083001 (2021).
  * McDonald, DC and Cordey, JG and Thomsen, K and Kardaun, OJWF and Snipes, JA and Greenwald, M and Sugiyama, L and Ryter, F and Kus, A and Stober, J et al. [2007] McDonald, DC and Cordey, JG and Thomsen, K and Kardaun, OJWF and Snipes, JA and Greenwald, M and Sugiyama, L and Ryter, F and Kus, A and Stober, J et al., “Recent progress on the development and analysis of the ITPA global H-mode confinement database,” Nuclear Fusion 47, 147 (2007).
  * Angioni  _et al._ [2023] C. Angioni, N. Bonanomi, E. Fable, P. A. Schneider, G. Tardini, T. Luda, G. M. Staebler, A. U. Team, _et al._ , “The dependence of tokamak L-mode confinement on magnetic field and plasma size, from a magnetic field scan experiment at ASDEX Upgrade to full-radius integrated modelling and fusion reactor predictions,” Nuclear Fusion 63, 056005 (2023).
  * Chen  _et al._ [2021] J. Chen, V. S. Chan, X. Jian, X. Zhang, Q. Ren, G. Li, C. Zhou, _et al._ , “Integrated modeling of CFETR hybrid scenario plasmas,” Nuclear Fusion 61, 046002 (2021).
  * Frank  _et al._ [2022] S. Frank, C. Perks, A. Nelson, T. Qian, S. Jin, A. Cavallaro, A. Rutkowski, A. Reiman, J. Freidberg, P. Rodriguez-Fernandez, _et al._ , “Radiative pulsed L-mode operation in ARC-class reactors,” Nuclear Fusion 62, 126036 (2022).
  * Miller  _et al._ [1998] R. Miller, M. Chu, J. Greene, Y. Lin-Liu, and R. Waltz, “Noncircular, finite aspect ratio, local equilibrium model,” Physics of Plasmas 5, 973–978 (1998).
  * Federici, G and Kemp, R and Ward, D and Bachmann, C and Franke, T and Gonzalez, S and Lowry, C and Gadomska, M and Harman, J and Meszaros, B et al. [2014] Federici, G and Kemp, R and Ward, D and Bachmann, C and Franke, T and Gonzalez, S and Lowry, C and Gadomska, M and Harman, J and Meszaros, B et al., “Overview of EU DEMO design and R&D activities,” Fusion Engineering and Design 89, 882–889 (2014).
  * Slendebroek  _et al._ [2022] T. Slendebroek, J. McClenaghan, A. M. Garofalo, L. Lao, and S. Ding, “The importance of core fueling for the Q>5 inductive mission for the Burning Experimental Superconducting Tokamak (BEST),” Bulletin of the American Physical Society (2022).


# External chemical gradient leads to efficient swimming of chiral swimmers

Ruma Maity  Department of Physics, Indian Institute Science Education and
Research Bhopal, Bhopal 462066, India  P.S. Burada<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology Kharagpur, Kharagpur
721302, India

###### Abstract

External gradients can strongly influence the collective behavior of
microswimmers. In this paper, we study the behavior of two hydrodynamically
interacting self-propelled chiral swimmers, in the low-Reynolds number regime,
under the influence of an external linear chemical gradient. We use the
generalized squirmer model called the chiral squirmer, a spherically shaped
body with an asymmetric surface slip velocity, to represent the swimmer. We
find that the external gradient favors the attraction between the swimmers
and, in some situations, leads to a bounded state in which the swimmers move
in a highly synchronous manner. Further, due to this cooperative motion,
swimmers efficiently reach the chemical target compared to the individual
swimmers. This study may help to understand the collective behavior of chiral
swimmers and to design synthetic microswimmers for targeted drug delivery.

##  I Introduction

Microorganisms are tiny, but their collective behavior can influence the
climate and human life in various ways. For example, massive plankton blooms
in the ocean, harmful red tides along the coastline, and bioconvection Platt
(1961). Individual microorganisms exhibit various propulsive mechanisms
depending on their size, shape, etc. Lauga (2020). When microorganisms, e.g.,
E-Coli, sperm cells, Paramesium, and Volvox, swim in a fluid, the viscous
force from the surrounding medium dominates over the inertia of the body
Purcell (1977). The corresponding Reynolds number is nearly zero for these
microswimmers. The hydrodynamic interaction among the microswimmers may lead
to long-range rotational order Simha and Ramaswamy (2002); Ramaswamy (2017).
It is due to the existence of dynamical instability in the system. Note that
it is possible to stabilize the instability by various means, for example,
chemical signaling Nejad and Najafi (2017). Biofilm formation is an apt
example of collective sensing of chemical stimulus by a bacteria colony Hall-
Stoodley, Costerton, and Stoodley (2004).

A popular squirmer model Lighthill (1952); Blake (1971), a sphere with an
axisymmetric surface slip, is used to understand the hydrodynamic behavior of
spherical or nearly spherical-shaped microorganisms. Note that the squirmer
can exhibit translational motion only. However, in general, many
microorganisms exhibit translational and rotational motion. Resulting in
swimming in helical paths and exhibiting chiral flows Crenshaw (1993); Burada,
Maity, and Jülicher (2022). Notably, the non-axisymmetric flow feature of a
swimmer is common in literature Pak and Lauga (2014); Ghose and Adhikari
(2014); Felderhof and Jones (2016). Therefore, the generalized squirmer model
called the chiral squirmer Burada, Maity, and Jülicher (2022); Pak and Lauga
(2014); Maity and Burada (2022a) is more applicable for studying the
collective behavior of spherical self-propelled bodies. In the chiral squirmer
model, the tangential slip velocity on the body’s surface considers both the
body’s translational and rotational degree of freedom. Also, note that a pair
of chiral swimmers exhibit a peculiar bounded state, in addition to commonly
observed states such as convergence and divergence states, due to the non-
axisymmetric nature of the body as reported in our earlier works Burada,
Maity, and Jülicher (2022); Maity and Burada (2022a). When the separation
distance between the swimmers is significant compared to the swimmer’s size,
the swimmer’s flow field affects the other swimmers’ migration velocity and
vice versa. However, when the swimmers approach each other, the lubrication
forces and torques arise and control the hydrodynamic behavior of the swimmers
Wang and Ardekani (2013); Di Leonardo  _et al._ (2011); Maity and Burada
(2022a).

In general, microorganisms sense the external stimulus in their vicinity and
respond to it by moving towards Kirchman (1994) or away from it. This
directional movement in the presence of external chemical stimulus is known as
chemotaxis. It plays a vital role in biological pattern formation Saintillan
and Shelley (2008), altering the viscosity of the medium Sokolov and Aranson
(2009); Haines  _et al._ (2009), generating turbulent flows Dunkel  _et al._
(2013), etc. Chemotaxis of sperm cells, E. coli, Dictyostelium, Paramecium,
Tetrahymena thermophila, Amoeba proteus, etc. is ubiquitous in nature Larsen,
Macnab, and Koshland (1974); Jikeli  _et al._ (2015); Sarvestani, Amir, and T
(2016). Artificial swimmers are also able to exhibit chemotaxis under
controlled experimental conditions. Unlike the swimming strategies of the
living microorganisms, different mechanisms like the Marangoni effect Jin,
Krüger, and Maass (2017), the active diffusion Hong  _et al._ (2007), the
interfacial tension Paxton  _et al._ (2004), and the extraction of work from
active medium Geiseler  _et al._ (2016) prevail for artificial swimmers to
respond to the chemical gradient.

In this article, we study the response of a pair of chiral swimmers to the
external chemical gradient by considering the near- and far-field hydrodynamic
interaction between them. The strength of the gradient solely controls the
response of a single swimmer to the external gradient. However, for a pair of
swimmers, the hydrodynamic interaction between them plays a crucial role in
the chemotaxis of the swimmers. The latter plays a constructive role in the
success of chemotaxis. The paper is organized as follows. The general chiral
squirmer model is briefly discussed in section II. In section IV, we analyze
the hydrodynamic behavior of two chiral swimmers in the presence of a linear
chemical gradient. The main conclusions are provided in section VI.

##  II Hydrodynamics of a chiral squirmer

The low Reynolds number swimmers Lighthill (1952); Purcell (1977) obey the
Stokes equation Happel and Brenner (1983) given by

$\eta\nabla^{2}\mathbf{u}=\nabla p\,,$ (1)

where $\eta$ is the viscosity of the fluid around the body, $\mathbf{u}$ is
the velocity field, and $p$ is the corresponding pressure field which plays
the role of a Lagrange multiplier to impose the incompressibility constraint
$\nabla\cdot\mathbf{u}=0$. In the chiral squirmer model, the effective slip
velocity Burada, Maity, and Jülicher (2022) on the surface of a spherical non-
deformable body of radius $a$ is prescribed by

$\displaystyle\mathbf{u}_{s}=\displaystyle{\sum\limits_{l=1}^{\infty}}\,\displaystyle{\sum\limits_{m=-l}^{l}}\Big{[}-\beta_{lm}\,{\bm{\nabla}}_{s}\left(P_{l}^{m}(\cos\theta)\,e^{im\phi}\right)+\gamma_{lm}\,{\bf\hat{r}}\times{\bm{\nabla}}_{s}\left(P_{l}^{m}(\cos\theta)\,e^{im\phi}\right)\Big{]}\,,$
(2)

where $\nabla_{s}$ =
$\mathbf{e_{\theta}}{\partial}/{\partial\theta}+\mathbf{e_{\phi}}({1}/{\sin\theta})({\partial}/{\partial\phi})$
is the surface gradient operator, ${\bf\hat{r}}$ is the unit vector in radial
direction, $P_{l}^{m}(\cos\theta)\,e^{im\phi}$ are non-normalized spherical
harmonics, where $P_{l}^{m}(\cos\theta)$ are the associated Legendre
polynomials of order $m$ and degree $l$. The complex coefficients $\beta_{lm}$
and $\gamma_{lm}$ are the surface slip velocity mode amplitudes. We introduce
the real and imaginary parts of these amplitudes as
$\beta_{lm}=\beta_{lm}^{r}+i\,m\,\beta_{lm}^{i}$ and
$\gamma_{lm}=\gamma_{lm}^{r}+i\,m\,\gamma_{lm}^{i}$ with complex conjugates
$\beta_{lm}^{\ast}=(-1)^{m}\beta_{l,-m}$ and
$\gamma_{lm}^{\ast}=(-1)^{m}\gamma_{l,-m}$, respectively.

Figure 1:  Chiral squirmer (a) with a prescribed surface slip velocity in the
body-fixed reference frame ${\bf n}$, ${\bf b}$, and ${\bf t}$. The
corresponding velocity and rotation of the swimmer are $\mathbf{V}$ and
$\bm{\Omega}$, respectively. The types of swimmers are shown in (b). For the
puller-type swimmer $\beta_{20}^{r}>0$ and for the pusher-type swimmer
$\beta_{20}^{r}<0$ (see the text for details). $\beta_{20}^{r}=0$ corresponds
to the neutral swimmer (not shown).

Interestingly, the swimming velocity $\mathbf{V}$ and the rotation rate
$\bm{\Omega}$ of the swimmer can be obtained directly from the surface slip
velocity Eq. 2, using the reciprocal theorem, without explicitly calculating
the corresponding flow- and stress-field (see Ref.Stone and Samuel (1996)).
They can be obtained in the body-fixed reference frame as
$\mathbf{V}=2(\beta_{11}^{r},\,\beta_{11}^{i},\,\beta_{10}^{r})/3$ and
${\bm{\Omega}}=(\gamma_{11}^{\,r},\,\gamma_{11}^{\,i},\,\gamma_{10}^{\,r})/a$
Burada, Maity, and Jülicher (2022). Without loss of generality, the frame
${\bf n},{\bf b}$, and ${\bf t}$ can be chosen such that ${\bf t}$ points in
the direction of velocity. Accordingly, we have
$\beta_{11}^{r}=\beta_{11}^{i}=0$ and write $\beta_{10}^{r}=3v/2$ such that
$v=|\mathbf{V}|$ is the speed of the swimmer. Thus, the velocity and rotation
rate of the chiral swimmer read Burada, Maity, and Jülicher (2022),

$\displaystyle\mathbf{V}$ $\displaystyle=v\,{\bf t}\ ,$ (3)
$\displaystyle{\bm{\Omega}}$ $\displaystyle=\frac{\gamma_{11}^{\,r}}{a}\,{\bf
n}+\frac{\gamma_{11}^{\,i}}{a}\,{\bf b}+\frac{\gamma_{10}^{\,r}}{a}\,{\bf
t}\,.$ (4)

Also, for simplicity, we choose that the swimmer has the rotation rate in the
${\bf n-t}$ plane only. With this choice, we have $\gamma_{11}^{\,i}=0$. In
addition, we choose the magnitude of the rotation as $|{\bm{\Omega}}|=v/a$
such that the components of the rotation rate can be expressed as
$\gamma_{11}^{\,r}/a=(v/a)\sin\chi\,\,,\gamma_{11}^{\,i}/a=0$ and
$\gamma_{10}^{\,r}/a=(v/a)\cos\chi$, where $\chi$ is the angle between
$\mathbf{V}$ and ${\bm{\Omega}}$. Fig. 1(a) depicts the chiral squirmer with
non-axisymmetric surface slip (Eq. 2), its swimming direction, and the
rotation rate. The corresponding equations of motion of the swimmer can be
obtained using the force- and torque-balance conditions as

$\displaystyle\dot{\textbf{q}}=\mathbf{V},\,\,\,\,\dot{\textbf{n}}=\bm{\Omega}\times\textbf{n},\,\,\,\,\dot{\textbf{b}}=\bm{\Omega}\times\textbf{b},\,\,\,\,\dot{\textbf{t}}=\bm{\Omega}\times\textbf{t}\,,$
(5)

where q is the swimmer’s position, and the dot represents the derivative with
respect to time. For $\mathbf{V}\parallel\bm{\Omega}$, we get $\chi=0$, and
the resulting swimming path is a straight line. In this case, the swimmer
rotates around the axis of motion. For $\chi=\pi/2$, the chiral swimmer moves
in a circular path in a plane. For other values of $\chi$, the path of the
chiral swimmer is helix (see Refs. Burada, Maity, and Jülicher (2022); Maity
and Burada (2022a)).

Using the boundary conditions, $\mathbf{u}=\mathbf{u}_{s}$ at the surface of
the chiral squirmer and $\mathbf{u}\to 0$ as $\mathbf{r}\to\infty$, i.e., in
the laboratory frame of reference, the Stokes equation (Eq. 1) can be solved
analytically to obtain the corresponding flow- and pressure-field Burada,
Maity, and Jülicher (2022); Maity and Burada (2022a). They read,

$\displaystyle{\bf u}_{\mathrm{lf}}({\bf r})$
$\displaystyle=\frac{3v}{2}\frac{a^{3}}{r^{3}}\left[P_{1}({\bf
t}\cdot{\bf\hat{r}})\,{\bf\hat{r}}-\frac{\bf
t}{3}\right]+3\,\beta_{20}^{r}\left(\frac{a^{4}}{r^{4}}-\frac{a^{2}}{r^{2}}\right)\,P_{2}({\bf
t}\cdot{\bf\hat{r}})\,{\bf\hat{r}}$
$\displaystyle+\beta_{20}^{r}\frac{a^{4}}{r^{4}}P_{2}^{\prime}\left({\bf
t}\cdot{\bf\hat{r}}\right)[({\bf t}\cdot{\bf\hat{r}}){\bf\hat{r}}-{\bf
t}]-\gamma_{20}^{\,r}\,\frac{a^{3}}{r^{3}}\,P_{2}^{\prime}\left({\bf
t}\cdot{\bf\hat{r}}\right){\bf t}\times{\bf\hat{r}}\,,$ (6) $\displaystyle
p_{\mathrm{lf}}({\bf r})$
$\displaystyle=-2\eta\,\beta_{20}^{r}\,\frac{a^{2}}{r^{3}}\,P_{2}\left({\bf
t}\cdot{\bf\hat{r}}\right)\,,$ (7)

where ${\bf t}$ is the swimming direction, $r$ is the distance from the center
of the swimmer where the flow field is determined, ${\bf\hat{r}}={\bf r}/r$ is
the unit radial vector, $P_{2}(x)$ denotes a second-order Legendre polynomial,
and $P_{2}^{\prime}=dP_{2}/dx$ with
$x=\mathbf{t}\cdot\mathbf{\hat{r}}=\cos\theta$. The flow field in the body-
fixed frame (bf) can be obtained from that in the lab frame (lf) as ${\bf
u}_{\mathrm{bf}}({\bf r})={\bf u}_{\mathrm{lf}}({\bf
r})-\mathbf{V}-\bm{\Omega}\times{\bf r}$. Note that in Eq. (II), we have only
written terms up to $l=2$ and do not consider higher-order contributions since
they decay more rapidly with $r$. Besides, to have a minimal model, we have
ignored $l=2$ modes with $m\neq 0$. However, it is straightforward to include
the additional terms in the analysis. The sign of $\beta_{20}^{r}$ decides the
nature of the swimming. For example, $\beta_{20}^{r}>0$ corresponds to puller-
type swimmer and $\beta_{20}^{r}<0$ corresponds to pusher-type swimmer (see
Fig. 1(b)). While pullers have an extensile force dipole, resulting, e.g.,
from the front part of the body, pushers have a contractile force dipole
stemming, e.g., from the rear part of the body Lauga (2020).

##  III Hydrodynamic interaction in the presence of lubrication forces

At low Reynolds numbers, microswimmers strongly influence their fluidic
environment, as do their nearby swimmers. When the swimmers approach each
other, the lubrication forces and torques arise, which play a vital role in
the hydrodynamic behavior of the swimmers Wang and Ardekani (2013); Yoshinaga
and Liverpool (2017); Maity and Burada (2022a). It has been reported recently
that a pair of chiral squirmers portray various behaviors, e.g., divergence,
convergence, and even a bounded state Burada, Maity, and Jülicher (2022) as a
result of their mutual hydrodynamic interaction. However, considering the
lubrication forces and torques (as in the dense suspension of swimmers), the
swimmers do not exhibit convergence behavior Maity and Burada (2022a). In this
work, we include the lubrication effects. The lubrication force acting on a
swimmer can be calculated by knowing the velocity field of the nearby swimmer
Wang and Ardekani (2013); Yoshinaga and Liverpool (2017). Due to this
lubrication force, the corresponding velocity contribution to a chiral swimmer
from the other swimmer, and vice versa, is Maity and Burada (2022a)

$\displaystyle U=-2a^{2}B\epsilon\ln\epsilon\,,$ (8)

where $\epsilon$ is half of the separation distance between the swimmers,
$B=\beta_{10}^{r(1)}t_{13}-\beta_{10}^{r(2)}t_{23}$, $\beta_{10}^{r(1)}$ and
$\beta_{10}^{r(2)}$ are $l=1$ and $m=0$ modes (see the text after Eq. 2) of
swimmer one and two, respectively,
$t_{13}=\mathbf{t}_{1}\cdot\mathbf{e}_{{}_{Z}}$,
$t_{23}=\mathbf{t}_{2}\cdot\mathbf{e}_{{}_{Z}}$, $\mathbf{e}_{{}_{Z}}$ is the
unit vector along the $Z$ direction, and $\mathbf{t}_{1}$ and $\mathbf{t}_{2}$
are the corresponding orientations of the swimmers (see Ref. Maity and Burada
(2022a) for more details). Note that the lubrication torques are of the order
$O(\epsilon^{1/2})$, and which can be neglected in the limit $\epsilon\to 0$.

The equations of motion of a swimmer in the presence of another swimmer, by
considering lubrication forces into account, read as

$\displaystyle{\bf\dot{q}}_{i}$
$\displaystyle=\mathbf{V}_{i}+\mathbf{U}^{\mathrm{lub}}_{i}+\sum^{2}\limits_{\begin{subarray}{c}j=1\,;\,i\neq
j\end{subarray}}\mathbf{u}_{j}(\mathbf{q}_{ij},\mathbf{n}_{j},\mathbf{b}_{j},\mathbf{t}_{j})\,$
$\displaystyle\left[\begin{array}[]{c}{\bf\dot{n}}_{i}\\\ {\bf\dot{b}}_{i}\\\
{\bf\dot{t}}_{i}\end{array}\right]$
$\displaystyle=\left[\bm{\Omega}_{i}+\sum^{2}\limits_{\begin{subarray}{c}j=1\,;\,i\neq
j\end{subarray}}\bm{\omega}_{j}(\mathbf{q}_{ij},\mathbf{n}_{j},\mathbf{b}_{j},\mathbf{t}_{j})\right]\times\left[\begin{array}[]{c}{\bf
n}_{i}\\\ {\bf b}_{i}\\\ {\bf t}_{i}\end{array}\right]\,,$ (15)

where
$\mathbf{U}^{\mathrm{lub}}_{i}=U(\cos\theta^{\prime}\mathbf{t}_{i}-\sin\theta^{\prime}\mathbf{n}_{i})$
is the additional velocity contribution arising due to the other swimmer in
the lubrication region,
$\theta^{\prime}=\cos^{-1}(\mathbf{t}_{i}\cdot\mathbf{e}_{{}_{Z}})$, and the
vorticity field $\bm{\omega}=({\bm{\nabla}}\times{\bf u})/2$. Note that
$\mathbf{U}^{\mathrm{lub}}_{i}=0$ for $R>2(a+\epsilon)$, and
$\mathbf{U}_{i}+\sum^{2}\limits_{\begin{subarray}{c}j=1\,;\,i\neq
j\end{subarray}}\mathbf{u}_{j}=\sum^{2}\limits_{\begin{subarray}{c}j=1\,;\,i\neq
j\end{subarray}}\bm{\omega}_{j}=0$ for $R\leq 2(a+\epsilon)$, where
$\epsilon\ll a$ and $R=|\textbf{q}_{ij}|$ is the radial distance between the
swimmers.

##  IV Influence of external linear chemical gradient

Consider that the two swimmers are subjected to an external chemical
concentration field $c(\mathbf{q})$; see the schematic diagram Fig. 2. As the
swimmers sense the gradient, they move from a lower- to a higher-concentration
region. Consequently, the local concentration field that the swimmers
experience changes with time. The stimulus level ($s=c(\mathbf{q}(t))$) on the
swimmer’s body is thus a function of time. A swimmer’s response to the
stimulus is related to the relative change in the chemoattractant
concentration ($\Delta c/c$) in the vicinity Sourjik and Wingreen (2012). The
binding of the chemoattractant to the body’s receptors, known as activation,
and the unbinding of that, known as deactivation, triggers the internal
chemotactic signaling system of the body. In our model systems, it results in
the modification of the slip coefficients (Eq. 2) of the swimmer, which leads
change in the velocity and rotation rate of the swimmer Friedrich (2009);
Böhmer  _et al._ (2005); Maity and Burada (2019, 2022b).

Figure 2:  Schematic diagram of two hydrodynamically interacting chiral
swimmers in the presence of a linear chemical gradient applied only along the
$x$ direction. The separation distance between the swimmers is $R$.

To study the chemotaxis of two swimmers, we use the Barkai–Leibler model
Barkai and Leibler (1997); Friedrich (2009); Friedrich and Jülicher (2007) for
the adaptation and relaxation mechanism of the swimmer to the external
chemical gradient. Note that this model is for information processing and
memory via the internal biochemical network, and it has been used widely for
bacterial and sperm chemotaxis Böhmer  _et al._ (2005); Friedrich (2009). The
ability of the system to adapt and respond to the external stimulus can be
captured by a simple dynamical system,

$\displaystyle\begin{split}\sigma\dot{a_{b}}&=p_{b}(s_{b}+s)-a_{b}\,,\end{split}$
(16a) $\displaystyle\begin{split}\mu\dot{p_{b}}&=p_{b}(1-a_{b})\,,\end{split}$
(16b)

where $\sigma$ is the relaxation time, $\mu$ is the adaptation time,
$a_{b}(t)$ is the dimensionless output variable, $p_{b}(t)$ is the dynamic
sensitivity related to adaptation, and $s_{b}(t)$ arises from the background
activity of receptors in the absence of the chemical stimulus. The
chemoattractant has a dimension of concentration. Note that these equations
are valid for weak concentration gradients only. Under a constant stimulus
$s(t)=S_{c}$, the system reaches a steady state for which $a_{b}=1$ and
$p_{b}=1/(s_{b}+S_{c})$. Therefore, the dynamic sensitivity ($p_{b}$)
maintains an inverse relation with $s_{b}$ and $s$. With increasing stimulus
level, $p_{b}$ decreases. The system here is totally adaptive as $a_{b}(t)$ is
independent of $S_{c}$.

In the presence of an external stimulus, the slip coefficients (Eq. 2) of the
swimmer are modified as,

$\displaystyle X=X^{(0)}+X^{(1)}\left(a_{b}(t)-1\right)\,,$ (17)

where $X=\left\\{{\beta_{lm},\gamma_{lm}}\right\\}$, the unperturbed slip
coefficients are
$X^{(0)}=\left\\{{\beta_{lm}^{(0)},\gamma_{lm}^{(0)}}\right\\}$, and the
perturbation due to the external gradient is
$X^{(1)}=\left\\{{\beta_{lm}^{(1)},\gamma_{lm}^{(1)}}\right\\}$. Here, we
consider a linear chemical gradient given by Friedrich (2009),

$c(\mathbf{q})=c_{0}+\mathbf{c_{1}}\cdot\mathbf{q}\,,$ (18)

where $\mathbf{c_{1}}=\nabla c(\mathbf{q})=c_{1}\mathbf{i}$ is the chemical
gradient applied in the $x$ direction and
$\mathbf{q}=\mathbf{i}x+\mathbf{j}y+\mathbf{k}z$ is the position vector. The
chemotactic stimulus is $S(t)=c(\mathbf{q(t)})$ Friedrich (2009).

Using Eq. (III), we numerically calculate the trajectories of a pair of chiral
swimmers, which are hydrodynamically interacting, in the presence of a linear
chemical gradient. Note that the velocity and the rotation rate of the chiral
squirmer are functions of the slip coefficients, which are affected by the
chemical gradient; see Eq. 17. As a result, the hydrodynamic interaction
between the swimmers is affected. In the present work, we consider chiral
swimmers having velocities of equal magnitudes, i.e.,
$|\mathbf{V}_{1}|=|\mathbf{V}_{2}|=v$. However, the rotation rates of the
swimmers are in general different and read,
$\bm{\Omega_{1}}=v(\sin\chi_{1},0,\cos\chi_{1})/a$ for swimmer one and
$\bm{\Omega_{2}}=v(\sin\chi_{2},0,\cos\chi_{2})/a$ for swimmer two. Chemical
and hydrodynamic perturbations in the velocity and rotation rate of the
swimmers change the corresponding torsion and curvature of the swimmers’
helical trajectories. Hence their movements are altered.

As discussed, a swimmer’s flow field influences the neighboring swimmer’s
motion. The hydrodynamic flow field of a swimmer (Eq. (II)) comprises $l=1,2$
modes as leading order terms. The $l=2$ mode is important as it plays a vital
role in determining the hydrodynamic interaction between the swimmers. The
slip coefficients corresponding to $l=2$ mode are assumed to have the
following form: $3\beta_{20}^{\,r}=3\gamma_{20}^{\,r}=\lambda_{1}$ for
squirmer one, and $3\beta_{20}^{\,r}=3\gamma_{20}^{\,r}=\lambda_{2}$ for
squirmer two. The former assumptions are made to minimize the parameters of
the problem. Note that for $\lambda_{1}\neq\lambda_{2}$, the flow fields of
the swimmers are different. Thus, variation in $\chi_{i}$ (angle between
$\mathbf{V}_{i}$ and $\bm{\Omega}_{i}$) and $\pm\lambda_{i}\,(i=1,2)$ (the
strength of the flow field), quantitatively measure the nature of the
interaction between the chiral squirmers and gives rise to several interesting
swimming characteristics. Notably, different initial configurations also
influence the interaction between the swimmers. We found that, out of various
possible initial configurations for the swimmers, swimmers interact with each
other for an extended period only in the planar configuration Burada, Maity,
and Jülicher (2022) as is the case in the present study. Here, at time $t=0$,
both the swimmers are on the $xy$ plane, with an initial orientation along the
$z-$direction and an initial separation distance $R_{0}$. Our previous study
(see Ref. Burada, Maity, and Jülicher (2022)) explored the interaction between
the two chiral swimmers with the former configuration, which exhibited some
new behavior (bounded state) in addition to the ones (convergence and
divergence state) displayed by axisymmetric swimmers. We mainly focus on how
the external chemical gradient influences these observed behaviors in the
present work.

###  IV.1 State diagram in the $\chi$ \- $\lambda$ plane

Figure 3: (a)-(d) Depicting swimming states of hydrodynamically interacting
chiral squirmers performing chemotaxis, where the strength of the chemical
gradient $c_{1}=0.1$. The symbols imply the following states. Square: bounded
state (BS), closed circles: monotonic divergence state (MD), open circles:
divergence state (D), plus: locked state, and cross symbols: parallel motion.
The two swimmers initially start at $\mathbf{q}_{1}=(9,9,0)a$ and
$\mathbf{q}_{2}=(3,3,0)a$. The corresponding initial velocities and rotation
rates of the swimmers are $\mathbf{V}_{1}=\mathbf{V}_{2}=v(0,0,1)$ and
${\bm{\Omega}}_{1}={\bm{\Omega}}_{2}=v(\cos\chi,0,\sin\chi)/a$, respectively.

Fig. 3 depicts the hydrodynamic behavior of two chiral swimmers in the
presence of an external chemical gradient. When
$\lambda_{1}=\lambda_{2}=\lambda$ and $\chi_{1}=\chi_{2}=\chi$, the swimmers
are identical (see the Fig. 3 caption). The swimmers portray various behaviors
for varying $\lambda/v$ and $\chi$. As reported earlier Maity and Burada
(2022a), in the absence of a chemical gradient, two chiral swimmers exhibit
divergence, monotonic divergence, and bounded states (see Fig. 4) only due to
the lubrication effects. These observed states arise due to pure hydrodynamic
interaction between the swimmers. However, the chemical gradient converts some
monotonic divergence states into divergence states and some divergence states
into bounded states (see Fig. 4). Indeed, the chemical stimulus favors the
attraction between the swimmers. The density of bounded states is more for the
odd combination of swimmers like pusher-puller or vice-versa (see fig. 3(c)
and (d)). For other combinations, swimmers rarely exhibit bounded states (see
fig. 3(a) and (b)). The occurrence of bounded motion in even combinations of
swimmers is due to the presence of the chemical gradient. Note that in the
case of axisymmetric swimmers (without chirality), the hydrodynamic behavior
of pusher-puller and puller-pusher-type swimmers is the same. However, in the
case of chiral swimmers, the hydrodynamic behavior of different swimmers
exhibit different behaviors Burada, Maity, and Jülicher (2022). It is because
the flow pattern of a puller-type chiral swimmer is different from a pusher-
type chiral swimmer Burada, Maity, and Jülicher (2022).

Note that for $\chi=0$ and $\lambda\neq 0$, swimmers move in a straight line
in the $z$ direction, attract each other, and converge to a locked state Maity
and Burada (2022a). In this state, the swimmers approach a minimum distance
and stay together. For $\chi=0$ (straight line motion) and $\chi=\pi/2$
(circular motion), the neutral swimmers, i.e., $\lambda=0$, move parallel to
each other without influencing each other Burada, Maity, and Jülicher (2022).
Not only the chemical gradient but the initial separation distance $R_{0}$ (at
$t=0$) between the swimmers also plays a vital role in the hydrodynamic
interaction between the swimmers. Interestingly, as $R_{0}$ increases, the
bounded state is more favorable. It is reported in detail in our earlier
studies (see Refs. Burada, Maity, and Jülicher (2022); Maity and Burada
(2022a)).

We also study the interaction between non-identical swimmers, i.e.,
$\lambda_{1}\neq\lambda_{2}$ and $\chi_{1}\neq\chi_{2}$ (see Appendix.A). For
a fixed $\chi=\pi/3$ and varying flowfield strengths
$(\lambda_{1},\lambda_{2})$, we observe mainly divergence states and a few
bounded states for the different combinations of the swimmers (see fig. 6(a)).
Identical or nearly identical swimmers mostly exhibit bounded motion primarily
for the puller-pusher combination (see fig. 6(a), top right panel). However,
other combinations also exhibit bounded motion when
$\lambda_{1}=\lambda_{2}=0$ (see fig. 6(a)). The former implies that non-
identical field strength is unfavorable for the bounded motion. On the other
hand, we can fix the flowfield strengths at $\lambda=v$ and vary relative
orientations of the swimmers $(\chi_{1},\chi_{2})$ (see fig. 6(b)). Non-
identical orientations are associated with divergence and monotonic divergence
states. See Appendix A for a detailed discussion.

Figure 4:  (a)-(c) Show the trajectories of two swimmers in the presence and
in the absence of a linear chemical gradient for monotonic divergence (MD),
divergence (D), and bounded states (BS). Here, the strength of the linear
gradient $c_{1}=0.05$. The other parameters are, $\chi=\pi/6,\lambda_{1}=-2$,
and $\lambda_{2}=2$ for D state, $\chi=\pi/3,\lambda_{1}=-3,\lambda_{2}=3$ for
MD state, and $\chi=\pi/3,\lambda_{1}=1,\lambda_{2}=-1$ for BS state. (d)-(f)
The corresponding radial separation distance $R$ between the swimmers as a
function of time. Here, lengths are scaled by the radius of the swimmer $a$
and time scaled by $\tau=v/a$.

Fig. 4 depicts the trajectories and the radial separation distance between the
swimmers in the presence and absence of the chemical gradient. As mentioned
before, the simultaneous effect of the hydrodynamic interaction and the
chemical stimulus converts some of the divergence states into bounded states
(see figs. 4 (b),(e)). Here, the stimulus increases the attraction between the
swimmers. As a result, for the bounded state, the swimmers tend to align in
the direction of the chemical gradient. The oscillation in the separation
distance ($R$) between the swimmers decreases (see figs. 4(b),(c),(e),(f))
when compared to the stimulus-free situation. Whereas for the monotonic
divergence state, the swimmers diverge at a slower rate than the stimulus-free
case (see figs. 4 (a),(d)). As reported earlier, in the case of a single
chiral swimmer, the body completely aligns its helix axis to the direction of
the gradient Maity and Burada (2019). In the two-swimmer case, none of the
swimmers can make their helix axis parallel to the direction of the gradient
because of the hydrodynamic flow field generated by the neighboring swimmer.
In particular, it is due to the dominating term in the hydrodynamic flow field
$o(1/r^{2})$ and the weak chemical gradient we apply.

##  V Dependence of relative success time on chemical gradient

Figure 5:  The relative success time $\Delta\tau=(t^{0}-t^{c})/t^{0}$ as a
function of the chemical gradient strength $c_{1}$ for various observed
states, (a) for bounded states and (b) for divergence states. The diamond
symbols for the single swimmer. Here, lines connect the data points.

A question may arise whether a pair of swimmers or a single swimmer is
efficient in swimming in the direction of the external chemical gradient. To
answer this, we measure the relative success time defined as
$\Delta\tau=(t^{0}-t^{c})/t^{0}$. Here, $t^{0}$ is the time a pair of swimmers
or a single swimmer takes to reach a particular reference point on the $x$
axis in the absence of the chemical gradient, and $t^{c}$ is the same in the
presence of the chemical gradient. Note that $\Delta\tau$ increases with
$c_{1}$ and saturates for higher $c_{1}$. The latter is due to the saturation
of the internal signaling network of the bodies Friedrich and Jülicher (2007).
Henceforth, the network cannot modify its translational and rotational
velocities further.

We find pair swimming is more efficient than individual swimming (see Fig. 5).
It is because the combined effect of the hydrodynamic interaction between the
swimmers and chemical perturbation increases the motility of the swimmers.
Interestingly, the bounded swimmers are more efficient than others (see Fig.
5(a)), especially for lower values of $\chi$. However, the influence of $\chi$
on chemotactic success rate is minimal for diverging motion (see Fig. 5(b)).
Note that, the range of $c_{1}$ is from $o(10^{-6})$ to $o(10^{-1})$. For weak
gradients, divergence motion persists. With increasing strength of $c_{1}$, as
mentioned earlier, some divergence states get converted into bounded states.
Note that for values of $\chi\sim\pi/2$, swimmers effectively move in circular
paths, and correspondingly the success rate is low (not shown).

##  VI Conclusions

In this work, we have studied the combined behavior of two chiral swimmers by
considering the lubrication effect in the presence of a linear chemical
gradient. We have used a generalized squirmer model called the chiral
squirmer. In general, a pair of chiral swimmers can exhibit bounded,
divergence, convergence, and monotonic convergence states without a chemical
gradient. However, only bounded, divergence, and monotonic divergence states
are observed in the current study. The convergence and monotonic convergence
states are absent because of the lubrication force between the swimmers. Due
to this lubrication force, swimmers repel each other when they approach each
other.

Interestingly, a chemical gradient converts some divergence- and monotonic-
divergence states into bounded states. The chemical gradient favors the
attraction between the swimmers. It leads to efficient swimming. We have
presented detailed state diagrams depicting the hydrodynamic behavior of two
chiral swimmers in the presence of a chemical gradient. We have computed the
relative success time of a pair of swimmers in aligning their movement toward
the chemical gradient, whose strength is controlled by the parameter $c_{1}$.
We find that pair swimming is more efficient than individual swimming as the
combined effect of the hydrodynamic interaction between the two swimmers and
chemical perturbation increases the motility of the swimmers. Interestingly,
the bounded swimmers are more efficient than others.

This study is not limited to chemotaxis but can be generalized to various
external gradients, e.g., phototaxis and gyrotaxis. For these cases, one needs
to identify how the slip coefficients of the chiral swimmer can be modified
according to the external gradients. Also, the present study helps design
artificial micro locomotors which can be used for drug delivery and targeted
applications.

##  VII Acknowledgements

The authors acknowledge the financial support of the Indian Institute of
Technology Kharagpur, India.

##  VIII Data Availability

The data supporting this study’s findings are available from the corresponding
author upon reasonable request.

##  Appendix A State diagrams of $\lambda_{1}$ \- $\lambda_{2}$ and $\chi_{1}$
\- $\chi_{2}$

Figure 6:  State diagrams corresponding to two hydrodynamically interacting
chiral swimmers in the presence of chemotaxis. (a) For varying field-strengths
$(\lambda_{1},\lambda_{2})$ and (b) for $(\chi_{{}_{1}},\chi_{{}_{2}})$. See
the main text for details. The initial conditions and symbol descriptions are
the same as in figure 3. We set $\chi_{1}=\chi_{2}=\pi/3$ for the left panel
and $|\lambda_{1}|=|\lambda_{2}|=|\mathbf{V}|=v$ for the right panel.

Fig. 6(a) depicts the $\lambda_{1}-\lambda_{2}$ state diagram of two chiral
swimmers, in the presence of a chemical gradient, for fixed $\chi_{1}$ and
$\chi_{2}$. The additive contributions of hydrodynamic interaction and
chemical perturbations drive the swimmers away from each other. Only the
swimmers with identical or nearly identical field strengths
$(|\lambda_{1}|=|\lambda_{2}|)$ exhibit bounded motion (see Fig. 6(a)).
Moreover, identical combinations of swimmers, i.e., puller-puller or pusher-
pusher, mostly display a divergence state with a rare occurrence of bounded
states. The probability of observing a bounded state is more if the first
swimmer is a puller and the second swimmer is a pusher. However, the opposite
combination does not work while $|\lambda_{1}|\neq 0,|\lambda_{2}|\neq 0$. The
reason is the asymmetry in the swimmer’s flow field about the orientation of
the swimmer. At different positions around a chiral swimmer, the neighboring
swimmer experiences different hydrodynamic forces. As a result, the
interaction between the swimmers is different for their different relative
positions (see Ref. Burada, Maity, and Jülicher (2022) for more details).

Fig. 6(b) shows the $\chi_{1}-\chi_{2}$ state diagram of two chiral swimmers,
in the presence of a chemical gradient, for fixed $\lambda_{1}$ and
$\lambda_{2}$. The swimmers do not exhibit bounded motion unless
$\chi_{1}=\chi_{2}$. As mentioned before, similar combinations of the
swimmers, i.e., puller-puller or pusher-pusher, exhibit bounded states rarely.
However, the opposite combinations, i.e., puller-pusher or pusher-puller, do
exhibit more bounded states (see Fig. 6 (b)). In the absence of a chemical
gradient, more monotonic divergence states are observed (see Ref. Burada,
Maity, and Jülicher (2022)). However, the presence of the chemical field
increases the attraction between the swimmers and converts some of the
monotonic divergence (MD) states into divergence (D) states (see Fig. 6(b)).

## References

  * Platt (1961) J. R. Platt, Science 133, 1766 (1961).
  * Lauga (2020) E. Lauga, in  _The fluid dynamics of cell motility_ (Cambridge University Press, Cambridge, 2020).
  * Purcell (1977) E. Purcell, Am. J. Phys. 45, 3 (1977).
  * Simha and Ramaswamy (2002) R. A. Simha and S. Ramaswamy, Physical review letters 89, 058101 (2002).
  * Ramaswamy (2017) S. Ramaswamy, Journal of Statistical Mechanics: Theory and Experiment 2017, 054002 (2017).
  * Nejad and Najafi (2017) M. R. Nejad and A. Najafi, arXiv preprint arXiv:1712.06004 (2017).
  * Hall-Stoodley, Costerton, and Stoodley (2004) L. Hall-Stoodley, J. W. Costerton, and P. Stoodley, Nature reviews microbiology 2, 95 (2004).
  * Lighthill (1952) M. Lighthill, Communications on pure and applied mathematics 5, 109 (1952).
  * Blake (1971) J. Blake, Journal of Fluid Mechanics 46, 199 (1971).
  * Crenshaw (1993) H. C. Crenshaw, Bulletin of Mathematical Biology 55, 231 (1993).
  * Burada, Maity, and Jülicher (2022) P. Burada, R. Maity, and F. Jülicher, Physical Review E 105, 024603 (2022).
  * Pak and Lauga (2014) O. S. Pak and E. Lauga, J. Eng. Math. 8, 1 (2014).
  * Ghose and Adhikari (2014) S. Ghose and R. Adhikari, J. Eng. Math. 112, 118102 (2014).
  * Felderhof and Jones (2016) B. Felderhof and R. Jones, Phys. Fluids 28, 073601 (2016).
  * Maity and Burada (2022a) R. Maity and P. Burada, Physical Review E 106, 054613 (2022a).
  * Wang and Ardekani (2013) S. Wang and A. Ardekani, Physical Review E 87, 063010 (2013).
  * Di Leonardo  _et al._ (2011) R. Di Leonardo, D. Dell’Arciprete, L. Angelani, and V. Iebba, Phys. Rev. Lett. 106, 038101 (2011).
  * Kirchman (1994) D. L. Kirchman, Microbial Ecology 28, 255 (1994).
  * Saintillan and Shelley (2008) D. Saintillan and M. J. Shelley, Physical Review Letters 100, 178103 (2008).
  * Sokolov and Aranson (2009) A. Sokolov and I. S. Aranson, Physical review letters 103, 148101 (2009).
  * Haines  _et al._ (2009) B. M. Haines, A. Sokolov, I. S. Aranson, L. Berlyand, and D. A. Karpeev, Physical Review E 80, 041922 (2009).
  * Dunkel  _et al._ (2013) J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär, and R. E. Goldstein, Physical review letters 110, 228102 (2013).
  * Larsen, Macnab, and Koshland (1974) S. Larsen, R. Macnab, and D. Koshland, Nature 249, 74 (1974).
  * Jikeli  _et al._ (2015) J. Jikeli, L. Alvarez, B. Friedrich, L. Wilson, R. Pascal, R. Colin, M. Pichlo, A. Rennhack, C. Brenker, and U. Kaupp, Nat. Commun. 6, 7985 (2015).
  * Sarvestani, Amir, and T (2016) A. N. Sarvestani, S. Amir, and A. M. T, Cell Biochem Biophys. 74, 241 (2016).
  * Jin, Krüger, and Maass (2017) C. Jin, C. Krüger, and C. C. Maass, Proceedings of the National Academy of Sciences 114, 5089 (2017).
  * Hong  _et al._ (2007) Y. Hong, N. M. Blackman, N. D. Kopp, A. Sen, and D. Velegol, Physical review letters 99, 178103 (2007).
  * Paxton  _et al._ (2004) W. F. Paxton, K. C. Kistler, C. C. Olmeda, A. Sen, S. K. St. Angelo, Y. Cao, T. E. Mallouk, P. E. Lammert, and V. H. Crespi, Journal of the American Chemical Society 126, 13424 (2004).
  * Geiseler  _et al._ (2016) A. Geiseler, P. Hänggi, F. Marchesoni, C. Mulhern, and S. Savel’ev, Physical Review E 94, 012613 (2016).
  * Happel and Brenner (1983) J. Happel and H. Brenner,  _Low Reynolds number hydrodynamics_ (Springer, 1983).
  * Stone and Samuel (1996) H. A. Stone and A. D. Samuel, Physical review letters 77, 4102 (1996).
  * Yoshinaga and Liverpool (2017) N. Yoshinaga and T. B. Liverpool, Physical Review E 96, 020603 (2017).
  * Sourjik and Wingreen (2012) V. Sourjik and N. S. Wingreen, Current opinion in cell biology 24, 262 (2012).
  * Friedrich (2009) B. Friedrich, Institut Fur Theoretische Physik Fakultat Fur Mathemtik und Naturwissenschaften Technische Universitat Dresden (2009).
  * Böhmer  _et al._ (2005) M. Böhmer, Q. Van, I. Weyand, V. Hagen, M. Beyermann, M. Matsumoto, M. Hoshi, E. Hildebrand, and U. B. Kaupp, The EMBO journal 24, 2741 (2005).
  * Maity and Burada (2019) R. Maity and P. Burada, The European Physical Journal E 42, 1 (2019).
  * Maity and Burada (2022b) R. Maity and P. Burada, J. Fluid Mech. 940, A13 (2022b).
  * Barkai and Leibler (1997) N. Barkai and S. Leibler, Nature 387, 913 (1997).
  * Friedrich and Jülicher (2007) B. M. Friedrich and F. Jülicher, Proceedings of the National Academy of Sciences 104, 13256 (2007).