# Coarse-grained gyrokinetics for the critical ion temperature gradient in
stellarators

G. T. Roberg-Clark<EMAIL_ADDRESS>G. G. Plunk  P. Xanthopoulos  Max-Planck-
Institut Für Plasmaphysik, D-17491, Greifswald, Germany

###### Abstract

We present a modified gyrokinetic theory to predict the critical gradient that
determines the linear onset of the ion temperature gradient (ITG) mode in
stellarator plasmas. A coarse-graining technique is applied to the drift
curvature, entering the standard gyrokinetic equations, around local minima.
Thanks to its simplicity, this novel formalism yields an estimate for the
critical gradient with a computational cost low enough for application to
stellarator optimization. On comparing against a gyrokinetic solver, our
results show good agreement for an assortment of stellarator designs. Insight
gained here into the physics of the onset of the ITG driven instability
enables us to devise a compact configuration, similar to the Wendelstein 7-X
device, but with almost twice the ITG linear critical gradient, an improved
nonlinear critical gradient, and reduced ITG mode transport above the
nonlinear critical gradient.

Introduction.– The excitation of the ion temperature gradient (ITG) mode in
magnetically confined fusion devices leads to turbulence that is responsible
for energy losses, which deteriorate plasma confinement. For instance, it has
been argued that, during operation of the Wendelstein 7-X (W7-X) stellarator
in electron heating scenarios, the ITG mode leads to so-called “ion
temperature clamping” (Beurskens  _et al._ , 2021), thus preventing the
heating of ions in the plasma core above 2 keV. In order to lessen the
negative effects of the ITG mode, a possible strategy to follow involves the
manipulation of the density and electron profiles, as implemented in W7-X
using pellet injections (Bozhenkov  _et al._ , 2020).

The problem might also be faced at the stage of stellarator design, by shaping
the magnetic field to reduce losses from core micro-turbulence. Such
strategies can already be applied to target electron-temperature gradient and
trapped-electron mode turbulence, the former suppressed by multiple field
periods (Plunk  _et al._ , 2014) and the latter by optimization for the
maximum-J property (Proll  _et al._ , 2012). Magnetic field shaping is the
motivation behind turbulence optimization studies, which target the rate of
the ITG transport increase (“stiffness”) as a function of the ion temperature
gradient (Mynick  _et al._ , 2010; Xanthopoulos  _et al._ , 2014; Hegna  _et
al._ , 2018; Nunami  _et al._ , 2013).

To complement this effort, one can identify another strategy: to target the
threshold or “critical gradient” of the mode itself. This should be an
effective strategy for such cases where transport is stiff, as in the
aforementioned ion clamping scenario of W7-X, such that the rapid onset of
turbulence relaxes the temperature gradient back to its threshold (Biglari
_et al._ , 1989), essentially determining the temperature profile allowed in
the plasma (Baumgaertel  _et al._ , 2013). With the actual onset of turbulence
expected to differ by a positive increment from the linear ITG threshold
(Dimits  _et al._ , 2000), the linear critical gradient can be used as a
lower-bound estimate of the nonlinear critical gradient, indeed, one that is
much simpler to model and compute.

Sustained efforts over the past four decades (Terry  _et al._ , 1982; Hahm and
Tang, 1989; Dominguez and Rosenbluth, 1989; Romanelli, 1989; Biglari  _et al._
, 1989; Kadomtsev and Pogutse, 1995; Zocco  _et al._ , 2018) have yielded an
understanding of the linear onset of the instability in toroidal geometry,
forming the basis of tokamak models Jenko  _et al._ (2001); Zocco  _et al._
(2018). In this Letter, we propose a model of the toroidal ITG critical
gradient in generic stellarator geometry. This is achieved by spatially
coarse-graining the gyrokinetic equations under the assumption of scale
separation, which, as inferred from our results, appears to be valid in a
large class of stellarator optimization lines. We then apply the model
findings to a novel configuration, generated by modifying the large scale
properties of the W7-X stellarator, that demonstrates a significant increase
in the critical gradient as well as reduced turbulent transport above the
critical gradient.

Definitions.– Following (Plunk  _et al._ , 2014), we use the standard
gyrokinetic system of equations (Brizard and Hahm, 2007) to describe
electrostatic fluctuations destabilized along a thin flux tube tracing a
magnetic field line. The ballooning transform (Dewar and Glasser, 1983) and
twisted slicing representation (Roberts and Taylor, 1965) are used to separate
out the fast perpendicular (to the magnetic field) scale from the slow
parallel scale. The magnetic field representation in field following (Clebsch)
representation reads,
$\mathbf{B}=\bm{\nabla}\psi\times\bm{\nabla}\mbox{$\alpha$}$, where $\psi$ is
a flux surface label and $\alpha$ labels the magnetic field line on the
surface. The perpendicular wave vector is then expressed as
$\mathbf{k_{\perp}}=k_{\alpha}\bm{\nabla}\alpha+k_{\psi}\bm{\nabla}\psi$,
where $k_{\alpha}$ and $k_{\psi}$ are constants, so the variation of
$\mathbf{k_{\perp}}(l)$ stems from that of the geometric quantities
$\bm{\nabla}\alpha$ and $\bm{\nabla}\psi$, with $\ell$ the field-line-
following (arc length) coordinate.

We assume Boltzmann-distributed (adiabatic) electrons, thus solving for the
perturbed ion distribution $g_{i}(v_{\parallel},v_{\perp},\mbox{$\ell$},t)$,
defined to be the non-adiabatic part of $\delta f_{i}$ ($\delta
f_{i}=f_{i}-f_{i0})$ with $f_{i}$ the ion distribution function and $f_{i0}$ a
Maxwellian. The electrostatic potential is $\phi(\mathbf{\mbox{$\ell$}})$, and
$v_{\parallel}$ and $v_{\perp}$ are the particle velocities parallel and
perpendicular to the magnetic field, respectively. The integral equation for
$g$ (Plunk  _et al._ , 2014) reads

$g(\ell)=-i\sigma(\omega-\omega_{*}^{T})f_{0}\int^{\ell}_{-\sigma\infty}d\mbox{$\ell^{\prime}$}\frac{J^{{}^{\prime}}_{0}}{|v_{\parallel}^{{}^{\prime}}|}\mbox{$\varphi$}(\mbox{$\ell^{\prime}$})\text{exp}(i\sigma
M(\ell,\mbox{$\ell^{\prime}$}))$ (1)

with

$M(\ell^{\prime},\ell)=\int_{\ell^{\prime}}^{\ell}\frac{\mbox{$\omega$}-\mbox{$\omega_{d}$}(\ell^{\prime\prime})}{\mbox{$v_{\mathrm{T}}$}x_{\parallel}}d\ell^{\prime\prime}$
(2)

where $\sigma=\text{sgn}(v_{\parallel})$, $\omega$ is the mode frequency,
$\omega_{*}^{T}=(Tk_{\alpha}/q)d\ln
T/d\psi\left(v^{2}/\mbox{$v_{\mathrm{T}}$}^{2}-3/2\right)$ is the diamagnetic
frequency, and
$J^{\prime}_{0}=J_{0}(k_{\perp}(\mbox{$\ell^{\prime}$})v_{\perp}/\Omega(\mbox{$\ell^{\prime}$}))$
is the Bessel function of zeroth order. The thermal velocity is
$\mbox{$v_{\mathrm{T}}$}=\sqrt{2T/m}$, the thermal ion Larmor radius is
$\rho=\mbox{$v_{\mathrm{T}}$}/(\Omega\sqrt{2})$, $n$ and $T$ are the
background ion density and temperature, $q$ is the ion charge,
$\varphi=q\phi/T$ is the normalized electrostatic potential, $\Omega=qB/m$ is
the cyclotron frequency, with $B=|\mathbf{B}|$, and
$v_{\parallel}^{\prime}=v_{\parallel}(\mbox{$\ell^{\prime}$})$, where
$\ell^{\prime}$ is a dummy variable of integration along the field line. The
magnetic drift frequency in the low $\beta$ approximation becomes
$\mbox{$\omega_{d}$}=(1/\Omega)(\mathbf{k_{\perp}}\cdot\mathbf{b}\times\bm{\kappa})(v_{\parallel}^{2}+v_{\perp}^{2}/2)$,
with $\bm{\kappa}=\mathbf{b}\cdot\bm{\nabla}\ \mathbf{b}$ and
$\mathbf{b}=\mathbf{B}/B$. We select the magnetic field line with $\alpha=0$,
thus setting $k_{\psi}=0$. Finally, we define a radial coordinate
$r=a\sqrt{\psi/\psi_{edge}}$, with $a$ the minor radius corresponding to the
flux surface at the edge, and $\psi_{edge}$ the toroidal flux at that
location. As usual, the temperature gradient scale length is measured relative
to the minor radius, $a/L_{T}=-(a/T)dT/dr$. To study the most unstable ITG
mode conditions, we neglect certain stabilizing factors such as the density
gradient (Proll  _et al._ , 2022; Jorge and Landreman, 2021) and plasma beta
(electromagnetic effects) (Pueschel  _et al._ , 2008).

Note that Eq. (1) is equivalent to the ion gyrokinetic equation under the
assumption of ballooning boundary conditions along the field line (Connor  _et
al._ , 1980). The system is completed by quasineutrality,

$\int d^{3}{\bf v}J_{0}g=n(1+\tau)\varphi,$ (3)

with $\tau=q_{e}T_{i}/(q_{i}T_{e})$, $T_{e}$ the electron temperature, and
$q_{e}$ the electron charge.

Scale separation and coarse-grained theory.– Spatial variation due to magnetic
geometry appears in two places in the above equations, namely $k_{\perp}$ in
the argument of $J_{0}$, and the magnetic drift frequency $\omega_{d}$, from
which a purely geometric quantity can be factored out,
$\mbox{$\omega_{d}$}\propto
K_{d}=a^{2}{\mathbf{\nabla}}\alpha\cdot\mathbf{b}\times\bm{\kappa}$, here
called the “drift curvature”. We exploit the fact that the drift curvature
contains large-scale variation (e.g., due to global magnetic shear, toroidal
variation of the normal curvature) as well as short-scale variation (e.g.,
helical ripple, local shear). We denote these macro- and microscales as $L$
and $\ell_{r}\ll L$, respectively. The macroscale is assumed to be comparable
to the linear parallel correlation length,
$\mbox{$L_{\parallel}$}=\mbox{$v_{\mathrm{T}}$}/\omega$, and the global
magnetic shear length $L_{s}=|qR/\hat{s}|$, where $R$ is the major radius,
$\hat{s}=(r/q)(dq/dr)$, and $q(r)$ is the safety factor.

We introduce a coarse-graining operator $\langle\cdot\rangle$ along the
parallel direction, which can be regarded as an average over a spatial
interval that is intermediate in size to $L$ and $\ell_{r}$. We can,
therefore, decompose any spatial function as
$\varphi=\left\langle\varphi\right\rangle+\delta\varphi$, so that
$\left\langle\delta\varphi\right\rangle=0$ by definition.

The averaging procedure arises in a straightforward manner in Eq. (1): We are
here interested in localized ITG modes, i.e., modes whose amplitudes
$|\varphi|$ peak within a confining well, defined roughly by a region of drift
curvature of destabilizing sign ($\omega_{*}^{T}\mbox{$\omega_{d}$}>1$), and
decay over a distance $\sim L$. Fixing $\ell$ to be in such a region, we see
that the behavior of the gyrokinetic system, Eqs. (1)-(3), is such that $M<1$
is satisfied within the confining well (for $|\ell-\ell^{\prime}|\lesssim L$)
and $M>1$ outside, so that the contribution to $\varphi$ is negligible under
velocity integration in Eq. (3) (see also the discussion in Roberg-Clark  _et
al._ (2021)). Furthermore, we assume $J_{0}\simeq 1$ within this region as its
argument is small there. Thus, for $|\ell-\ell^{\prime}|<L$, the integral in
Eq. (2) acts to average over small scale variation in $\omega_{d}$, except for
an asymptotically small interval in $\ell^{\prime}$, where
$|\ell-\ell^{\prime}|\sim\ell_{r}$, allowing us to take $M\approx\left\langle
M\right\rangle+\delta M$, with $\delta M$ a small correction. Noting that
$M\sim 1$, we write $\exp(i\sigma M)\approx\exp(i\sigma\left\langle
M\right\rangle)(1+i\sigma\delta M)$, and thus find that we can express
$g=\left\langle g\right\rangle+\delta g$, where $\delta g\ll\left\langle
g\right\rangle$.

Applying the coarse-graining operator to the integral equation (1), introduces
the averaged quantities $\left\langle g\right\rangle$ ,
$\left\langle\varphi\right\rangle$, and
$\left\langle\mbox{$\omega_{d}$}\right\rangle$ which vary on the scale $L$.
The same is true when applied to the quasineutrality condition, Eq. (3), and
the substitution of Eq. (1) into Eq. (3) produces a single integral equation,
identical in form to full linear gyrokinetic theory, except in coarse-grained
variables. We conclude that solutions of the linear gyrokinetic equation using
a coarse-grained drift curvature should give comparable answers to that of the
original geometry, with errors on the order
$\mathcal{O}(\mbox{$\ell$}_{r}/L)$.

Figure 1: Simulation of critical gradients versus model predictions on a log-
log scale for flux tubes taken from several stellarator geometries. (a) GENE
critical gradients versus the model equation (4) prediction for a series of
stellarator shapes, indicated by various shapes and colors. The theoretical
line
$a/L_{T,\text{crit}}=2.66(a/\mbox{$R_{\mathrm{eff}}$}+8.00a/\mbox{$L_{\parallel}$})$
[Eqn. (4)] is overlaid. (b) Predicted critical gradients from the model,
including (green circles) and excluding (empty circles) the term
$\propto\mbox{$L_{\parallel}$}^{-1}$, as compared with critical gradients
obtained from GENE simulations (red diamonds), using a series of flux tubes
taken from a negative shear circular tokamak. Each value of $\hat{s}/q$ on the
horizontal axis corresponds to a distinct flux tube.

Critical gradient model.– The coarse-grained theory is now applied to
calculate the ITG critical gradient. Our starting point is the drift resonance
condition, $\omega_{*}^{T}\sim\mbox{$\omega_{d}$}$ (Biglari  _et al._ , 1989),
which describes a balance between the drive from the temperature gradient
$\omega_{*}^{T}\sim k_{\alpha}\rho v_{T}/L_{T}$ and resonant damping from
drift curvature $\mbox{$\omega_{d}$}\sim k_{\alpha}\rho
v_{T}/\mbox{$R_{\mathrm{eff}}$}$, where $R_{\mathrm{eff}}$ is an effective
radius of curvature. It expresses the simple fact that increasing the
magnitude of drift curvature will raise the threshold gradient for instability
Baumgaertel  _et al._ (2013). In the case of a circular-cross-section tokamak
with large aspect ratio, the radius of curvature can be replaced by the major
radius $R$, as was used by Jenko et al. Jenko  _et al._ (2001). For ITG modes,
ignoring shear and finite aspect ratio effects, the Jenko formula reads
$R/L_{T,\text{crit}}=2.66$, assuming equal ion and electron temperatures and
charges. For general geometry, such a formula cannot be expected to be
accurate, but we argue that, in cases where localized modes set the threshold,
and our coarse-grained theory applies, the replacement of $R$ with a modeled
$R_{\mathrm{eff}}$ should yield a good estimate.

To calculate $R_{\mathrm{eff}}$, we take a drift curvature profile along a
thin magnetic flux tube, which extends a number of poloidal turns needed to
adequately sample a given magnetic field geometry, and find each point where
the sign of the drift curvature $K_{d}$ changes. This defines a series of
local “drift wells” along the flux tube. Coarse-graining in $\ell$ is
implemented through a least-squares quadratic fit of $K_{d}(\ell)$, of the
form
$K_{d,n}(\ell)=\mbox{$R_{\mathrm{eff}}$}_{,n}^{-1}(1-(\mbox{$\ell$}-\mbox{$\ell$}_{c})^{2}/\mbox{$\ell$}_{\text{n}}^{2})$,
for the drift curvature profile within the $n$th well, where
$\mbox{$\ell$}_{c}=(\mbox{$\ell$}_{\text{min}}+\mbox{$\ell$}_{\text{max}})/2$
is the point at the center of the drift well and
$\mbox{$R_{\mathrm{eff}}$}_{,n}$ and $\mbox{$\ell$}_{\text{n}}$ are free
fitting parameters (see Fig. 2). We exclude inverted wells, i.e., those where
the sign of the drift curvature corresponds to stable “good” curvature wells.
Thus, the effective radius of curvature corresponds to the peak value of the
coarse-grained “bad” drift curvature.

Figure 2: Fitting quadratic curves (dashed red curves) to the drift curvature
profile $-K_{d}a$ versus $\mbox{$\ell$}/a$ (solid black curves) used in the
coarse-grained model for selected stellarator geometries.

We find that an estimate for the critical gradient based purely on the above
calculation of $R_{\mathrm{eff}}$ accurately describes most stellarator
configurations, except for cases where the ITG mode is significantly
stabilized by a finite parallel correlation length (Jenko  _et al._ , 2000;
Hahm and Tang, 1989). A typical estimate of this length is the distance
between regions of good and bad curvature, i.e., the parallel connection
length, but this does not account for necessity of global magnetic shear,
which causes secular growth of drift curvature along the field line. We thus
propose an alternative definition
$\mbox{$L_{\parallel}$}^{-1}\sim\mbox{$R_{\mathrm{eff}}$}^{-1}\Theta(p_{1})\Theta(\mbox{$R_{\mathrm{eff}}$}/R_{\text{good}}-1)/(\mbox{$R_{\mathrm{eff}}$}/R_{\text{good}}+1)\times
0.5(1+\tanh[{20\>p_{2}}])$, where
$L_{w}=\mbox{$\ell$}_{\text{max}}-\mbox{$\ell$}_{\text{min}}$ and
$R_{\text{good}}$ is the scale obtained by the same fitting method as for
$R_{\mathrm{eff}}$, but applied to the neighboring region of good curvature.
The quantity $p_{1}=\mbox{$R_{\mathrm{eff}}$}/L_{w}-0.20$ is calculated to
check for a sufficiently narrow drift well ($\Theta(x)=xH(x)$ with $H$ a
Heaviside step function) and $p_{2}=\mbox{$R_{\mathrm{eff}}$}/L_{s}-0.15$ must
be positive so that the mode cannot extend beyond the neighboring good
curvature well, which occurs when the growth of the drift curvature is only
transient. The term has the expected physical behavior
$1/\mbox{$L_{\parallel}$}\rightarrow 0$ as
$\mbox{$R_{\mathrm{eff}}$}/R_{\text{good}}\rightarrow 1$ (periodic limit of
Floquet modes Bhattacharjee  _et al._ (1983); Zocco  _et al._ (2018)) and
$1/\mbox{$L_{\parallel}$}\sim 1/L_{w}$ as
$\mbox{$R_{\mathrm{eff}}$}/R_{\text{good}}\rightarrow 0$ (perfect confinement
of the mode inside the region of good curvature Plunk  _et al._ (2014)).
Taking the above considerations into account, the final model with the
correction is given by

$\frac{a}{L_{T,crit}}=2.66\left(\frac{a}{\mbox{$R_{\mathrm{eff}}$}}+8.00\frac{a}{\mbox{$L_{\parallel}$}}\right).$
(4)

where $R_{\mathrm{eff}}$ and $L_{\parallel}$ are calculated from the drift
well which minimizes $a/L_{T}$. The constant is obtained for the parallel term
by rough calibration with a series of flux tubes taken from a negative shear
circular tokamak with aspect ratio $3.87$, in which $q=2.22$ is held fixed and
$\hat{s}$ is varied from $-0.20$ to $-1.20$ by scanning over the plasma
radius, see Fig. 1(b).

The method described above, in accordance with the electrostatic limit of
gyrokinetics, is applied to vacuum flux-tube geometries generated for an
assortment of stellarator designs, including HSX (Talmadge  _et al._ , 2008),
CFQS (Shimizu  _et al._ , 2022), NCSX (Zarnstorff  _et al._ , 2001), LHD
(Takeiri  _et al._ , 2017), QPS (Nelson  _et al._ , 2003), the high-mirror
configuration of Wendelstein 7-X (Pedersen  _et al._ , 2016), a quasi-
axisymmetric stellarator NAS (Plunk and Helander, 2018), and a modified high-
mirror W7-X (KJM) case that we call W7-K, described in the next section. In
Fig. 1(a) we plot the estimated critical gradients from the formula (4) along
the horizontal axis versus the critical gradient extracted from linear flux-
tube gyrokinetic calculations (Beer  _et al._ , 1995) using the GENE code
(Jenko  _et al._ , 2001). Flux tubes are constructed (Xanthopoulos  _et al._ ,
2009) with 8 poloidal turns except for the case of NCSX (4 turns) to obtain
convergence of the critical gradient, as short flux tubes tend to artificially
destabilize modes which extend far along the field line (Faber  _et al._ ,
2018), especially in low-shear geometries. The chosen flux tube geometries are
centered at the outboard midplane of the devices ($\alpha=0$, zero toroidal
angle), and radially located at $r^{2}/a^{2}=0.5$. The critical gradient is
obtained from the GENE simulations by tracking the most unstable mode and
extrapolating its growth rate to $\gamma=0$ with two data points very close to
marginality, i.e. $\gamma\simeq 10^{-3}\sqrt{T/m}/a$ or smaller for each data
point. The marginally unstable modes tend to have values of the poloidal
wavenumber that fall in the range $k_{y}\rho\simeq 0.2$ to $0.4$. The LHD
heliotron has the highest critical gradient of all cases considered, as a
result of significant drift curvature $a/\mbox{$R_{\mathrm{eff}}$}\simeq
0.45$. Interestingly, its highly favorable comparison with W7-X (a factor of
$\sim 7$ greater) is in line with the difference in nonlinear heat fluxes
obtained via gyrokinetic simulations (Warmer  _et al._ , 2021).

Several cases of the fitted drift curvature are depicted in Fig. 2 to
demonstrate how the method behaves for four different stellarator geometries.
The case of CFQS, Fig. 2(b), is an ideal application of the coarse-grained
approach, showing that small-scale ripples in the drift curvature can be
removed and revealing the underlying parabolic structure of the well. Other
cases [Fig. 2 (c)-(d)] show that the fit can strongly deviate from the exact
profile – indeed the fluctuations need not be small according to our theory.

Figure 3: Comparing the W7-X high mirror configuration with W7-K. (Upper)
Surfaces ($r/a=1$) showing $B$ in Tesla, with a single magnetic field line
(solid black line) overlaid. (Lower) Cuts of magnetic surfaces at specified
toroidal angles mapped onto the cylindrical coordinates $Z$ versus $R$. The
cuts at $r/a=1$ are plotted in green while those at $r/a=0.5$ are in red. By
construction the surface $r/a=1$ has the same shape for both configurations
aside from a 5/3 rescaling of the toroidal angle.

Increasing the critical gradient.– Applying principles gleaned from the model
equation (4) we devise two novel configurations. The first is generated by
modifying two geometric parameters of the high-mirror configuration of the
Wendelstein 7-X stellarator to improve the critical gradient. First, we
directly control the effective radius of curvature, as compared to the minor
radius, by reducing the aspect ratio of the configuration. Specifically, we
create a configuration whose outer flux surface shape is the same as the W7-X
configuration, except for a rescaling of the major radius from $5.51$ m down
to $2.00$ m. This involves only modifying the $m=n=0$ Fourier component of the
radial component of surface shape specification in cylindrical cooordinates,
i.e. $R_{00}$ in the input file to the equilibrium code VMEC (Hirshman and
Whitson, 1983). Second, we reduce the field period number of the configuration
from $5$ down to $3$. The field period reduction broadens the parallel
connection length $L_{w}$, which we find stabilizes more extended ITG modes.
As Fig. 1 shows, the resulting low sheared configuration, “Wendelstein 7
Kompakt” (W7-K), has twice the ITG critical gradient of the high-mirror W7-X
configuration. The primary reason is a near doubling of the drift curvature.
Figure 3 shows the shape of the outermost flux surfaces for W7-X and W7-K,
colored according to $B$, as well as surface cuts at evenly spaced toroidal
angles. Furthermore, W7-K has values of $\epsilon_{\text{eff}}$ below $0.01$
out to the surface with $r^{2}/a^{2}=0.5$. This configuration thus
demonstrates that the ITG mode threshold for a given stellarator configuration
can be increased without spoiling neoclassical confinement. It also contains a
reduction in the average value of $|\bm{\nabla}s|^{2}$, implying smaller
linear growth rates for the ITG mode above marginality (Xanthopoulos  _et al._
, 2014).

The benefits of W7-K, however, are not limited to linear physics effects when
considering ITG mode transport. We find two significant improvements in W7-K
over W7-X via nonlinear GENE calculations of the time-averaged ion heat flux
$Q_{i}$, shown in Fig. 4. Firstly, the nonlinear threshold for instability
$a/L_{T,\text{crit}}$, inferred by extrapolation, is larger for W7-K, so the
strategy of improving the linear critical gradient can lead to increases in
its nonlinear counterpart. Secondly, the heat fluxes above the nonlinear
threshold are consistently lower in the case of W7-K – at a gradient
$a/L_{T,i}=2$, the heat flux is nearly a factor of two smaller. Thus a
reduction in transport as well as an increase in the nonlinear threshold has
been simultaneously achieved. We attribute both improvements to the enhanced
nonlinear zonal flow drive present in W7-K, which is outside the scope of the
present discussion but will be explored in a future publication.

Figure 4: Nonlinear heat fluxes computed by the GENE code (adiabatic
electrons, zero density gradient, $T_{e}=T_{i}$) for flux tube simulations at
$\alpha=0$ (bean cross section, toroidal angle $0$) comparing the high mirror
W7-X configuration to W7-K. The time-averaged heat steady-state heat flux
$Q_{i}/Q_{GB}$ is plotted on a log-linear scale versus $a/L_{T,i}$, where
$Q_{gB}=(\rho^{2}_{s}/a^{2})c_{s}P_{i}$ is the gyro-Bohm heat flux,
$c_{s}=\sqrt{T_{e}/m_{i}}$, $\rho_{s}=c_{s}/\Omega_{i}$, and $P_{i}=nT$.

A second “Nearly Axi-Symmetric” (NAS) configuration is generated by using the
method described in Plunk and Helander (2018), yielding a weakly non-
axisymmetric vacuum field that satisfies quasi-axisymmetry (QA) approximately,
with an aspect ratio of $3$. As shown in Fig. 1(a), NAS enjoys the largest
critical gradient of the QA configurations studied here. This is due to the
relatively large and uniform curvature on the outboard side of the device, in
contrast to NCSX, whose shaping actually leads to relatively weak curvature,
which is however compensated by strong negative shear.

Conclusions.– A coarse-graining gyrokinetic model has been implemented to
estimate the critical gradient for the ITG mode in a large variety of
stellarator configurations. The comparison against calculations with the GENE
code confirms the robustness of this approach, even for configurations whose
geometry is not ideally represented by the “smoothed” drift curvature.

We note that the coarse-graining method becomes less accurate when scale
separation is absent, as in the case of QPS [Fig. 2(d)], or when global shear
becomes extremely small. The latter limit is not approached in the cases
analyzed here, as we have not taken flux tubes at inner radii
$r^{2}/a^{2}<0.5$. Our investigation of the low-shear regime, to be published
later, confirms that the critical gradient is set by Floquet-like, extended
modes (Bhattacharjee  _et al._ , 1983; Zocco  _et al._ , 2018; Faber  _et al._
, 2018) in such cases.

The main benefit of this novel technique is twofold: i) A swift prediction of
the critical gradient without the need for tedious gyrokinetic numerical
calculations, and ii) the potential to construct new configurations that
combine neoclassical optimization with a high ITG critical gradient. Two such
configurations found here are a compact version of the Wendelstein 7-X, and a
low aspect ratio quasi-axisymmetric configuration. Our results suggest that,
using the ITG critical gradient as a figure of merit, a smaller aspect ratio
can be beneficial, but there is ample room among other degrees of freedom for
improvement in future designs of fusion reactors of the stellarator type.

The authors thank P. Helander, A. Zocco, and M. Landreman for helpful
conversations. This work was supported by a grant from the Simons Foundation
(No. 560651, G. T. R.-C.). Computing resources at the Cobra cluster at IPP
Garching and the Marconi Cluster were used to perform the simulations. This
work has been carried out within the framework of the EUROfusion Consortium,
funded by the European Union via the Euratom Research and Training Programme
(Grant Agreement No 101052200 — EUROfusion). Views and opinions expressed are
however those of the author(s) only and do not necessarily reflect those of
the European Union or the European Commission. Neither the European Union nor
the European Commission can be held responsible for them.

## References

  * Beurskens  _et al._ (2021) M. N.A. Beurskens _et al._ , “Ion temperature clamping in Wendelstein 7-X electron cyclotron heated plasmas,” Nuclear Fusion 61, 116072 (2021).
  * Bozhenkov  _et al._ (2020) S. A. Bozhenkov _et al._ , “High-performance plasmas after pellet injections in Wendelstein 7-X,” Nuclear Fusion 60, 066011 (2020).
  * Plunk  _et al._ (2014) G. G. Plunk, P. Helander, P. Xanthopoulos, and J. W. Connor, “Collisionless microinstabilities in stellarators. III. the ion-temperature-gradient mode,” Physics of Plasmas 21 (2014), 10.1063/1.4868412.
  * Proll  _et al._ (2012) J. H.E. Proll, P. Helander, J. W. Connor, and G. G. Plunk, “Resilience of quasi-isodynamic stellarators against trapped-particle instabilities,” Physical Review Letters 108, 7–10 (2012).
  * Mynick  _et al._ (2010) H. E. Mynick, N. Pomphrey, and P. Xanthopoulos, “Optimizing stellarators for turbulent transport,” Physical Review Letters 105, 1–4 (2010).
  * Xanthopoulos  _et al._ (2014) P. Xanthopoulos, H. E. Mynick, P. Helander, Y. Turkin, G. G. Plunk, F. Jenko, T. Görler, D. Told, T. Bird, and J. H.E. Proll, “Controlling turbulence in present and future stellarators,” Physical Review Letters 113, 1–4 (2014).
  * Hegna  _et al._ (2018) C. C. Hegna, P. W. Terry, and B. J. Faber, “Theory of ITG turbulent saturation in stellarators: Identifying mechanisms to reduce turbulent transport,” Physics of Plasmas 25 (2018), 10.1063/1.5018198.
  * Nunami  _et al._ (2013) M. Nunami, T. H. Watanabe, and H. Sugama, “A reduced model for ion temperature gradient turbulent transport in helical plasmas,” Physics of Plasmas 20 (2013), 10.1063/1.4822337.
  * Biglari  _et al._ (1989) H. Biglari, P. H. Diamond, and M. N. Rosenbluth, “Toroidal ion-pressure-gradient-driven drift instabilities and transport revisited,” Physics of Fluids B 1, 109–118 (1989).
  * Baumgaertel  _et al._ (2013) J. A. Baumgaertel, G. W. Hammett, and D. R. Mikkelsen, “Comparing linear ion-temperature-gradient-driven mode stability of the National Compact Stellarator Experiment and a shaped tokamak,” Physics of Plasmas 20 (2013), 10.1063/1.4791657.
  * Dimits  _et al._ (2000) A. M. Dimits _et al._ , “Comparisons and physics basis of tokamak transport models and turbulence simulations,” Physics of Plasmas 7, 969–983 (2000).
  * Terry  _et al._ (1982) P. Terry, W. Anderson, and W. Horton, “Kinetic effects on the toroidal ion pressure gradient drift mode,” Nuclear Fusion 22, 487–497 (1982).
  * Hahm and Tang (1989) T. S. Hahm and W. M. Tang, “Properties of ion temperature gradient drift instabilities in H-mode plasmas,” Physics of Fluids B 1, 1185–1192 (1989).
  * Dominguez and Rosenbluth (1989) R. R. Dominguez and M. N. Rosenbluth, “Local kinetic stability analysis of the ion temperature gradient mode,” Nuclear Fusion 29, 844–848 (1989).
  * Romanelli (1989) F. Romanelli, “Ion temperature-gradient-driven modes and anomalous ion transport in tokamaks,” Physics of Fluids B 1, 1018–1025 (1989).
  * Kadomtsev and Pogutse (1995) B. B. Kadomtsev and O. P. Pogutse,  _Reviews of Plasma Physics_, Vol. 5 (1995) p. 249.
  * Zocco  _et al._ (2018) A. Zocco, P. Xanthopoulos, H. Doerk, J. W. Connor, and P. Helander, “Threshold for the destabilisation of the ion-temperature-gradient mode in magnetically confined toroidal plasmas,” Journal of Plasma Physics 84 (2018), 10.1017/S0022377817000988.
  * Jenko  _et al._ (2001) F. Jenko, W. Dorland, and G. W. Hammett, “Critical gradient formula for toroidal electron temperature gradient modes,” Physics of Plasmas 8, 4096–4104 (2001).
  * Brizard and Hahm (2007) A. J. Brizard and T. S. Hahm, “Foundations of nonlinear gyrokinetic theory,” Reviews of Modern Physics 79, 421–468 (2007).
  * Dewar and Glasser (1983) R. L. Dewar and A. H. Glasser, “Ballooning mode spectrum in general toroidal systems,” Physics of Fluids 26, 3038–3052 (1983).
  * Roberts and Taylor (1965) K. V. Roberts and J. B. Taylor, “Gravitational resistive instability of an incompressible plasma in a sheared magnetic field,” Physics of Fluids 8, 315–322 (1965).
  * Proll  _et al._ (2022) J.H.E. Proll, G.G. Plunk, B.J. Faber, T. Görler, P. Helander, I.J. McKinney, M.J. Pueschel, H.M. Smith, and P. Xanthopoulos, “Turbulence mitigation in maximum-J stellarators with electron-density gradient,” Journal of Plasma Physics 88 (2022), 10.1017/s002237782200006x.
  * Jorge and Landreman (2021) R. Jorge and M. Landreman, “Ion-temperature-gradient stability near the magnetic axis of quasisymmetric stellarators,” Plasma Physics and Controlled Fusion 63 (2021), 10.1088/1361-6587/abfdd4, arXiv:2102.12390 .
  * Pueschel  _et al._ (2008) M. J. Pueschel, M. Kammerer, and F. Jenko, “Gyrokinetic turbulence simulations at high plasma beta,” Physics of Plasmas 15, 102310 (2008).
  * Connor  _et al._ (1980) J. W. Connor, R. J. Hastie, and J. B. Taylor, “Stability of general plasma equilibria, III,” Plasma Physics 22, 757–769 (1980).
  * Roberg-Clark  _et al._ (2021) G. T. Roberg-Clark, G. G. Plunk, and P. Xanthopoulos, “Calculating the linear critical gradient for the ion-temperature-gradient mode in magnetically confined plasmas,” Journal of Plasma Physics (2021), 10.1017/S0022377821000507.
  * Jenko  _et al._ (2000) F. Jenko, W. Dorland, M. Kotschenreuther, and B. N. Rogers, “Electron temperature gradient driven turbulence,” Physics of Plasmas 7, 1904–1910 (2000).
  * Bhattacharjee  _et al._ (1983) A. Bhattacharjee, J. E. Sedlak, P. L. Similon, M. N. Rosenbluth, and D. W. Ross, “Drift waves in a straight stellarator,” Physics of Fluids 26, 880–882 (1983).
  * Talmadge  _et al._ (2008) J. N. Talmadge, F. S. B. Anderson, D. T. Anderson, C. Deng, W. Guttenfelder, K. M. Likin, J. Lore, J. C. Schmitt, and K. Zhai, “Experimental Tests of Quasisymmetry in HSX,” Plasma and Fusion Research 3, S1002–S1002 (2008).
  * Shimizu  _et al._ (2022) A. Shimizu _et al._ , “Recent developments in engineering design for the quasi-axisymmetric stellarator CFQS,” Nuclear Fusion 62, 016010 (2022).
  * Zarnstorff  _et al._ (2001) M. C. Zarnstorff _et al._ , “Physics of the compact advanced stellarator NCSX,” Plasma Physics and Controlled Fusion 43 (2001), 10.1088/0741-3335/43/12A/318.
  * Takeiri  _et al._ (2017) Y Takeiri _et al._ , “Extension of Operational Regime of LHD towards Deuterium Experiment,” Nuclear Fusion 57, 102023 (2017).
  * Nelson  _et al._ (2003) B. E. Nelson _et al._ , “Design of the quasi-poloidal stellarator experiment (QPS),” Fusion Engineering and Design 66-68, 205–210 (2003).
  * Pedersen  _et al._ (2016) T. Sunn Pedersen _et al._ , “Confirmation of the topology of the Wendelstein 7-X magnetic field to better than 1:100,000,” Nature Communications 7, 13493 (2016).
  * Plunk and Helander (2018) G. G. Plunk and Per Helander, “Quasi-axisymmetric magnetic fields: Weakly non-axisymmetric case in a vacuum,” Journal of Plasma Physics 84, 1–18 (2018).
  * Beer  _et al._ (1995) M. A. Beer, S. C. Cowley, and G. W. Hammett, “Field-aligned coordinates for nonlinear simulations of tokamak turbulence,” Physics of Plasmas 2, 2687–2700 (1995).
  * Xanthopoulos  _et al._ (2009) P. Xanthopoulos, W. A. Cooper, F. Jenko, Yu Turkin, A. Runov, and J. Geiger, “A geometry interface for gyrokinetic microturbulence investigations in toroidal configurations,” Physics of Plasmas 16 (2009), 10.1063/1.3187907.
  * Faber  _et al._ (2018) B. J. Faber, M. J. Pueschel, P. W. Terry, C. C. Hegna, and J. E. Roman, “Stellarator microinstabilities and turbulence at low magnetic shear,” Journal of Plasma Physics 84, 1–28 (2018).
  * Warmer  _et al._ (2021) Felix Warmer _et al._ , “Impact of Magnetic Field Configuration on Heat Transport in Stellarators and Heliotrons,” Physical Review Letters 127, 225001 (2021).
  * Hirshman and Whitson (1983) S. P. Hirshman and J. C. Whitson, “Steepest-descent moment method for three-dimensional magnetohydrodynamic equilibria,” Physics of Fluids 26, 3553–3568 (1983).


Exploration of different parameter optimization algorithms within the context
of ACTS software framework

Rocky Bala Garg1, Heather Gray2, Elyssa Hofgard1, Lauren Tompkins1

1Department of Physics, Stanford University, USA
2Department of Physics, University of California at Berkeley and Lawerence
Berkeley National Laboratory, USA

> ABSTRACT
>
>
>
>
> Particle track reconstruction, in which the trajectories of charged
> particles are determined, is a critical and time consuming component of the
> full event reconstruction chain. The underlying software is complex and
> consists of a number of mathematically intense algorithms, each dealing with
> a particular tracking sub-process. These algorithms have many input
> parameters that need to be supplied in advance. However, it is difficult to
> determine the configuration of these parameters that produces the best
> performance. Currently, the input parameter values are decided on the basis
> of prior experience or by the use of brute force techniques. A parameter
> optimization approach that is able to automatically tune these parameters
> for high performance is greatly desirable. In the current work, we explore
> various machine learning based optimization methods to devise a suitable
> technique to optimize parameters in the complex track reconstruction
> environment. These methods are evaluated on the basis of a metric that
> targets high efficiency while keeping the duplicate and fake rates small. We
> focus on derivative free optimization approaches that can be applied to
> problems involving non-differentiable loss functions. For our studies, we
> consider the tracking algorithms defined within A Common Tracking Software
> (ACTS) framework. We test our methods using simulated data from ACTS
> software corresponding to the ACTS Generic detector and the ATLAS ITk
> detector geometries.

##  1 Introduction

The reconstruction of charged particles’ trajectories, commonly known as
tracking, is one of the most important but computationally extensive tasks of
the full event reconstruction chain. Precise and efficient measurements of
charged particle trajectories are critical for the entire physics program:
from reconstructing low momentum tracks in order to identify primary vertices,
through using medium momentum tracks to identify heavy flavor decays and
calculate jet energies, to accurately measuring high momentum messengers of
electroweak particles or new physics. The underlying tracking software
consists of a number of mathematically complex algorithms, each dealing with a
particular tracking sub-process. These algorithms depend on a number of input
parameters whose values need to be determined. Great efforts are put into
optimizing these input parameters as these efforts yields both resource
savings and physics performance improvements. However, it is difficult to know
the configuration of these parameters that can produce the best performance.
The values of these input parameters are highly dependent on the underlying
tracking geometry and material configuration. They are also affected by other
factors such as the number of simultaneous pp collisions known as pile-up ,
center-of-mass energy of collision, target physics process, etc. A different
set of these parameter values is needed if underlying tracking geometry or
material configuration changes. Currently, most of these optimizations are
performed using the physicist’s previous experience or specialized brute-force
techniques. Approaches that can automatically tune these parameters for high
performance are greatly desirable. Hyperparameter tuning algorithms used in
machine learning offer a number of different approaches in which the optimal
parameters can be learned. These approaches have the potential to reduce the
effort spent in finding the configurations, as well as they can scan a larger
range of parameters, allowing for potentially more efficient solutions. With
the upcoming High Luminosity Large Hadron Collider (HL-LHC) upgrade, we expect
more complex tracking environment which would require more refined parameter
tuning.

In this study, we have explored five different optimization approaches, which
include: Evolutionary Algorithms [1], Optuna Hyperparameter optimization [2],
Lipschitz Optimization [3], Scikit Optimization [4] and Orion Hyperparameter
tuning [5]. These approaches are described in more detail later. Initial
studies on this topic were presented in [6]. This work expands on that initial
work, include a systematic study of several optimization algorithms and two
components of charged particle reconstruction.

##  2 Problem Summary

###  2.1 Testing Framework

We have used “A Common Tracking Software (ACTS)” [7] for our studies which is
a community-driven tracking software suite consisting of high-level track
reconstruction modules. ACTS provides track reconstruction algorithms within a
generic, experiment-independent open-source software toolkit. It includes data
structures and algorithms for performing track reconstruction in addition to a
tool for fast track simulation. The algorithms are designed to be inherently
thread-safe to support parallel code execution. The implementation is designed
to be fully agnostic to detector technologies, and the event processing
framework, so that it can be used with different experiments. A number of
detector geometries has already been implemented within ACTS. Some of these
are: the Belle-II Detector [8], the ATLAS ITk geometry [9], the sPHENIX
Detector [10] and the LDMX Detector [11]. The underlying tracking algorithms
need to be tuned as per the specific detector geometry in order to achieve
ultimate physics performance. Therefore, ACTS provides ideal environment for
our studies and allows us to tune same parameters for different detector
geometries.

ACTS provides Generic detector geometry which is equivalent to a typical full
silicon LHC detector. The Generic Detector is a hermetic detector with silicon
based tracker consisting of high granularity pixel layers surrounded by outer
tracking layers. A schematic of Generic detector has been shown in Figure 1.
In order to perform our studies, we have used 14TeV $t\overline{t}$ dataset
with pile-up = 140 and 200. The data is generated using Pythia8 [12] and
simulated through Generic detector using the ACTS fast simulation algorithm.
The input particles are required to have $P_{T}$ $>$ 1GeV and $|\eta|$ $<$
2.5.

Figure 1: A schematic of ACTS Generic detector used in this study. This figure
is taken from [13].

###  2.2 Tracking Algorithms

Two different tracking algorithms from the ACTS tracking framework were
considered for these studies:

  * •

Track finding/fitting: The combinatorial kalman filter (CKF) [14] algorithm
implemented within ACTS framework is used. CKF combines both track finding and
fitting in a tree-search based algorithm. In ACTS, the main bottleneck in CKF
performance is the track seeding algorithm, the algorithm that provides
initial track seeds to CKF. There are around 20 parameters in track seeding
algorithm. Many of these parameters are dependent on the underlying detector
geometry. However, other geometry independent parameters can significantly
impact track seeding performance. We have chosen to optimize the following
geometry independent parameters:

    1. 1.

maxPtScattering: Upper $P_{T}$ limit for scattering angle***the angle by which
particle scatters after colliding with detector material. calculations.

    2. 2.

impactMax: maximum value for impact parameter.

    3. 3.

deltaRMin: minimum distance in $r$ between two measurements within one seed.

    4. 4.

deltaRMax: maximum distance in $r$ between two measurements within one seed.

    5. 5.

sigmaScattering: number of sigma used for scattering angle calculations.

    6. 6.

radLengthPerSeed: average radiation lengths of material on the length of a
seed.

    7. 7.

maxSeedsPerSpM: number of 3-D space-points in top and bottom layers considered
for compatibility with middle space-point.

    8. 8.

cotThetaMax: maximum cotTheta angle between two space-points in a seed to be
considered compatible.

  * •

Vertex finding/fitting: The Adaptive Multi Vertex Finder (AMVF) [15] algorithm
implemented within ACTS framework is used. AMVF performs simultaneous
vertex†††Actual proton-proton collision points. finding and multi-vertex
fitting with the help of a multi-vertex fitter. We have considered the
following 5 parameters of AMVF for optimization:

    1. 1.

tracksMaxZinterval: maximum $z$-interval used for adding tracks to multi-
vertex fit.

    2. 2.

maxVertexChi2: maximum $\chi^{2}$ value for tracks to be compatible with
fitted vertex.

    3. 3.

maxMergeVertexSignificance: maximum significance on the distance between two
vertices to allow merging.

    4. 4.

minWeight: minimum weight assigned to the track for track to be considered
compatible with vertex candidate.

    5. 5.

maximumVertexContamination: maximum vertex contamination value.

###  2.3 Optimization Algorithms Explored

We have performed optimization of CKF and AMVF parameters using 5 different
optimization algorithms. Optimized parameter performance corresponding to each
optimization algorithm has been compared with the default parameter
performance of ACTS.

Each algorithm starts with a random initial parameter configuration chosen
from the input parameter range. The algorithms run for a number of iterations
and, at each iteration, evaluate the score of the current parameter
configuration where the score/objective function is based on the tracking
algorithm performance. At each new iteration, the algorithms try to maximize
or minimize the score by providing better parameter configuration using
dedicated parameter estimation methods. The algorithms use different
estimation methods. They typically keep track of previous configurations and
scores to get better configurations at new iterations. An overview of these
optimization algorithm is provided below:

  * •

Evolutionary Algorithms (EA) [1]: This algorithm is based on a process
analogous to the genetic evolution in living organisms. It starts by
initializing a population of individuals where each individual refers to one
configuration of the input parameter set. We considered a population of size
50 with each individual assigned to same initial configuration. Then, at each
iteration, called a generation, the fitness of each individual is computed
using a customized score function based on the tracking or vertexing algorithm
performance. After the fitness evaluation, simultaneous selections and
mutations are performed on randomly chosen individuals from the population
such that better performing ones are more likely to stay for next generation,
and worse performing ones are more likely to be removed. We used a total of 16
generations for the evolution of our parameters. This algorithm maximizes the
score function and provides the configuration with the best score in the end.

  * •

Optuna HyperParameter Optimization (Optuna) [2]: Optuna is an optimization
framework that consists of multiple hyper-parameter sampling algorithms for
parameter selection at each trial. For our studies, we have used Tree-
structured Parzen Estimator (TPE) which fits one Gaussian Mixture Model (GMM)
$l(x)$ to the set of parameter values associated with the best objective
values, and a second GMM $g(x)$ to the remaining parameter values. It chooses
the parameter value $x$ that maximizes the quasi-likelihood ratio
$\frac{l(x)}{g(x)}$.

  * •

Lipschitz Optimization (LIPO) [3]: Lipschitz optimization is based on a simple
parameter-free mathematical model for finding the best $x$ that maximizes the
score function $f(x)$. The key idea is to maintain a piecewise linear upper
bound of $f(x)$ and use that to decide which $x$ to evaluate at the next step.
It defines the upper bound of the function $f(x)$ as $U(x)$ = mini
$(f(x_{i})+k.||x-x_{i}||_{2})$ such that $U(x)$ $\geq$ $f(x)$ $\forall$ $x$,
$k$ is the Lipschitz constant. This algorithm picks new points at random and
checks if the upper bound at the new point is better than the best point seen
so far, and if so, it selects this point as the next step to evaluate.

  * •

Scikit Optimization (Skopt) [4]: This is an optimization framework built
inside Python’s sklearn library [16]. It consists of a number of different
parameter optimization techniques. For this study, we have used the
forest_optimize algorithm which performs sequential optimization using
decision trees. Skopt runs over a number of iterations and minimizes the score
function to provide the best parameter configuration.

  * •

Or$\acute{\textrm{\i}}$on HyperParameter Tuning (Orion) [5]:
Or$\acute{\textrm{\i}}$on is an asynchronous framework for black-box function
optimization. It can be used as a command-line interface as well as python
interface. It also implements a number of different optimization methods
including random search and grid search. For our studies, we have used the
random search method based on the uniform probability distribution provided to
the search parameters.

###  2.4 Performance Evaluation: Score/Objective Function

The score or objective function is one of the most important component of
these optimization studies. Performance and outcome of the optimization
algorithms are highly dependent on the form of score function used. The score
function is constructed using the performance metrics of the underlying
tracking algorithm. Positive weights are given to quantities that we want to
increase while negative weights are given to quantities that we want to
decrease. Likewise, higher weights are given to more important quantities.

####  2.4.1 Construction of Score Function for the CKF

The important metrics determining the performance of CKF algorithm are:

  * •

Track Reconstruction Efficiency: fraction of generated particles that have
created at least 9 measurements on the traversed detectors and are matched
with reconstructed tracks with greater than 50% probability.

  * •

Fake Rate: Fraction of reconstructed tracks that are not associated with any
truth particle.

  * •

Duplicate rate: Fraction of reconstructed tracks associated with same truth
generated particle.

  * •

Run Time: Time taken in running the CKF algorithm.

Based on these performance metrics, we constructed the following score
function for measuring the CKF performance for different parameter
configurations:

$\textrm{Score/Objective
Function}=\textrm{Efficiency}-\left(\textrm{FakeRate}+\frac{\textrm{DuplicateRate}}{\textrm{k}}+\frac{\textrm{RunTime}}{\textrm{k}}\right),\textrm{k}=7$

where $k$ is tuned for optimal tradeoff between different quantities.

####  2.4.2 Construction of Score Function for the AMVF

The important metrics determining the performance of AMVF algorithm are:

  * •

Eff_total: Number of vertices reconstructed by AMVF algorithm out of total
detector accepted vertices.

  * •

Eff_clean: Number of reconstructed vertices associated to a single true
generated particle out of all vertices within the detector acceptance.

  * •

Eff_split: fraction of reconstructed vertices where more than one
reconstructed vertices are associated with same truth particle.

  * •

Eff_merge: fraction of reconstructed vertices where one reconstructed vertex
is associated to more than one truth particle.

  * •

Eff_fake: fraction of reconstructed vertices that are not associated to any
truth particle.

  * •

Resolution:

$\frac{{\Delta}{R}}{R}=\sqrt{\sum_{\textrm{vertices}}\frac{((x_{\textrm{reco}}-x_{\textrm{truth}})^{2}+(y_{\textrm{reco}}-y_{\textrm{truth}})^{2}+(z_{\textrm{reco}}-z_{\textrm{truth}})^{2})}{x_{\textrm{truth}}^{2}+y_{\textrm{truth}}^{2}+z_{\textrm{truth}}^{2}}}$

Based on these performance metrics of the AMVF, the following score function
has been constructed:

$\textrm{Score/Objective
Function}=(\textrm{Eff\\_total}+2\textrm{Eff\\_clean})-(\textrm{Eff\\_split}+\textrm{Eff\\_merge}+\textrm{Eff\\_fake}+\textrm{Resolution})$

###  2.5 Integration of two frameworks

In order to perform the optimization studies, we have integrated the ACTS
tracking framework with our parameter tuning framework as represented in
Figure 2. The tuning framework provides input parameter values to the ACTS
tracking framework. Next, the tracking algorithm runs with these values and
returns the tracking performance to the tuning framework. The tuning framework
computes a score based on the tracking performance and runs a parameter
estimation method to evaluate new configuration of input parameters and pass
it again to the ACTS framework. This continues for a number of iterations
until optimization algorithm achieves satisfactory score and corresponding
parameter configuration with high performance.

Figure 2: A representation of integration of ACTS tracking software framework
with parameter optimization framework.

##  3 Results

Performance metrics for the CKF and the AMVF have been evaluated corresponding
to different optimized parameter configurations and compared with the default
parameter configuration.

###  3.1 CKF Optimization Results

The track reconstruction efficiency and duplicate rate as a function of track
$P_{T}$ and $\eta$ for different parameter configurations is presented in
Figure 3. The black solid circles correspond to the default configuration
while the colored markers represent optimized configurations. A clear
improvement in performance is observed with optimized parameters over all the
$P_{T}$ and $\eta$ range, specially in high $\eta$ range. The table 1 shows
the improvements in terms of the average values of the metrics for different
configurations. The fake rate is negligible for both the default and optimized
parameters.

Figure 3: Comparison of CKF performance between different optimized configurations and default configuration for tracking efficiency vs. $\eta$ (upper left), tracking efficiency vs $P_{T}$ (upper right), duplicate rate vs. $\eta$ (lower left) and duplicate rate vs. $P_{T}$ (lower right). Black solid circles represent performance corresponding to default configuration while colored markers represent optimized configurations.  | Default | EA | Optuna | LIPO | Skopt | Orion
---|---|---|---|---|---|---
Score | 69.06 | 78.88 | 75.38 | 77.28 | 77.05 | 77.79
Efficiency | 93.6% | 96.5% | 96.7% | 96.8% | 96.5% | 96.3%
Duplicate Rate | 72.6% | 56.8% | 59.8% | 60.3% | 56.6% | 58.7%
Fake Rate | 5.56E-03% | 6.2E-03% | 5.2E-03% | 6.8E-03% | 5.7E-03% | 8.8E-03%
time/event (sec) | 50.2 | 31.1 | 46.8 | 37.2 | 40.2 | 33.9
Table 1: CKF performance metrics for default and optimized configurations.

The parameter values obtained using different optimization algorithms are
found to differ from each other as can be seen in Figure 4. This variation
clearly states that there are multiple local minima in the parameter space of
CKF performance.

Figure 4: Parameter configurations obtained for CKF from different
optimization algorithms.

###  3.2 AMVF Optimization Results

The number of reconstructed vertices tagged as total, clean, fake, split and
merged corresponding to optimized parameters are compared with the default
configuration in Figure 5. With the optimized parameters there are more clean
and fewer fake vertices compared to the default. However, there is a trade-off
between the total and split vertices. The position resolution in $x$, $y$ and
$z$ are compared in Figure 6 and shows more centered $z$-resolution for
optimized configurations. The parameter values obtained by different
optimization algorithms are presented in Figure 7, again emphasizing on the
presence of multiple local minima in AMVF space.

Figure 5: Comparison of AMVF performance between different optimized and
default configurations for total (upper left), clean (upper middle), fake
(upper right), split (lower left) and merged (lower right) reconstructed
vertices. Black solid circles represent performance corresponding to the
default configuration while colored ones represent performance corresponding
to the optimized configurations.

Figure 6: Comparison of AMVF resolution in $x$, $y$ and $z$ for different
optimized configurations and default configuration. Black lines represent
resolution corresponding to the default configuration while colored lines
represent the optimized configurations.  Figure 7: Parameter configurations
obtained for the AMVF from the different optimization algorithms.

##  4 CKF optimization using the ATLAS ITk geometry within ACTS

We used the same optimization techniques for the ATLAS ITk geometry [17],
which is shown in Figure 8, within the ACTS framework. We optimized the same
CKF parameters using $t\overline{t}$ data with pile-up = 200 having $P_{T}$
$>$ 1GeV and $|\eta|$ $<$ 4.0. We used the same CKF score function with $k$ =
5.

Figure 8: Diagram showing ATLAS ITk geometry.

The comparison of performance for track reconstruction efficiency and
duplicate rate as a function of $P_{T}$ and $\eta$ is shown in Figure 9. The
black solid circles represent the performance corresponding to the default
unoptimized configuration in ACTS while colored ones represent optimized
performance. The optimization algorithms are able to improve over default
configuration for this geometry, specially in the high-$\eta$ regions.

Figure 9: Comparison of CKF performance for ATLAS ITk geometry between
different optimized configurations and default configuration for tracking
efficiency vs. $\eta$ (upper left), tracking efficiency vs $P_{T}$ (upper
right), duplicate rate vs. $\eta$ (lower left) and duplicate rate vs. $P_{T}$
(lower right). Black solid circles represent performance corresponding to
default configuration while colored marker represent optimized configurations.

##  5 Conclusions

We have presented the implementation and performance of five different
derivative-free optimization algorithms within the context of ACTS software
framework and demonstrated our results using ACTS Generic detector and ATLAS
ITk detector, both integrated within ACTS framework. We have performed our
studies using Combinatorial Kalman Filter and Adaptive Multi Vertex Finder
tracking algorithms. The performance of optimized parameters have been
compared with default configuration currently present within ACTS framework.
Our results show that the different optimization algorithms are able to
automatically find better performing parameter configurations. These studies
demonstrate the potential of such algorithms to be used for the automatic
tuning of tracking algorithms to different detector geometries.

ACKNOWLEDGEMENTS

This work was supported in part by the National Science Foundation under
Cooperative Agreement OAC-1836650.

## References

  * [1] Evolutionary Algorithms, “https://deap.readthedocs.io/en/master/index.html”.
  * [2] Optuna HyperParamter Optimization Framework, “https://optuna.org/”.
  * [3] Cedric Malherbe, Nicolas Vayatis, “Global optimization of Lipschitz functions” [https://arxiv.org/abs/1703.02628].
  * [4] Scikit Optimization, “https://scikit-optimize.github.io/stable/”.
  * [5] Orion HyperParameter Tuning Framework, “https://orion.readthedocs.io/en/stable/index.html”.
  * [6] Peter Chatain, Rocky Garg and Lauren Tompkins, “Evolutionary Algorithms for Tracking Algorithm Parameter Optimization” EPJ Web Conf., 251 (2021) 03071 [https://doi.org/10.1051/epjconf/202125103071].
  * [7] Ai, X., Allaire, C., Calace, N. et al., “A Common Tracking Software Project”, Comput Softw Big Sci 6, 8 (2022), [https://doi.org/10.1007/s41781-021-00078-8].
  * [8] T. Abe, I. Adachi, K. Adamczyk, S. Ahn, H. Aihara, K. Akai, M. Aloi, L. Andricek, K. Aoki, Y. Arai et al., Belle II Technical Design Report (2010) [1011.0352].
  * [9] “Expected tracking and related performance with the updated ATLAS Inner Tracker layout at the High-Luminosity LHC”, ATL-PHYS-PUB-2021-024.
  * [10] Sarah Campbell (for the sPHENIX Collaboration), “sPHENIX: The next generation heavy ion detector at RHIC”, J. Phys.: Conf. Ser. 832 012012 (2017) [https://arxiv.org/abs/1611.03003].
  * [11] T. Akesson, A. Berlin, N. Blinov, O. Colegrove, G. Collura, V. Dutta, B. Echenard, J. Hiltbrand, D.G. Hitlin, J. Incandela et al., Light Dark Matter eXperiment (LDMX) (2018) [1808.05219].
  * [12] Pythia8 Manual, “https://pythia.org/”.
  * [13] M. Kiehn et al, “The TrackML high-energy physics tracking challenge on Kaggle”, EPJ Web of Conferences 214, 06037 (2019), [https://doi.org/10.1051/epjconf/201921406037].
  * [14] P. Billoir, “Progressive track recognition with a Kalman-like fitting procedure”, Comput. Phys. Commun. 57 390-394 (1989) [https://doi.org/10.1016/0010-4655(89)90249-X].
  * [15] G Piacquadio et al, “Primary vertex reconstruction in the ATLAS experiment at LHC” J. Phys.: Conf. Ser. 119 032033 (2008)
  * [16] Pedregosa, F. and Varoquaux, G. et al., “Scikit-learn: Machine Learning in Python”, Journal of Machine Learning Research 12 2825-2830 (2011) [https://scikit-learn.org/].
  * [17] ATLAS Collaboration, “Technical Design Report for the ATLAS Inner Tracker Pixel Detector” CERN-LHCC-2017-021, ATLAS-TDR-030 [https://cds.cern.ch/record/2285585].