††thanks: ††thanks:

# Frequency-robust Mølmer-Sørensen gates via balanced contributions of
multiple motional modes

Brandon P. Ruzic<EMAIL_ADDRESS>Sandia National Laboratories, Albuquerque,
New Mexico 87185, USA  Matthew N. H. Chow<EMAIL_ADDRESS>Sandia National
Laboratories, Albuquerque, New Mexico 87185, USA  Department of Physics and
Astronomy, University of New Mexico, Albuquerque, New Mexico 87106, USA
Center for Quantum Information and Control, University of New Mexico,
Albuquerque, New Mexico 87131, USA  Ashlyn D. Burch  Daniel Lobser  Melissa C.
Revelle  Joshua M. Wilson  Christopher G. Yale  Susan M. Clark  Sandia
National Laboratories, Albuquerque, New Mexico 87185, USA

(August 27, 2024)

###### Abstract

In this work, we design and implement frequency-robust Mølmer-Sørensen gates
on a linear chain of trapped ions, using Gaussian amplitude modulation and a
constant laser frequency. We select this frequency to balance the entanglement
accumulation of all motional modes during the gate to produce a strong
robustness to frequency error, even for long ion chains. We demonstrate this
technique on a three-ion chain, achieving < 1% reduction from peak fidelity
over a 20 kHz range of frequency offset, and we analyze the performance of
this gate design through numerical simulations on chains of two to 33 ions.

††preprint: AIP/123-QED

##  I Introduction

Linear chains of trapped ions are one of the leading platforms for quantum
computation in the near term. The application of Mølmer-Sørensen (MS) gates
[1] on these systems has achieved some of the highest two-qubit entanglement
fidelities to date, reaching above 99.9% while targeting the axial motional
modes of a two-ion chain [2, 3]. To implement powerful quantum algorithms,
like digital quantum simulation [4] and quantum error correction [5, 6, 7],
one must extend these high-fidelity gates to systems of many physical qubits
by, for example, increasing the length of the chain and individually
addressing each ion [8, 9]. In this approach, the MS gates provide all-to-all
connectivity between ion pairs, but the gate fidelity can suffer due to the
residual spin-motion entanglement after the gate in the increased number of
spectator motional modes [1].

There have been many successful demonstrations of high-fidelity MS gates by
modulating the amplitude [8, 9, 10, 11, 12, 13, 14], frequency [15, 16],
amplitude and frequency [17, 18], or phase [19, 20, 21] of the laser beams.
These approaches have achieved 97% to 99.5% fidelity when targeting the radial
modes of a two-ion chain, for which the tighter confinement than in the axial
direction allows better cooling, less heating, and faster gates. The
modulation techniques improve gate performance by eliminating the residual
spin-motion entanglement for ideal experimental conditions and by adding
robustness to this quantity in the presence of motional frequency error. For
example, simulations of frequency-modulated gates maintain a 99% fidelity with
a motional frequency error of $\pm$1.5 kHz for a two-ion chain [15], and
optimizing over a distribution of gate parameters improves this level of
robustness to at least $\pm 5$ kHz [22].

Nevertheless, motional frequency error remains an important error source in MS
gates and their applications. Modulated MS gates attempt to minimize the
sensitivity of the residual spin-motion entanglement to frequency error, and
as a result, the amount of spin entanglement accumulated during the gate also
gains robustness to this error. However, significant errors in the amount of
accumulated spin entanglement can remain and create a purely coherent rotation
error in spin space, which is especially damaging to the performance of
quantum algorithms that involve many gates [23]. This sensitivity to rotation
error was recently demonstrated by the repeated application of MS gates with a
frequency offset on two-ion and four-ion chains [16].

For longer chains, the sensitivity to frequency error increases due to the
higher density of motional modes. Further, the majority of frequency-robust
gate designs become more difficult to implement due to more stringent
experimental requirements, including the need to account for all modes by
linearly increasing the number of optimized pulse-shape parameters with the
number of ions [15]. Robust gate designs exist that reduce this requirement by
only targeting closely spaced ions or a reduced set of motional modes [18],
but the experimental requirements to implement these techniques can still grow
with longer chains. Modulated gates on longer ion chains can require larger
laser powers [16] and generally have a higher sensitivity to drift in the
calibrated model parameters (e.g. motional frequencies, ion separation, laser
power, and gate duration) that are used during the optimization of pulse-shape
parameters [22].

In this paper, we develop and implement an MS gate with an analytic pulse
shape that does not require optimizing a large set of pulse-shape parameters
yet is still broadly robust to motional frequency error, even for long ion
chains. We perform amplitude modulation during our gate with a simple,
Gaussian time dependence that strongly suppresses residual displacement errors
in all modes, as long as the detuning from each mode remains sufficiently
large. While many studies have demonstrated error suppression using amplitude
modulation, including modulation that resembles a Gaussian [17, 18, 14], we
also select a specific, constant detuning that balances the amount of
entanglement accumulation during the gate from all motional modes and provides
robustness to this source of coherent gate error. With the ability to adjust
the detuning without significantly impacting displacement errors, we are free
to tune the laser frequency to a point where the derivative in the
entanglement accumulation with respect to frequency goes to zero. This
produces a gate that is first-order insensitive to frequency error, resulting
in regions of broad robustness to this error. Our protocol is simple to
realize experimentally, as we can optimize performance by calibrating only two
pulse-shape parameters: the constant detuning and the peak Rabi rate. As a
result, our gate design has a low classical computational overhead,
facilitating its adoption on other trapped ion quantum processors and making
it suitable for systems suffering from moderate amounts of drift. We
demonstrate the frequency robustness of our gate on a three-ion chain and
analyze this robustness in numerical simulations for chains of up to 33 ions.

This work is done on the Quantum Scientific Computing Open User Testbed,
QSCOUT. We use qubits encoded in the hyperfine clock states of 171Yb+ ions
trapped in a linear chain on a surface trap. Gates are site-selectively driven
with an optical Raman transition. Details of the apparatus are described in
previous work [24].

##  II Gate Design

###  II.1 MS Gate Model

We model the application of an MS gate on two ions that are part of a linear
chain of ions in a surface trap using the Hamiltonian,

$H(t)=-\Omega(t)\sum_{k}S_{y,k}a_{k}e^{i\delta_{k}t}+h.c.,$ (1)

which is in a rotating frame with respect to the atomic and trap degrees of
freedom. The collective spin operator $S_{y,k}$ has the form:
$S_{y,k}=(\eta_{1,k}\sigma_{y,1}+\eta_{2,k}\sigma_{y,2})/2$, where
$\sigma_{y,j}$ is the $y$ Pauli spin operator for the $j$-th ion targeted by
the gate. The Lamb-Dicke parameter $\eta_{j,k}$ can differ for each ion and
each motional mode, and $\Omega(t)$ is the Rabi rate of the carrier transition
for both ions. In this work, $\Omega(t)$ is a time dependent parameter of the
drive field, while $\delta_{k}$ is effectively held constant in time for each
mode. The operators $a_{k}^{\dagger}$ and $a_{k}$ are the raising and lowering
operators, respectively, for a harmonic oscillator that represents the
motional mode of the ion chain with angular frequency $\nu_{k}$. During the
gate, a dual-tone laser illuminates the ions with detunings
$\pm\delta_{k}=\pm(\delta_{\text{c}}-\nu_{k})$ from their blue and red
motional sideband transitions, respectively, where the parameter
$\delta_{\text{c}}$ is the detuning of the blue-detuned laser tone from the
carrier transition. For simplicity, we have made the Lamb-Dicke approximation:
$e^{i\eta(a_{k}+a_{k}^{\dagger})}\approx 1+i\eta(a_{k}+a_{k}^{\dagger})$. We
have also neglected the carrier transition and the far-off-resonant sideband
transitions with detunings larger than $|\delta_{\text{c}}|$.

Since the Hamiltonian $H(t)$ acts on each motional mode independently, we can
write the propagator $U(t)$ as a product over motional modes:

$U(t)=\Pi_{k}U_{k}(t),$ (2)

and the exact analytic solution for $U_{k}(t)$ is [1, 25],

$\displaystyle
U_{k}(t)=e^{-i\mathcal{B}_{k}(t)S_{y,k}^{2}}D(S_{y,k}\alpha_{k}(t)),$
$\displaystyle\mathcal{B}_{k}(t)=\frac{i}{2}\int_{0}^{t}\left(\frac{\mathrm{d}\alpha_{k}(t^{\prime})}{\mathrm{d}t^{\prime}}\alpha_{k}^{*}(t^{\prime})-\alpha(t^{\prime})\frac{\mathrm{d}\alpha_{k}^{*}(t^{\prime})}{\mathrm{d}t^{\prime}}\right)\mathrm{d}t^{\prime}.$
(3)

The displacement operator
$D(S_{y,k}\alpha_{k}(t))=\exp\left[S_{y,k}(\alpha_{k}(t)a_{k}^{\dagger}-\alpha_{k}^{*}(t)a_{k})\right]$
is conditioned on the spin state of the targeted ions, and $\alpha_{k}(t)$
describes the phase-space trajectory of the ion chain. The phase
$\eta_{1,k}\eta_{2,k}\mathcal{B}_{k}(t)$, which governs the amount of spin
entanglement accrued during the gate, is real and positive (negative) for
clockwise (counter-clockwise) trajectories.

To gain an intuitive picture of the gate dynamics, we express the phase-space
trajectory of each motional mode in terms of the parameters of $H(t)$,

$\displaystyle\alpha_{k}(t)=i\int_{0}^{t}\Omega(t^{\prime})e^{-i\delta_{k}t^{\prime}}\mathrm{d}t^{\prime},$
(4)

where $t=\tau$ corresponds to the end of the gate. From this equation, we see
that $\alpha_{k}(\tau)$ is proportional to the Fourier transform of
$\Omega(t)$ evaluated at $\delta_{k}$, assuming that $\Omega(t)$ is zero
before ($t<0$) and after ($t>\tau$) the gate. This is a key insight that will
aid our choice of pulse shape for frequency-robust gates, as discussed in
section II.3.

In this study, we focus on the robustness of MS-gate performance to a
frequency error $\delta\omega$ that is applied to both laser tones and moves
them symmetrically with respect to the carrier transition, resulting in new
carrier detunings:
$\pm\delta_{\text{c}}^{\prime}=\pm(\delta_{\text{c}}+\delta\omega)$ and
sideband detunings: $\pm\delta_{k}^{\prime}=\pm(\delta_{k}+\delta\omega)$, for
the blue-detuned and red-detuned tones, respectively. Equivalently, this
frequency error can be interpreted as a common change in the motional
frequency of each mode: $\nu_{k}^{\prime}=\nu_{k}-\delta\omega$. Although
other error sources can affect gate performance, such as laser power
fluctuations and anomalous heating [26, 27], we choose to focus on frequency
error due to the high sensitivity of gate performance to this error [25],
especially in the context of long ion chains with many closely spaced motional
modes.

###  II.2 Performance Metrics

We use the state fidelity $\mathcal{F}$ as the figure of merit for gate
performance, which can be computed by wavefunction overlap:

$\mathcal{F}=|\braket{\Phi}{\Psi(\tau)}|^{2},$ (5)

where $\ket{\Psi(\tau)}$ is the wavefunction of the ion chain at the end of
the gate and $\ket{\Phi}$ is the target state. We assume that we perfectly
initialize the ions in the ground spin state $\ket{00}$ and laser cool them to
reach the motional ground state $\ket{0}$. The state of the ion chain after
the gate is then
$\ket{\Psi(\tau)}=\Pi_{k}e^{i\mathcal{B}_{k}(\tau)S_{y,k}^{2}}D(S_{y,k}\alpha_{k}(\tau))\ket{00,0}$,
and we choose to target the state
$\ket{\phi}=1/\sqrt{2}(\ket{00}+i\ket{11})\ket{0}$, a maximally entangled spin
state and the motional ground state, as any residual displacement after the
gate leads to spin-motion entanglement.

For ideal gate performance $(\mathcal{F}=1)$, we require the propagator at the
end of the gate ($t=\tau)$ to take the following form:
$U(\tau)=e^{-i\sigma_{y,1}\sigma_{y,2}\theta/2}$, where $\theta$ is the
rotation angle of the gate. For our choice of initial and final states, the
ideal gate is accomplished for,

$\displaystyle\alpha_{k}(\tau)=0\quad\text{for each $k$},$ (6a)
$\displaystyle\theta=\sum_{k}\eta_{1,k}\eta_{2,k}\mathcal{B}_{k}(\tau)=\pi/2.$
(6b)

Although $\mathcal{F}$ is a sufficient metric for gate performance, we find
that it is illustrative in this work to decompose the state fidelity into two
contributing terms: displacement error $\epsilon_{d}$ and rotation-angle error
$\epsilon_{r}$, which arise from inequalities in equations (6a) and (6b),
respectively. In the next two sections, we derive the contribution to the
state infidelity $1-\mathcal{F}$ from each error separately, where the sum of
these errors,

$\epsilon_{s}=\epsilon_{d}+\epsilon_{r},$ (7)

is approximately equal to the state infidelity: $\epsilon_{s}\approx
1-\mathcal{F}$, for small errors ($\epsilon_{d}\ll 1$ and $\epsilon_{r}\ll
1$).

####  II.2.1 Displacement Error

Displacement error occurs when the coherent displacement at the end of the
gate is non-zero, $\alpha_{k}(\tau)\neq\alpha(0)=0$, and leads to residual
spin-motion entanglement. In the phase-space trajectory picture, a gate with
no displacement error for a particular mode will produce a closed curve, and
any residual displacement (or open curve) contributes to the gate error. We
call this contribution the displacement error $\epsilon_{d,k}$ for mode $k$.
As accomplished in several previous works, including [13, 15, 22, 21], we
design gates that are robust to frequency errors for this error mechanism by
shaping the laser pulse.

We derive an expression for the displacement error from each mode by acting
the displacement operator for the gate on the initial state and computing the
wavefunction overlap with the target state:
$\mathcal{O}_{d,k}=\bra{\Phi}D(S_{y,k}\alpha_{k}(\tau))\ket{00,0}$. In the
spin basis
{$\ket{+}=(\ket{0}+\ket{1})/\sqrt{2},\ket{-}=(\ket{1}-\ket{0})/\sqrt{2}$}, the
operator $S_{y,k}$ has four eigenvalues:
$\lambda_{++}=(\eta_{1,k}+\eta_{2,k})/2$,
$\lambda_{+-}=(\eta_{1,k}-\eta_{2,k})/2$,
$\lambda_{-+}=(\eta_{2,k}-\eta_{1,k})/2$, and
$\lambda_{--}=-(\eta_{1,k}+\eta_{2,k})/2$, corresponding to the spin states
$\ket{++}$, $\ket{+-}$, $\ket{-+}$, and $\ket{--}$, respectively. In terms of
these eigenvalues,

$\mathcal{O}_{d,k}=\frac{1}{4}\sum_{\lambda}e^{-|\lambda\alpha_{k}(\tau)|^{2}/2},$
(8)

and the total displacement error is,

$\displaystyle\epsilon_{d}=\sum_{k}\epsilon_{d,k},$ (9a)
$\displaystyle\epsilon_{d,k}=1-\left|\frac{1}{4}\sum_{\lambda}e^{-|\lambda\alpha_{k}(\tau)|^{2}/2}\right|^{2}.$
(9b)

If we assume the Lamb-Dicke parameters are equal for each ion and mode
involved in the gate ($\eta_{1,k}=\eta_{2,k}=\eta$), then the eigenvalues of
$S_{y,k}$ are $\lambda=\\{\eta,0,0,-\eta\\}$, and from the equation above, we
find that $\epsilon_{d,k}\approx\eta^{2}|\alpha_{k}(\tau)|^{2}/2$ for small
errors.

####  II.2.2 Rotation-angle Error

Rotation-angle error $\epsilon_{r}$ is error in the entangling phase
accumulated at the end of the gate. In the phase-space trajectory picture,
this can be visualized as the area enclosed by the curves. In previous works,
most gate designs do not explicitly target solutions that are robust to
rotation-angle error caused by motional frequency drift. In this work, we
utilize the contributions from multiple motional modes in order to derive a
gate in which the entangling phase is independent to first order in symmetric
detuning offset.

We derive an expression for the rotation error by acting the spin-entangling
operator for the gate on the initial state and computing the wavefunction
overlap with the target state:
$\mathcal{O}_{r}=\bra{\Phi}\Pi_{k}e^{i\mathcal{B}_{k}(\tau)S_{y,k}^{2}}\ket{00,0}$.
Since the spin-entangling operator is independent of the motional state, we
only require that the phase accrued during the gate
$\theta=\sum_{k}\eta_{1,k}\eta_{2,k}\mathcal{B}_{k}(\tau)$ is equal to $\pi/2$
to achieve an MS gate, as shown in equation (6b). Therefore, any additional
phase $\delta\theta=\theta-\pi/2$ results in the rotation $\delta
U=e^{i\delta\theta\sigma_{y,1}\sigma_{y,2}/2}$ of the target state, up to an
arbitrary global phase, and the rotation error is,

$\epsilon_{r}=1-\left|\bra{\Phi}\delta
U\ket{\Phi}\right|^{2}=\frac{1}{4}\left|\delta\theta\right|^{2}.$ (10)

###  II.3 Gaussian Pulse Shape

The robustness of our gate design to motional frequency error relies on the
specific shape of the pulse amplitude during the gate, which we describe by a
time-dependent Rabi rate $\Omega(t)$. From section II.2.1, we see that the
state infidelity grows linearly with the residual displacement of each mode
$|\alpha_{k}(\tau)|^{2}$. To reduce these contributions, we note that
$\alpha_{k}(\tau)$ is proportional to the Fourier transform of $\Omega(t)$
evaluated at $\delta_{k}$, as seen in section II.1, and we choose to employ a
truncated-Gaussian pulse shape,

$\Omega(t)=\begin{cases}\Omega_{0}e^{-(t-\tau/2)^{2}/2z^{2}},&\text{if }0\leq
t\leq\tau,\\\ 0,&\text{otherwise},\end{cases}$ (11)

where $\Omega_{0}$ is the peak Rabi rate, and $z$ is the Gaussian width. To
understand the effect of this pulse shape, we assume that $z$ is sufficiently
small ($z^{2}\ll\tau^{2}/8$) such that we can ignore the truncation of
$\Omega(t)$ and extend the limits of integration in equation (4) for
$\alpha(\tau)$ to all times,

$\displaystyle\alpha_{k}(\tau)\approx
i\int_{-\infty}^{\infty}\Omega^{\prime}(t^{\prime})e^{-i\delta_{k}t^{\prime}}\mathrm{d}t^{\prime},$
(12a) $\displaystyle|\alpha_{k}(\tau)|^{2}\approx
2\pi\Omega_{0}^{2}z^{2}e^{-\delta_{k}^{2}z^{2}},$ (12b)

where $\Omega^{\prime}(t)=\Omega_{0}e^{-(t-\tau/2)^{2}/2z^{2}}$ for all $t$.
Here, we see that the Gaussian-like pulse shape guarantees that the
displacement error from each mode will exponentially decay with
$\delta_{k}^{2}$. For the appropriate choice of detuning, this can strongly
suppress contributions to $\epsilon_{d}$ from the modes primarily driving the
gate and almost entirely eliminate contributions to $\epsilon_{d}$ from far-
detuned "spectator" modes.

###  II.4 Mode Balancing

Many MS-gate implementations, including the standard gate with a constant
detuning and Rabi rate [1], suffer from the displacement error $\epsilon_{d}$
caused by spectator modes, and the gate performance degrades when the detuning
of the mode primarily targeted by the gate has a similar magnitude as the
detuning from the other modes. However, since the Gaussian pulse shape
strongly suppresses $\epsilon_{d}$, as discussed in section II.3, we have the
freedom to detune close to multiple modes without suffering from large
displacement errors. In this section, we take advantage of this freedom and
determine a choice of detuning $\delta_{c}$ that not only maintains the
suppression of $\epsilon_{d}$ but also reduces the rotation-angle error
$\epsilon_{r}$.

While staying sufficiently far away from all modes ($\delta_{k}\gg z^{-1}$) to
suppress $\epsilon_{d}$, we aim to balance the contributions to $\theta$ from
all modes, such that frequency errors that increase the contribution to
$\theta$ from some modes are cancelled to first-order by the decrease in the
contribution to $\theta$ from the other modes. In general, the detunings that
balance the contributions to $\theta$ from all modes can be found by solving
the following equation for $\delta_{\text{c}}$,

$\frac{\mathrm{d}\theta}{\mathrm{d}\delta_{\text{c}}}=\sum_{k}\eta_{1,k}\eta_{2,k}\frac{\mathrm{d}\mathcal{B}_{k}(\tau)}{\mathrm{d}\delta_{\text{c}}}=0.$
(13)

This equation can be solved numerically to implement gates for which
$\epsilon_{r}$ is first-order robust to frequency error. At least for mono-
species ion chains, solutions to equation (13) exist for all pairs of ions.

Although solving equation (13) is straightforward and numerically efficient,
we can search for approximate solutions to further simplify the design and
implementation of this gate. In the region of $\delta_{c}$ between two
neighboring modes $k_{1}$ and $k_{2}$, we can find an approximate solution by
neglecting the contribution to $\mathrm{d}\theta/\mathrm{d}\delta_{\text{c}}$
from all other modes because the $k$-th term in the sum becomes large near
$\delta_{k}=0$ and tends to dominate all other terms. We also find that
$\mathrm{d}\mathcal{B}_{k_{1}}/\mathrm{d}\delta_{\text{c}}$ and
$\mathrm{d}\mathcal{B}_{k_{2}}/\mathrm{d}\delta_{\text{c}}$ have the same sign
in this region, such that solutions to equation (13) are likely to be found
when the products $\eta_{1,k_{1}}\eta_{2,k_{1}}$ and
$\eta_{1,k_{2}}\eta_{2,k_{2}}$ have opposite signs. To the extent that these
products have similar magnitudes, the two modes make similar contributions to
the gate when $|\delta_{k_{1}}|\approx|\delta_{k_{2}}|$, and we can expect to
find a solution to equation (13) that is close to the midpoint between modes
$k_{1}$ and $k_{2}$: $\delta_{\text{c}}=(\nu_{1}+\nu_{2})/2$. In practice, one
can use this value of $\delta_{c}$ as a starting point and experimentally
calibrate $\delta_{c}$ as we describe in section III.3.

###  II.5 Gate Comparison

(a)

(b)

(c)

Figure 1: (a) Displacement error $\epsilon_{d}$, (b) rotation-angle error
$\epsilon_{r}$, and (c) state infidelity
$\epsilon_{s}=\epsilon_{d}+\epsilon_{r}$ for a three-ion chain with three
different pulse shapes. From bottom to top at most detunings, the curves
correspond to (solid blue) a Gaussian pulse with a symmetric detuning of
$\delta_{0}/2\pi=37.2$ kHz, (dashed orange) a Gaussian pulse with
$\delta_{0}/2\pi=-40$ kHz, and (solid green) a square pulse with
$\delta_{0}/2\pi=-40$ kHz. The detuning $\delta_{0}=0$ is resonant with the
lowest frequency motional mode, and the next two higher motional modes are at
$\delta_{0}/2\pi=94.7$ kHz and $\delta_{0}/2\pi=160$ kHz, corresponding to the
peaks in both error metrics. Vertical dash-dotted lines are drawn at (orange)
$\delta_{0}/2\pi=-40$ kHz and (blue) $\delta_{0}/2\pi=37.2$ kHz. The frequency
error $\delta\omega/2\pi$ is defined relative to the optimal value of
$\delta_{0}/2\pi$ for each pulse.

To demonstrate the robustness to frequency error of the balanced Gaussian
gate, we simulate MS gates with three different pulse shapes subject to the
symmetric frequency error $\delta\omega$ and compute the contributions to the
state infidelity from displacement error $\epsilon_{d}$ and rotation-angle
error $\epsilon_{r}$ for each pulse. We compare the balanced Gaussian pulse,
which has a detuning that solves equation (13), with a standard Gaussian pulse
and a square pulse (constant Rabi rate during the gate) that are both detuned
below the lowest motional mode. The duration of each gate is 200 $\mu$s and
the width of each Gaussian pulse is $z=25$ $\mu$s. Being a small fraction of
the pulse length, this choice of $z$ creates Gaussian-like pulses with small
truncation effects.

We provide a concrete example for this comparison by simulating each gate on
the outer ions of a three-ion chain with an ion separation of $4.5$ $\mu$m and
radial trapping frequencies of $\omega_{\text{trap, ra}}/2\pi=2.52$ MHz and
$\omega_{\text{trap,rb}}/2\pi=2.19$ MHz in the orthogonal radial-a and
radial-b directions, respectively. These trapping frequencies correspond to
the motional frequency of the center-of-mass (highest) mode in each direction
and have values that are typical for the surface trap of the QSCOUT testbed.
We also assume the laser $k$-vector is aligned at a 45$\degree$ angle between
the two radial directions, and we neglect any excitation of the axial motional
modes, which will have a large detuning compared to the radial modes during
the gate and are orthogonal to the $k$-vector of the laser.

We design the balanced Gaussian gate by targeting the lowest two radial modes:
the zig-zag ($k=0$) and tilt ($k=1$) radial-b modes, which have a splitting of
$\Delta\nu_{10}/2\pi=(\nu_{1}-\nu_{0})/2\pi=94.7$ kHz. The targeted ions have
Lamb-Dicke parameters that are equal in the zig-zag mode, while being equal
and opposite in the tilt mode, such that the products $\eta_{1,0}\eta_{2,0}$
and $\eta_{1,1}\eta_{2,1}$ have opposite signs and are approximately equal in
magnitude. We then solve equation (13) in the region between these modes and
obtain a detuning of $\delta_{0}/2\pi=37.2$ kHz above the lowest frequency
mode (zig-zag). For the standard Gaussian and square gates, we choose a
detuning of $\delta_{0}/2\pi=-40$ kHz. As a final step, we select the peak
Rabi rate $\Omega_{0}$ for each gate separately to guarantee $\theta=\pi/2$ at
the chosen detuning. The frequency error $\delta\omega$ is defined relative to
these nominal detunings for each gate.

Figure 1 shows the simulated values of $\epsilon_{d}$ and $\epsilon_{r}$ for
the three different pulse shapes as a function of the detuning
$\delta_{0}/2\pi$ from the lowest motional mode. The Gaussian pulse shape
(dashed orange and solid blue) suppresses displacement error as long as the
gate is performed sufficiently far from all motional modes. By contrast, the
square pulse gate (solid green) has narrow minima, and relatively high
displacement error persists at all detunings in the range displayed. A broad
dip in the rotation-angle error is apparent for the balanced Gaussian gate
(solid blue) when the derivative with respect to detuning of the sum of
contributions of multiple modes goes to zero. Intuitively, this can be
described by the fact that as the detuning moves in one direction away from
the optimal point the contribution from one mode becomes smaller but the
contribution from the other mode becomes larger. By contrast, the standard
Gaussian gate (dashed orange) and the square gate (solid green) only have
narrow dips in rotation-angle error and thus have a small detuning range over
which the target phase is accumulated. For small frequency errors, the state
infidelity $\epsilon_{s}=\epsilon_{d}+\epsilon_{r}$ for each Gaussian pulse is
strongly dominated by the contribution from the rotation-error $\epsilon_{r}$.
Figure 1c shows $\epsilon_{s}$ for each pulse shape over a range of
experimentally relevant frequency errors, demonstrating the improved
robustness to $\delta\omega$ when using the balanced Gaussian pulse shape.

##  III Gate Implementation

###  III.1 Gate Parameter Selection

In addition to choosing the correct detuning to balance the contributions of
multiple motional modes, the time-domain standard deviation, $z$, of the
truncated Gaussian pulse shape is a free parameter that may be tuned to
optimize the gate performance. This may be done numerically or empirically
given sufficient intuition about the contributions to gate error. In
particular, we find that optimal value of $z$ for robustness to frequency
error depends on the detuning from the closest two motional modes. On the one
hand, the truncation of the Gaussian lying outside of the pulse duration
$\tau$ leads to some amount of square pulse character with an abrupt pulse
turn on. This truncation effect and square pulse character lead to some
residual displacement error, especially when the Gaussian width becomes
comparable to the gate time ($z/\tau\approx\mathcal{O}(1)$). This error can be
seen as the oscillating floor in the plot of $\epsilon_{d}$ in Fig. 1a. On the
other hand, displacement error can become large when $z$ is small enough that
$e^{-\delta^{2}_{k}z^{2}}$ becomes significant. This can be understood as
Fourier broadening of the pulse as it becomes narrow in time. To optimize
robustness, a value of $z$ can be chosen that balances the infidelity
contribution from cutoff effects and Fourier broadening over a given range of
frequency error. One must also take into consideration that as $z$ is reduced,
the peak intensity of the pulse can become an experimental challenge because
the Rabi rate must be scaled up (by increasing $\Omega_{0}$) to achieve the
proper rotation angle of $|\theta|=\pi/2$.

Additionally, the choice of modes targeted by the gate plays an important role
in its robustness to frequency error. In general, motional modes with a larger
frequency splitting $\Delta\nu/2\pi$ provide more robustness by allowing
larger detunings ($\delta_{k_{1}}$ and $\delta_{k_{2}}$) through which the
displacement error $\epsilon_{d}$ becomes more strongly suppressed. Also,
larger detunings reduce the magnitude of
$\mathrm{d}\theta/\mathrm{d}\delta_{c}$ away from its zero crossing, providing
a reduction in the rotation-angle error $\epsilon_{r}$ for finite frequency
errors. For typical ion chains in a harmonic well, the lowest two radial modes
will have the largest splitting and therefore the maximum robustness to
frequency error. This situation is reversed in the case of axial modes, for
which the highest two modes typically have the largest splitting.

Because the gate parameters $\delta_{c}$, $z$, and $\Omega_{0}$ depend on the
specific distribution of motional modes, we provide a concrete example of
optimal parameter selection for the model of a three-ion chain described in
section II.5, which targets the outer ions and the lowest two radial modes. To
determine the optimal gate parameters and assess their robustness for this
model, we numerically vary the Gaussian width $z$ for a fixed pulse duration
of $\tau=200\,\mu$s and show the performance of the gate over a range of
$\delta\omega=\pm 10\,$kHz in Fig. 2. At large values of $z$ in this plot, the
pulse resembles a square pulse with an approximately constant magnitude during
the gate and a strong truncation effect. In this regime, we recover an
infidelity proportional to ${\rm sinc}(\delta\omega\tau)$, as we expect for
the Fourier transform of a square pulse. Toward the lower end of $z$ on the
plot, the pulse resembles a Gaussian with a small truncation effect, and we
find that there is a broad region of $z$ and $\delta\omega$ where the gate
performs well. For example, the state infidelity $\epsilon_{s}$ remains below
$10^{-3}$ over the range: $-7.8$ kHz $\leq\delta_{\omega}/2\pi\leq 8.5$ kHz
for $z=25$ $\mu$s and over the range: $13$ $\mu$s $\leq z\leq 44$ $\mu$s for
$\delta\omega=0$.

Figure 2: State infidelity $\epsilon_{s}$ as a function of the Gaussian width
$z$ and symmetric detuning offset $\delta\omega$ for a truncated, balanced
Gaussian MS gate on a three-ion chain. As the gate duration is $\tau=200$
$\mu$s, the gate varies from delta-function to square-pulse character from
left to right.

(a)

(b)

(c)

Figure 3: (a) A symmetric detuning scan between two motional modes (marked
with dashed lines at 0 kHz and $\approx$110 kHz) shows a broad region around
$\delta_{0}/2\pi=46$ kHz where the balanced Gaussian gate performs well
(shaded region). Lines top to bottom at $\delta_{0}/2\pi<-25\,$kHz are the
detection probabilities $P(00)$, $P(11)$, and $P(01)+P(10)$ after preparing in
$\ket{00}$ and applying the gate. (b) Parity oscillations are measured
(points) by applying a single-qubit rotation with variable phase ($\phi$) and
fit (lines) to a sine. From light to dark, the symmetric detuning offset
increases over the range $\delta_{\omega}/2\pi=\pm$10 kHz for the balanced
Gaussian gate. (c) Measured entangling gate fidelity ($\mathcal{F}$) as a
function of symmetric detuning offset ($\delta_{\omega}$) for the balanced
Gaussian (blue circles) and unbalanced Gaussian (orange squares). Fidelity is
calculated according to Eq. 14 using the parity scan results and population
measurements. Uncertainty is derived from a Wilson score interval on the
population measurement and fitting uncertainty on the parity scan.

###  III.2 Experimental Implementation

We implement the derived gate on a chain of ${N=3}$ 171Yb+ ions to measure
robustness to symmetric detuning offsets. We use the hyperfine ground, “clock”
states as the qubit levels: $\ket{F=0,m_{F}=0}\equiv\ket{0}$ and
$\ket{F=1,m_{F}=0}\equiv\ket{1}$. In all experiments, all ions are initialized
in $\ket{0}$ and the gate under study is applied to two target ions. A global,
single-qubit $\pi/2$ rotation is then applied in parity scan measurements
(described in III.3). Finally, in all measurements, the population of each
qubit state is determined by fluorescence detection. The relevant two-ion
state is labelled as ${\rm
Tr}_{c}(\ket{a}\otimes\ket{b}\otimes\ket{c})\equiv\ket{ab}$ where $a$ and $b$
are the states of the two target ions, and $c$ is the state of the third
“spectator” ion in the chain, which is ignored.

We set the principal axes of the trap to be at a 45$\degree$ angle from the
effective Raman $k$-vector to allow for Raman sideband cooling on all radial
modes, and thus observe a total of $2N$ radial motional frequencies. For the
fidelity measurements presented here, we have radial motional frequencies of
$\nu_{\text{ra}}/2\pi=\\{2.134,2.229,2.296\\}$ MHz and
$\nu_{\text{rb}}/2\pi=\\{1.832,1.941,2.017\\}$ MHz, and an axial center-of-
mass frequency of $\nu_{\text{axial}}/2\pi=$ 0.52 MHz. As modelled previously,
we operate between the lowest two radial-b modes.

The Gaussian pulse shape is approximated by a natural cubic spline with 13
amplitude knots that are passed to our custom Radio Frequency System-on-Chip
(RFSoC) hardware, “Octet” [24]. Octet generates the RF waveform that drives
the acousto-optic modulators (AOMs) which perform the required RF to optical
transduction. The spline knots are equally spaced along the square root of a
Gaussian pulse shape and applied to both the individual addressing beams and
counterpropagating global beam AOMs, thus producing a Gaussian temporal
profile in the two-photon Raman Rabi rate of the pulse. The first and last
knots of the spline are non-zero, and truncate parts of the infinite Gaussian
lying outside of the pulse duration.

We choose to implement a gate with
$z/\tau=26.5\,\mu\text{s}/200\,\mu\text{s}\approx 0.13$ on the experiment to
maintain reasonably low truncation error while still avoiding demanding Rabi
rates. A gate time of 200 $\mu$s is chosen empirically to minimize the effects
of heating, truncation, and Fourier broadening. The “balance point” detuning
($\delta_{0}$) is found by scanning the symmetric detuning and finding the
point of zero slope in the $\ket{11}$ population, and the peak Rabi rate
($\Omega_{0}$) is set by varying the scaling of the Gaussian pulse shape and
finding the point of equal $\ket{00}$ and $\ket{11}$ populations. For the
balanced Gaussian gate demonstrated here, $\delta_{0}/2\pi=46\,$kHz and
$\Omega_{0}\approx 2\pi\times 190\,$kHz.

Although the wavelength (355 nm) of the Raman laser is chosen to approximately
balance the AC Stark Shift on the qubit transition, there remains a residual
differential AC Stark Shift on the order of 1 kHz. To compensate for this
shift, we dynamically apply a virtual frame rotation at a rate proportional to
the intensity of the laser applied to each ion. The magnitude of the frame
rotation is found by preparing in the $\ket{00}$ state, applying two MS gates
back to back, and maximizing the resultant population in the $\ket{11}$ state.

###  III.3 Experimental Results

We measure the performance of the balanced Gaussian MS gate and compare it to
the unbalanced Gaussian MS gates. While our heating rate is high and ultimate
fidelity of both gates appears to be dominated by incoherent heating errors,
there is clear difference in response to intentionally applied symmetric
detuning offsets. At several points in the measured range of detuning offsets,
we measure both the even parity population after an MS gate and the amplitude
of a parity oscillation acquired by applying a single-qubit $\pi/2$ pulse with
variable phase after the MS gate. We then estimate the fidelity
($\mathcal{F}$) according to [28, 29, 30, 31],

$\mathcal{F}=\frac{1}{2}(\rho_{00}+\rho_{11})+\frac{1}{2}A_{\pi},$ (14)

where $\rho_{ab}$ is the population of the $\ket{ab}$ state after the MS gate
and $A_{\pi}$ is the amplitude of the parity oscillation.

We apply the gate to the outer two ions on a three-ion chain, operating
between the zig-zag and tilt modes and characterize its performance, as
presented in Fig. 3. In a scan of the symmetric detuning offset, we prepare
the ions in $\ket{00}$, the gate is applied, and the populations are read out.
The even parity populations serve as an indicator for gate angle, and show the
correct gate angle is achieved in a range of $\approx$ 20 kHz centered around
$\delta_{0}/2\pi=46$ kHz relative to the zig-zag mode, when ${\rho_{00}\approx
0.5\approx\rho_{11}}$. The odd parity population ($\rho_{01}+\rho_{10}$) is an
indicator of displacement error, and shows good performance (population is
near zero) as long as the detuning is sufficiently far from the motional
modes. We then estimate the fidelity of the gate at various symmetric detuning
offsets by taking a parity curve and a population measurement with 4000 shots
at each detuning. We find that the fidelity of the balanced Gaussian MS gate
drops by < 1% over $\delta_{\omega}/2\pi\leq\pm 10\,$kHz, indicating broad
robustness to symmetric detuning offset. By contrast, the fidelity of an
unbalanced Gaussian gate (operated below the lowest mode,
$\delta_{0}/2\pi=-25.3\,$kHz 111The laser power required to perform the
unbalanced Gaussian MS gate is significantly higher since there is only strong
contribution from one motional mode instead of two. As such, we are limited to
a detuning of -25.3 kHz to stay within the bounds of our laser system. We also
use the center ion and an edge ion in the unbalanced case to take advantage of
the stronger zig-zag coupling from the center ion.) drops by > 3% even from an
offset of $\pm 1$ kHz. While the peak fidelity of the balanced Gaussian is
only 97.7${}^{+0.3}_{-0.4}$%, we currently have relatively high heating rates
(1000s of quanta/s) compared to contemporary systems [33, 34], and we believe
that our peak fidelity is limited mostly by incoherent heating errors.

###  III.4 Extension to Larger Numbers of Ions

In addition to numerically and experimentally verifying performance of the
balanced Gaussian gate on a chain of three ions, we also numerically explore
the simulated MS gate performance on chains of variable length from $N=2$ up
to $N=33$. As described in section III.1, we target the lowest two motional
modes ($k=0$ and $k=1$) of the $N$-ion chain because we expect to achieve more
robustness to frequency error $\delta\omega$ by targeting the neighboring
modes with the largest frequency splitting. In the region between these modes,
we will find a detuning $\delta_{c}$ that solves equation (13) and balances
the contributions to $\epsilon_{r}$ from all modes, which is possible when the
products $\eta_{1,0}\eta_{2,0}$ and $\eta_{1,1}\eta_{2,1}$ have opposite
signs.

To maintain consistency in the sign and magnitude of the Lamb-Dicke parameters
between chains of different $N$, we target the ions to the immediate left and
right of center, for which we find a solution to equation (13) in each case.
The targeted ions are the two center ions of the chain for even $N$ and the
center ion and one of its neighbors for odd $N$. For this choice of ions,
targeting the lowest two modes leads to a balanced gate with a value of
$\delta_{c}$ that solves equation (13); however, $\theta$ has the opposite
sign for even $N$ than for odd $N$. To maintain the same value of the target
propagator $U(\tau)=e^{-i\sigma_{y,1}\sigma_{y,2}\theta/2}$ (and hence the
same target gate) in our simulations for both even and odd $N$, we impose a
differential laser phase of $\pi$ between the two ions for even $N$. We model
this change in phase by $\sigma_{y,2}\rightarrow-\sigma_{y,2}$, which modifies
the handedness of the phase-space trajectories and has the same effect on
$U(\tau)$ as $\theta\rightarrow-\theta$.

As $N$ grows, we relax the axial frequency of the harmonic potential such that
at all $N$, the separation between the equilibrium positions of the two center
ions (the center ion and its neighbors) for even (odd) $N$ is a fixed value.
We call this quantity the center ion separation parameter, $\Delta x_{0}$. We
make this choice to keep a fixed center ion separation parameter at the
expense of varying axial frequency in order to maintain consistency with an
experimental apparatus designed for a fixed individual addressing beam
separation. This choice comes at the cost of an increasingly dense radial
motional mode sideband spectrum as $N$ grows because the spacing of the radial
modes is set by the ratio of the axial to radial trap frequencies. We also
note that at any particular $N$, the separation of ions will increase from
center to edge of the chain, as we model a simple harmonic potential in the
axial direction. While anharmonic potentials have been considered in other
works to keep spacing constant within a chain [35, 36, 37], we deem the
implementation of a non-harmonic potential to be outside the scope of this
investigation. We also use the same values of $\tau=200$ $\mu$s and $z=25$
$\mu$s for each $N$ and $\Delta x_{0}$.

(a)

(b)

(c)

(d)

Figure 4: (a and c) The rotation-angle error $\epsilon_{r}$ and (b and d) the
displacement error $\epsilon_{d}$ for a center ion separation parameter of
$\Delta x_{0}=3.0$ $\mu$m. The sum of $\epsilon_{r}$ and $\epsilon_{d}$ is the
state infidelity $\epsilon_{s}=\epsilon_{d}+\epsilon_{r}$. From light to dark,
the curves correspond to (a and b) even $N$ from 2 to 32 in steps of 6 and to
(c and d) odd $N$ from 3 to 33 in steps of 6. For most $N$ shown here,
$\epsilon_{r}$ dominates as the main limiting factor for gate robustness to
detuning error. For large $N$, the displacement error $\epsilon_{d}$
eventually lifts off the Gaussian truncation floor as the motional modes bunch
closer together and significantly reduce error suppression from the factor of
$e^{-\delta^{2}_{k}z^{2}}$.

Figure 4 shows the contributions $\epsilon_{r}$ and $\epsilon_{d}$ to the
state infidelity $\epsilon_{s}=\epsilon_{d}+\epsilon_{r}$ over an
experimentally relevant range of frequency errors $\delta\omega/2\pi$ for ion
chains of length $N=2$ to $N=33$ with $\Delta x_{0}=3$ $\mu$m. Here, we have
reduced the value of $\Delta x_{0}$ from the previous sections to achieve
frequency-robust gates for all $N$. We can see that the robustness to
$\delta\omega$ generally decreases with increasing $N$, but the gate error has
a slightly different dependence on $\delta\omega$ for even and odd $N$. This
difference comes from the Lamb-Dicke parameters of the lowest two modes for
each chain, which cause the value of $\delta_{c}$ that solves equation (13) to
lie much closer to the second lowest mode ($k=1$) for even $N$ than for odd
$N$. We can see that while the robustness of $\epsilon_{r}$ is approximately
equal for similar length chains of even and odd $N$, the smaller value of
$\delta_{1}$ for even $N$ causes $\epsilon_{d}$ over this range of
$\delta\omega$ to be dramatically larger than for odd $N$.

Figure 4 also shows that $\epsilon_{r}$ strongly dominates $\epsilon_{d}$ when
$N\lesssim 20$ (for even $N$) and when $N\lesssim 27$ (for odd $N$). For $N$
greater than these values, $\delta_{1}$ becomes small enough that
$\epsilon_{d}$ makes a significant contribution to $\epsilon_{s}$.
Consequently, although $\epsilon_{r}$ is minimized at $\delta\omega=0$, the
minimum of $\epsilon_{s}$ vs. $\delta\omega$ is shifted away from
$\delta\omega=0$ as $N$ is increased. This effect is consistent with our gate
design since we chose to maximize the robustness of $\epsilon_{r}$ with
respect to $\delta\omega$ instead of minimizing $\epsilon_{s}$. Although the
sensitivity to frequency error grows with increasing $N$, the balanced
Gaussian gate achieves a remarkable robustness to frequency error with
$\epsilon_{s}$ remaining below the $10^{-2}$ level for $\delta\omega=\pm 10$
kHz for ion chains of length $N=2$ to 33.

(a)

(b)

Figure 5: The sensitivity to frequency error, as quantified by the maximum
state infidelity $\epsilon_{s,\text{max}}$ over a range of $\pm 3$ kHz from
the minimum of $\epsilon_{s}$ vs. $\delta\omega$ for (a) even $N$ from 2 to 32
and for (b) odd $N$ from 3 to 33. The motional frequency splitting
$\Delta\nu_{10}$ decreases monotonically with increasing $N$. From light to
dark in each plot, the center ion separation parameter is $\Delta x_{0}=3$
$\mu$m to $4.5$ $\mu$m in steps of $0.5$ $\mu$m. The green diamonds indicate
(a) $N=2$ and (b) $N=3$ for each $\Delta x_{0}$. These plots show that
$\Delta\nu_{10}$ is a good predictor of the sensitivity to motional frequency
error for the range of parameters considered in this work.

To quantify the sensitivity to frequency error for each length of chain, we
compute the maximum state infidelity $\epsilon_{s,\text{max}}$ over a $\pm 3$
kHz range from the frequency that minimizes $\epsilon_{s}$. Figure 5 shows
this measure of sensitivity for even and odd $N$, separately, and for $\Delta
x_{0}=3.5\,\mu$m to $4.5\,\mu$m. For both even and odd $N$,
$\epsilon_{s,\text{max}}$ is approximately a function of only the splitting
between the lowest two modes $\Delta\nu_{10}$ for each value of $\Delta
x_{0}$, despite the different motional frequencies, Lamb-Dicke parameters,
number of ions $N$, and pulse shape parameters $\delta_{c}$ and $\Omega_{0}$.
Since $\Delta\nu_{10}$ appears to be a good predictor of balanced Gaussian
gate performance for the chains considered in this work, we note that designs
of new experiments intending to implement these gates should consider the mode
spacing when choosing a value of $\Delta x_{0}$ and $N$. The value of
$\Delta\nu_{10}$ decreases monotonically with larger $N$ (for fixed $\Delta
x_{0}$) and increases for smaller $\Delta x_{0}$ (at fixed $N$).

The approximate correspondence between $\epsilon_{s,\text{max}}$ and
$\Delta\nu_{10}$ arises because $\epsilon_{s}$ is often dominated by the
contribution from $\epsilon_{r}$ for the chains we have considered. When
$\epsilon_{r}$ dominates, the sensitivity depends strongly on the magnitude of
$\mathrm{d}^{2}\theta/\mathrm{d}\delta\omega^{2}$ at $\delta\omega=0$, which
is set by $\Delta\nu_{10}$. This approximate correspondence breaks down when
$\epsilon_{d}$ becomes significant, as $\epsilon_{d}$ is set by the detuning
$\delta_{k}$ instead of $\Delta\nu_{10}$. In such scenarios,
$\epsilon_{s,\text{max}}$ can be substantially different for chains with
different $\Delta x_{0}$ but the same $\Delta\nu_{10}$. For example, for
$\Delta x_{0}=3\,\mu$m, the lightest curves in Fig. 5, we find increased
$\epsilon_{s,\text{max}}$ for long chains (low $\Delta\nu_{10}$) as compared
to chains with $\Delta x_{0}=3.5\,\mu$m to $4.5\,\mu$m for the same values of
$\Delta\nu_{10}$. This effect is much more pronounced for even $N$ than for
odd $N$, consistent with the greater contribution from $\epsilon_{d}$ for even
$N$.

In this extension to larger $N$, we have reduced the ion separation parameter
$\Delta x_{0}$ from our experimental implementation for $N=3$ to increase
$\Delta\nu_{10}$ and improve the performance of MS gates at experimentally
relevant frequency errors for large $N$. Although this improves the robustness
of our simulated gates, a reduction in $\Delta x_{0}$ could have detrimental
effects on quantum circuit performance that are not included in our
simulations by, for example, increasing the level of crosstalk between
neighboring qubits [38]. However, there has been promising research on
applying spin-echo pulses to cancel out the effects of crosstalk when
individual ions cannot be perfectly resolved [39]. Alternatively, one could
decrease the radial trap frequency to increase $\Delta\nu_{10}$, but one needs
to balance this increased robustness to frequency noise with the additional
amount of anomalous heating experienced by lower frequency motional modes.

##  IV Conclusion and Outlook

In summary, we have designed an MS gate that shows broad robustness to
motional frequency errors, addressing a key error mechanism in trapped ion
entangling gates. Further, our gate is performed with a simple analytic pulse
shape, and therefore requires little in the way of computational overhead and
technical complexity. In conjunction with this low cost, the prospects for
scaling to a larger number of ions present a promising outlook for
implementation of this gate on contemporary and next generation trapped ion
systems.

##  V Acknowledgments

This research was supported by the U.S. Department of Energy, Office of
Science, Office of Advanced Scientific Computing Research Quantum Testbed
Program. Sandia National Laboratories is a multimission laboratory managed and
operated by National Technology & Engineering Solutions of Sandia, LLC, a
wholly owned subsidiary of Honeywell International Inc., for the U.S.
Department of Energy’s National Nuclear Security Administration under contract
DE-NA0003525. This paper describes objective technical results and analysis.
Any subjective views or opinions that might be expressed in the paper do not
necessarily represent the views of the U.S. Department of Energy or the United
States Government.

B.P.R. and M.N.H.C. contributed equally to this work.

## References

  * Sørensen and Mølmer [2000] A. Sørensen and K. Mølmer, Phys. Rev. A 62, 022311 (2000).
  * Ballance  _et al._ [2016] C. J. Ballance, T. P. Harty, N. M. Linke, M. A. Sepiol, and D. M. Lucas, Phys. Rev. Lett. 117, 060504 (2016).
  * Gaebler  _et al._ [2016] J. P. Gaebler, T. R. Tan, Y. Lin, Y. Wan, R. Bowler, A. C. Keith, S. Glancy, K. Coakley, E. Knill, D. Leibfried, and D. J. Wineland, Phys. Rev. Lett. 117, 060505 (2016).
  * Lanyon  _et al._ [2011] B. P. Lanyon, C. Hempel, D. Nigg, M. Müller, R. Gerritsma, F. Zähringer, P. Schindler, J. T. Barreiro, M. Rambach, G. Kirchmair, M. Hennrich, P. Zoller, R. Blatt, and C. F. Roos, Science 334, 57 (2011).
  * Shor [1995] P. W. Shor, Phys. Rev. A 52, R2493 (1995).
  * Kitaev [2003] A. Kitaev, Annals of Physics 303, 2 (2003).
  * Aharonov and Ben-Or [2008] D. Aharonov and M. Ben-Or, SIAM Journal on Computing 38, 1207 (2008).
  * Debnath  _et al._ [2016] S. Debnath, N. M. Linke, C. Figgatt, K. A. Landsman, K. Wright, and C. Monroe, Nature 536, 63 (2016).
  * Wright  _et al._ [2019] K. Wright, K. M. Beck, S. Debnath, J. M. Amini, Y. Nam, N. Grzesiak, J.-S. Chen, N. C. Pisenti, M. Chmielewski, C. Collins, K. M. Hudek, J. Mizrahi, J. D. Wong-Campos, S. Allen, J. Apisdorf, P. Solomon, M. Williams, A. M. Ducore, A. Blinov, S. M. Kreikemeier, V. Chaplin, M. Keesan, C. Monroe, and J. Kim, Nature Communications 10, 5464 (2019).
  * Zhu  _et al._ [2006] S.-L. Zhu, C. Monroe, and L.-M. Duan, Europhysics Letters (EPL) 73, 485 (2006).
  * Roos [2008] C. F. Roos, New Journal of Physics 10, 013002 (2008).
  * Benhelm  _et al._ [2008] J. Benhelm, G. Kirchmair, C. F. Roos, and R. Blatt, Nature Physics 4, 463 (2008).
  * Choi  _et al._ [2014] T. Choi, S. Debnath, T. A. Manning, C. Figgatt, Z.-X. Gong, L.-M. Duan, and C. Monroe, Phys. Rev. Lett. 112, 190502 (2014).
  * Tinkey  _et al._ [2022] H. N. Tinkey, C. R. Clark, B. C. Sawyer, and K. R. Brown, Phys. Rev. Lett. 128, 050502 (2022).
  * Leung  _et al._ [2018] P. H. Leung, K. A. Landsman, C. Figgatt, N. M. Linke, C. Monroe, and K. R. Brown, Phys. Rev. Lett. 120, 020501 (2018).
  * Wang  _et al._ [2020] Y. Wang, S. Crain, C. Fang, B. Zhang, S. Huang, Q. Liang, P. H. Leung, K. R. Brown, and J. Kim, Phys. Rev. Lett. 125, 150505 (2020).
  * Leung and Brown [2018] P. H. Leung and K. R. Brown, Phys. Rev. A 98, 032318 (2018).
  * Landsman  _et al._ [2019] K. A. Landsman, Y. Wu, P. H. Leung, D. Zhu, N. M. Linke, K. R. Brown, L. Duan, and C. Monroe, Phys. Rev. A 100, 022332 (2019).
  * Green and Biercuk [2015] T. J. Green and M. J. Biercuk, Phys. Rev. Lett. 114, 120502 (2015).
  * Lu  _et al._ [2019] Y. Lu, S. Zhang, K. Zhang, W. Chen, Y. Shen, J. Zhang, J.-N. Zhang, and K. Kim, Nature 572, 363 (2019).
  * Milne  _et al._ [2020] A. R. Milne, C. L. Edmunds, C. Hempel, F. Roy, S. Mavadia, and M. J. Biercuk, Phys. Rev. Applied 13, 024022 (2020).
  * Kang  _et al._ [2021] M. Kang, Q. Liang, B. Zhang, S. Huang, Y. Wang, C. Fang, J. Kim, and K. R. Brown, Phys. Rev. Applied 16, 024039 (2021).
  * Iverson and Preskill [2020] J. K. Iverson and J. Preskill, New Journal of Physics 22, 073066 (2020).
  * Clark  _et al._ [2021] S. M. Clark, D. Lobser, M. C. Revelle, C. G. Yale, D. Bossert, A. D. Burch, M. N. Chow, C. W. Hogle, M. Ivory, J. Pehr, B. Salzbrenner, D. Stick, W. Sweatt, J. M. Wilson, E. Winrow, and P. Maunz, IEEE Transactions on Quantum Engineering 2, 1 (2021).
  * Ruzic  _et al._ [2022] B. P. Ruzic, T. A. Barrick, J. D. Hunker, R. J. Law, B. K. McFarland, H. J. McGuinness, L. P. Parazzoli, J. D. Sterk, J. W. Van Der Wall, and D. Stick, Phys. Rev. A 105, 052409 (2022).
  * Bruzewicz  _et al._ [2015] C. D. Bruzewicz, J. M. Sage, and J. Chiaverini, Phys. Rev. A 91, 041402 (2015).
  * Boldin  _et al._ [2018] I. A. Boldin, A. Kraft, and C. Wunderlich, Phys. Rev. Lett. 120, 023201 (2018).
  * Sackett  _et al._ [2000] C. A. Sackett, D. Kielpinski, B. E. King, C. Langer, V. Meyer, C. J. Myatt, M. Rowe, Q. A. Turchette, W. M. Itano, D. J. Wineland, and C. Monroe, Nature 404, 256 (2000).
  * Kim  _et al._ [2009] K. Kim, M.-S. Chang, R. Islam, S. Korenblit, L.-M. Duan, and C. Monroe, Phys. Rev. Lett. 103, 120502 (2009).
  * Manning [2014] T. A. Manning,  _Quantum Information Processing with Trapped Ion Chains_ , Ph.D. thesis, University of Maryland, College Park (2014).
  * Figgatt  _et al._ [2019] C. Figgatt, A. Ostrander, N. M. Linke, K. A. Landsman, D. Zhu, D. Maslov, and C. Monroe, Nature 572, 368 (2019).
  * Note [1] The laser power required to perform the unbalanced Gaussian MS gate is significantly higher since there is only strong contribution from one motional mode instead of two. As such, we are limited to a detuning of -25.3 kHz to stay within the bounds of our laser system. We also use the center ion and an edge ion in the unbalanced case to take advantage of the stronger zig-zag coupling from the center ion.
  * Egan [2021] L. N. Egan,  _Scaling Quantum Computers with Long Chains of Trapped Ions_ , Ph.D. thesis, University of Maryland, College Park (2021).
  * Brown  _et al._ [2021] K. R. Brown, J. Chiaverini, J. M. Sage, and H. Häffner, Nature Reviews Materials 6, 892 (2021).
  * Lin  _et al._ [2009] G.-D. Lin, S.-L. Zhu, R. Islam, K. Kim, M.-S. Chang, S. Korenblit, C. Monroe, and L.-M. Duan, EPL (Europhysics Letters) 86, 60004 (2009).
  * Johanning [2016] M. Johanning, Applied Physics B 122, 71 (2016).
  * Xie  _et al._ [2017] Y. Xie, X. Zhang, B. Ou, T. Chen, J. Zhang, C. Wu, W. Wu, and P. Chen, Phys. Rev. A 95, 032341 (2017).
  * Crain  _et al._ [2014] S. Crain, E. Mount, S. Baek, and J. Kim, Applied Physics Letters 105, 181115 (2014).
  * Fang  _et al._ [2022] C. Fang, Y. Wang, S. Huang, K. R. Brown, and J. Kim, arXiv 10.48550/ARXIV.2206.02703 (2022).


# Local clustering coefficient based on three-way partial correlations in
climate networks as a new marker of tropical cyclone

Mikhail Krivonosov  Lobachevsky State University of Nizhny Novgorod, Nizhny
Novgorod, Russia  Olga Vershinina  Lobachevsky State University of Nizhny
Novgorod, Nizhny Novgorod, Russia  Anna Pirova  Lobachevsky State University
of Nizhny Novgorod, Nizhny Novgorod, Russia  Shraddha Gupta  Potsdam Institute
for Climate Impact Research, PO Box 601203, 14412 Potsdam, Germany  Department
of Physics, Humboldt University, Berlin, Germany  Oleg Kanakov  Lobachevsky
State University of Nizhny Novgorod, Nizhny Novgorod, Russia  Jürgen Kurths
Lobachevsky State University of Nizhny Novgorod, Nizhny Novgorod, Russia
Potsdam Institute for Climate Impact Research, PO Box 601203, 14412 Potsdam,
Germany  Department of Physics, Humboldt University, Berlin, Germany

###### Abstract

We introduce a new network marker for climate network analysis. It is based
upon an available special definition of local clustering coefficient for
weighted correlation networks, which was previously introduced in the
neuroscience context and aimed at compensating for uninformative correlations
caused by indirect interactions. We modify this definition further by
replacing Pearson’s pairwise correlation coefficients and Pearson’s three-way
partial correlation coefficients by the respective Kendall’s rank
correlations. This reduces statistical sample size requirements to compute the
correlations, which translates into the possibility of using shorter time
windows and hence into shorter response time of the real-time climate network
analysis. We compare this proposed network marker to the conventional local
clustering coefficient based on unweighted networks obtained by thresholding
the correlation matrix. We show several examples where the new marker is found
to be better associated to tropical cyclones than the unweighted local
clustering coefficient.

##  I Introduction

Network analysis is recognized as a powerful tool in climate science nowadays,
in essence being a specialization of the general correlation networks approach
Langfelder and Horvath (2008), which is also well established in systems
biology Bruggeman and Westerhoff (2007); Van der Greef  _et al._ (2007);
Friedman and Alm (2012); Rosato  _et al._ (2018) and neuroscience Fransson and
Marrelec (2008); Bassett and Sporns (2017). Correlation networks in climate
science are constructed based upon cross-correlations in multivariate time
series of chosen climate variables between nodes of a spatial grid (e.g. a
geographic coordinate grid on Earth surface). Such networks are then used to
compute various graph metrics, which in turn become input features for
subsequent analysis to identify patterns of climate dynamics Tsonis  _et al._
(2006); Donges  _et al._ (2009); Gozolchiani  _et al._ (2011); Dijkstra  _et
al._ (2019a).

It was suggested in Feng  _et al._ (2016) and in follow-up publications (see
Dijkstra  _et al._ (2019b) for review) to use climate network metrics as input
features to machine learning algorithms, and this approach was successfully
applied to predicting El Niño phases. From the machine learning perspective,
network analysis can be viewed as a preprocessing stage to reduce the
dimensionality of input data. Without such preprocessing, the input data at
any specified moment of observation would include all data in a time section
of specific duration (sliding time window preceding the current moment) of the
multivariate time series. The dimensionality of such data is the product of
the number of spatial grid points and the number of time samples within the
window. Network analysis reduces this dataset to a few graph metrics. In case
of local metrics, which are computed for a specific geographic location
associated to a network node (e.g. node degree or local clustering
coefficient), machine learning methods can be applied also locally, with the
set of local network metrics at a specified geographic point taken as input
features, while the metrics themselves are computed, generally speaking, based
on full multivariate input data. Thus, network analysis in the context of
machine learning allows to cut down the dimensionality of the learning problem
by orders of magnitude, while still retaining the dependence of the outcome
upon all available data.

It is essential for the efficiency of this approach that as much as possible
of information contained in the initial data be retained in the feature set
produced by network analysis. A common practice in the literature on climate
networks is to construct an unweighted graph based upon the correlation matrix
of the input data by simple thresholding: network nodes are associated to the
spatial grid nodes of the input data, and a link between two specific nodes is
assumed to exist if the corresponding component of the correlation matrix
(i.e. the correlation coefficient of the corresponding variables) exceeds a
chosen threshold. The threshold value is a free parameter of the method, and
generally allows tuning to maximize performance (in whatever quantitative
sense), but it is also common to fix the threshold value by specifying the
edge density (defined as the fraction which linked pairs of nodes constitute
among the total number of node pairs), e.g. at 5%, without any further
optimization. Regardless of whether the threshold value is optimized or not,
thresholding inevitably leads to information loss.

A known way to mitigate the information loss introduced by thresholding is to
construct and analyze an ensemble of several unweighted graphs at once, by
using a series of threshold values, which has been successfully applied to
genetic network analysis Bockmayr  _et al._ (2013). It is also possible to
eliminate the thresholding operation at all, by analyzing a full weighted
graph whose edge weights are determined by (in the simplest case, taken equal
to) the respective pairwise correlation coefficients. The transition from
unweighted to weighted networks calls for the respective extension of graph
metrics definitions. For many metrics such extensions are available Newman
(2016), but may be not unique, thus giving rise to an additional problem of
choosing the best weighted-network modification of a particular graph metric.
An important consideration to guide this choice is that the networks of
interest are correlation networks.

This approach was successfully followed in the neuroscience context Masuda
_et al._ (2018), where an improved formulation of local clustering coefficient
(LCC) specially focused on weighted correlation networks was proposed and
demonstrated to outperform the conventional thresholding approach in terms of
an illustrative neuroscience problem (revealing the age dependence of human
brain network structure based on functional magnetic resonance data). The key
distinctive feature of the special definition of LCC for weighted correlation
networks Masuda  _et al._ (2018) consists in accounting for three-way partial
correlations in order to compensate for spurious correlations caused by
indirect interactions, which otherwise disguise the true interaction structure
in the correlation matrix.

No literature is currently available to implement this approach in climate
science. The present study aims at filling this gap by adapting the
correlations-focused weighted-network LCC definition from Masuda  _et al._
(2018) to climate network analysis. We further modify the LCC definition of
Masuda  _et al._ (2018) by switching from (parametric) Pearson’s to (non-
parametric) Kendall’s correlations (which moves from capturing linear to more
general monotonic dependencies between variables and imposes weaker
requirements to sample size as compared to both Pearson’s and Spearman’s
correlation Bonett and Wright (2000)) and relate the LCC of mean sea level
pressure correlation networks to available data on tropical cyclones. We show
that the LCC modification based on three-way partial correlations outperforms
the unweighted thresholded network LCC as a marker of tropical cyclone.

##  II Methods

###  II.1 Cimate data

As a source of multivariate time series for constructing climate correlation
networks we use ERA5 reanalysis data Hersbach  _et al._ (2020), which is
essentially a model-based interpolation of available climatic observational
data to a regular grid over geographic coordinates and time. We use time
series of mean sea level pressure (MSLP) on a coordinate grid with step of
0.75° over latitude and longitude, and time grid with step of 3 hours. The
data are taken for the sea surface within the range 5°N to 30°N latitude and
50°E to 100°E longitude (northern part of the Indian ocean), the time range
covers 38 years from 1982 to 2019.

We compute pairwise correlations of MSLP anomalies between the nodes of the
geographic coordinate grid. Anomaly is defined as the deviation of the ERA5
reanalysis data from local daily climate normals, which in turn are computed
by averaging the observable (here, MSLP) at the specific location (grid node)
over all daily time samples for the particular date of year and over all years
of observation (for simplicity, data for 29 February of leap years are
discarded).

We relate our findings to the Best Track data Mohapatra  _et al._ (2012) on
tropical cyclones in the region, which are obtained from the Regional
Specialized Meteorological Centre for Tropical Cyclones Over North Indian
Ocean (India) and contain information on the position and strength of all
registered cyclones in the region over time.

###  II.2 Correlation analysis

To quantify correlation, we use the Kendall’s rank correlation coefficient
$\tau$ Kendall (1938) (more precisely, the Goodman-Kruskal version thereof,
see below), which is defined for a joint bivariate sample (taken for
definiteness on a grid of observation times $\\{t_{i}\\}$) of two random
variates $(x_{1}(t_{i}),x_{2}(t_{i}))$ as a normalized difference between the
counts of concordant and discordant pairs of bivariate observations111A pair
of bivariate observations $(x_{1}(t_{i}),x_{2}(t_{i}))$ and
$(x_{1}(t_{j}),x_{2}(t_{j}))$ is called concordant (discordant), if the
product $(x_{1}(t_{i})-x_{1}(t_{j}))(x_{2}(t_{i})-x_{2}(t_{j}))$ is positive
(negative), implying that this particular pair of observations shows an
increasing (decreasing) dependence between $x_{1}$ and $x_{2}$.. Normalization
is meant to ensure that $\tau\in[-1,1]$, and in particular, $\tau=1$
($\tau=-1$) implies a deterministic increasing (decreasing) functional
dependence between $x_{1}$ and $x_{2}$. Different formulations for the
normalizing denominator are available in the literature, depending on the
chosen way to account for ties (i.e. equality cases with
$x_{k}(t_{i})=x_{k}(t_{j})$), which is still a matter of research Amerise and
Tarsitano (2015). As long as ties are negligible in our problem due to the
continuous nature of climate variables, the method of their resolution is not
actually significant. For definiteness, we opt to discard tied pairs if such
occur, and define the normalizing denominator as the total of concordant and
discordant pairs, as suggested in Adler (1957), (Kruskal, 1958, Eq. (15.2) and
below therein), the resultant quantity also known as the Goodman-Kruskal gamma
coefficient (Gibbons and Chakraborti, 2003, Eq. (2.39) therein), named after
authors who applied it to variates taking on finite sets of values (Goodman
and Kruskal, 1954, Eq. (21) therein).

Our choice of the Kendall’s rank correlation coefficient over other available
measures of association is due to the following considerations: (i) it imposes
weaker requirements on sample size compared to both Pearson’s and Spearman’s
correlation coefficients Bonett and Wright (2000), which translates into
shorter length of the sliding time window, and hence better time resolution;
(ii) like any rank-order statistics (including Spearman’s correlation), it is
invariant to monotone nonlinear transformations of variables (Gibbons and
Chakraborti, 2003, Sec. 5.5), which is not the case with Pearson’s correlation
being only invariant to linear transformations (essentially, Kendall’s and
Spearman’s correlation assess the proximity of a bivariate sample to an
arbitrary monotone dependence between variates, instead of proximity to linear
dependence as in Pearson’s correlation); (iii) it is more robust to outliers
(large-amplitude noise) than both Pearson’s and Spearman’s correlation
coefficients Xu  _et al._ (2013).

Correlation coefficient attributed to a particular moment of time is computed
over a time window of 10 days preceding the given moment. We developed an
optimized online algorithm to compute Kendall’s rank correlation coefficient
over a sliding time window, which benefits from the reuse of computation
results obtained at the previous position of the sliding window.

###  II.3 Network metrics

We consider both weighted and unweighted networks, the former represented by
the full unchanged matrix of Kendall’s pairwise rank correlation coefficients
$(\tau_{ij})$, and the latter constructed by thresholding the correlation
matrix: a pair of nodes are considered connected when their correlation
coefficient exceeds a threshold, which in turn is chosen so that the fraction
of connected nodes among the total number of node pairs (edge density) is a
specified quantity (taken equal to 5% for definiteness).

Local clustering coefficient (LCC) for a given node $i$ on an unweighted graph
is defined as the ratio of the number of closed triangles (connected triplets
of nodes) containing the node $i$ to the total number of node triplets
constructed from this node and its adjacent nodes Watts and Strogatz (1998).
LCC may be expressed in terms of the adjacency matrix $a_{ij}$ of an
unweighted graph (by definition, $a_{ij}=1$ if a link is present between the
nodes $i$ and $j$, and $a_{ij}=0$ otherwise) as

$LCC_{i}^{\text{unw}}=\frac{\sum_{1\leq j<l\leq N,(j,l\neq i)}a_{ij}\cdot
a_{il}\cdot a_{jl}}{\sum_{1\leq j<l\leq N(j,l\neq i)}a_{ij}\cdot a_{il}},$ (1)

where $N$ is the total number of nodes in the network, and indices $i$, $j$,
$l$ enumerate the nodes. Note that the denominator in (1) is essentially the
number of unordered pairs (2-combinations) among the nodes adjacent to $i$ and
can be equivalently expressed as $k_{i}(k_{i}-1)/2$, where $k_{i}$ is the
degree of the node $i$ (the number of the nodes adjacent to $i$). We show the
formulation (1) due to its straightforward correspondence with the weighted-
network extension of the LCC definition (see Eq. (2) below).

Several extensions of the LCC definition to account for weighted networks are
available in the literature Masuda  _et al._ (2018); Saramäki  _et al._
(2007); Rubinov and Sporns (2010, 2011); Wang  _et al._ (2017). The version of
Masuda  _et al._ (2018) was specifically focused on correlation networks, as
opposed to other LCC definitions adapted to weighted graphs in general. The
idea behind this definition is to compensate for the impact that indirect
interaction paths (e.g. node $j$ interacting with node $i$, which in turn
interacts with node $l$) inevitably have upon the correlation matrix. Namely,
the nodes $j$ and $l$ in the example above are expected to exhibit some
correlation even in the absence of actual direct interaction between them,
implying that such indirect correlations disguise the true interaction
structure in the correlation matrix. In order to minimize the influence of
such indirect correlations on the resultant LCC value, the definition of LCC
in Masuda  _et al._ (2018) makes use of the three-way partial Pearson’s
product-moment correlation coefficient $\rho^{\text{part}}_{jl|i}$ Whittaker
and Whittaker (1990), which essentially indicates the surplus correlation
between the nodes $j$ and $l$ beyond their indirect correlation through the
node $i$.

In our study we modify the LCC definition from Masuda  _et al._ (2018) by
replacing all Pearson’s product-moment correlation coefficients by the
respective Kendall’s rank correlation coefficients, which produces the
expression

$LCC_{i}^{\text{wei}}=\frac{\sum_{1\leq j<l\leq N,(j,l\neq
i)}|\tau_{ij}\cdot\tau_{il}\cdot\tau^{\text{part}}_{jl|i}|}{\sum_{1\leq
j<l\leq N(j,l\neq i)}|\tau_{ij}\cdot\tau_{il}|},$ (2)

where the three-way partial Kendall’s rank correlation coefficient
$\tau^{\text{part}}_{jl|i}$ is defined according to (Gibbons and Chakraborti,
2003, Section 12.6) as

$\tau^{\text{part}}_{jl|i}=\frac{\tau_{jl}-\tau_{ij}\cdot\tau_{il}}{\sqrt{1-\tau_{ij}^{2}}\sqrt{1-\tau_{il}^{2}}}.$
(3)

Notably, this expression is analogous to the definition of the three-way
partial Pearson’s correlation coefficient $\rho^{\text{part}}_{jl|i}$ Masuda
_et al._ (2018); Whittaker and Whittaker (1990) with all pairwise Pearson’s
correlations replaced by the respective Kendall’s correlations.

As long as the expression (3) contains a division by zero whenever
$|\tau_{ij}=1|$ or $|\tau_{il}=1|$, such values of the indices $j$ and $l$ are
excluded from summation both in the numerator and in the denominator of (2),
whenever this occurs (such cases turn out to be rare but not impossible in our
actual computations). Essentially, as mentioned in Sec. II.2, unity absolute
value of Kendall’s $\tau$ implies that the bivariate statistics within such a
pair of variables (between nodes $i$ and $j$, or $i$ and $l$) is
indistinguishable by available data from a deterministic monotone functional
dependence. This may be seen as a justification for our decision to exclude
such pairs of nodes from summation, which may be interpreted as them being
temporarily lumped together into a single node (for a particular position of
the time window), as soon as they turn out to behave as a single node anyway.
That said, excluding unity-correlated pairs of nodes inevitably introduces a
methodical perturbation into the obtained LCC value, and the final
justification for this operation relies mostly on our observation that such
perturbations, being rare and relatively small, do not affect the capability
of thus computed LCC as a climate marker.

As it is pointed out in Masuda  _et al._ (2018), the use of three-way partial
correlations, which naturally appear in (2) due to the triangle-based nature
of the LCC definition, allows to avoid the computation of full partial
correlation matrices (taking account of all indirect correlation paths between
each pair of nodes), which is a difficult task in terms of both computational
burden and poor tractability for statistical estimation Brier  _et al._
(2015).

Note also that due to the absolute value taken in (2), positive and negative
correlations are assumed to contribute equally to the LCC value Masuda  _et
al._ (2018).

It was shown in Masuda  _et al._ (2018) that LCC defined in terms of three-way
partial correlations better reveals the age dependence of brain connectivity
by functional magnetic resonance data as compared to the LCC definition (1)
for unweighted networks (constructed by thresholding the correlation matrix
with two threshold values corresponding to fixed edge densities 10% and 20%),
and as well to other available LCC definitions for weighted networks Saramäki
_et al._ (2007); Rubinov and Sporns (2010, 2011); Wang  _et al._ (2017). The
present study is the first to apply the LCC definition (2) to climate
networks.

##  III Results

We have computed snapshots of LCC (namely, its respective formulations (1) or
(2)) for each position of the sliding time window using both the full weighted
correlation network and the unweighted network with edge density 5% obtained
by thresholding the correlation matrix. We focus on the relation of the
computed LCC value (both versions thereof) to the available data on tropical
cyclones in the region of interest, as per the Best Track database Mohapatra
_et al._ (2012). We represent the results in the form of a color geographic
map for LCC (sea area only) for each given position of the time window,
superimposed with a track of a cyclone if the current time of observation
(defined as the ending time of the time window used to construct the network)
falls within the existence period of this cyclone.

We aim at demonstrating that the LCC formulation (2) based on a weighted
network may deliver information associated with a tropical cyclone above and
beyond the information delivered by the unweighted LCC formulation (1). In the
present study we limit ourselves to substantiating this statement by examples
where it is evident by a visual side-by-side comparison of the respective
color maps. Several such examples are presented in the Figure 1, where each
row corresponds to a particular moment of observation (indicated above each
map) during a tropical cyclone (the name and the time period of each cyclone
is indicated above the cyclone track), color encodes the unweighted (left
column) and weighted (right column) versions of LCC at the same observation
time. Each cyclone track is marked with a series of circles, whose radius
visualizes the cyclone strength in the respective location. The current
position of a cyclone (at the time of observation) is marked with a red
circle.

The comparison of the respective maps in the Figure 1 on the left and on the
right leads to an observation that the weighted version of LCC (on the right)
produces a more pronounced local peak in the vicinity of the current location
of the cyclone than the unweighted version thereof (on the left). Note that
the Figure 1 intentionally shows only examples supporting the statement above.
In other cases the weighted LCC version may show little or no improvement over
the unweighted version. Nevertheless, when addressing a question, whether a
new quantitative indicator (here, the weighted-network LCC formulation) should
be included in the toolset of climate network analysis to get the most
information of it, even an individual positive example is sufficient to draw a
positive conclusion.

##  IV Conclusion

We have found several particular cases where the new proposed network marker
defined by Eqs. (2) and (3) turns out to be better associated to tropical
cyclones than the unweighted local clustering coefficient. This finding alone
is sufficient to recommend the usage of this marker among the toolset of
climate network analysis, in particular, in application to predicting extreme
climate events. Up to this point we have not identified any specific
predictive criterion (e.g. for tropical cyclones) based on the proposed
marker. Formulating such a criterion, which remains a matter of further study,
would additionally make possible a quantitative comparison of the considered
network measures in terms of prediction accuracy.

Figure 1: Color maps of local clustering coefficient calculated from
unweighted (left) and weighted (right) correlation networks of mean sea level
pressure anomaly during several selected tropical cyclones in the northern
Indian Ocean. Each row corresponds to a particular position of the sliding
time window used to construct the networks. Observation time (i.e. the ending
time of the window) is indicated above the map. Tracks of tropical cyclones
are marked with a series of circles, whose radius denotes the cyclone
strength. Name and time period of each cyclone is indicated above the track.
Current position of each cyclone at the time of observation is marked with a
red circle.

## References

  * Langfelder and Horvath (2008) P. Langfelder and S. Horvath, BMC bioinformatics 9, 1 (2008).
  * Bruggeman and Westerhoff (2007) F. J. Bruggeman and H. V. Westerhoff, TRENDS in Microbiology 15, 45 (2007).
  * Van der Greef  _et al._ (2007) J. Van der Greef, S. Martin, P. Juhasz, A. Adourian, T. Plasterer, E. Verheij, and R. McBurney, Journal of proteome research 6, 1540 (2007).
  * Friedman and Alm (2012) J. Friedman and E. J. Alm, PLoS Comput Biol 8, e1002687 (2012).
  * Rosato  _et al._ (2018) A. Rosato, L. Tenori, M. Cascante, P. R. D. A. Carulla, V. A. M. Dos Santos, and E. Saccenti, Metabolomics 14, 1 (2018).
  * Fransson and Marrelec (2008) P. Fransson and G. Marrelec, Neuroimage 42, 1178 (2008).
  * Bassett and Sporns (2017) D. S. Bassett and O. Sporns, Nature neuroscience 20, 353 (2017).
  * Tsonis  _et al._ (2006) A. A. Tsonis, K. L. Swanson, and P. J. Roebber, Bulletin of the American Meteorological Society 87, 585 (2006).
  * Donges  _et al._ (2009) J. F. Donges, Y. Zou, N. Marwan, and J. Kurths, The European Physical Journal Special Topics 174, 157 (2009).
  * Gozolchiani  _et al._ (2011) A. Gozolchiani, S. Havlin, and K. Yamasaki, Physical review letters 107, 148501 (2011).
  * Dijkstra  _et al._ (2019a) H. A. Dijkstra, E. Hernández-García, C. Masoller, and M. Barreiro,  _Networks in Climate_ (Cambridge University Press, 2019).
  * Feng  _et al._ (2016) Q. Y. Feng, R. Vasile, M. Segond, A. Gozolchiani, Y. Wang, M. Abel, S. Havlin, A. Bunde, and H. A. Dijkstra, Geoscientific Model Development (2016).
  * Dijkstra  _et al._ (2019b) H. Dijkstra, E. Hernandez-Garcia, C. Lopez, _et al._ , Frontiers in Physics 7, 153 (2019b).
  * Bockmayr  _et al._ (2013) M. Bockmayr, F. Klauschen, B. Györffy, C. Denkert, and J. Budczies, BMC systems biology 7, 1 (2013).
  * Newman (2016) M. E. J. Newman, in  _The New Palgrave Dictionary of Economics_ (Palgrave Macmillan UK, London, 2016) pp. 1–8.
  * Masuda  _et al._ (2018) N. Masuda, M. Sakaki, T. Ezaki, and T. Watanabe, Frontiers in Neuroinformatics 12, 7 (2018).
  * Bonett and Wright (2000) D. G. Bonett and T. A. Wright, Psychometrika 65, 23 (2000).
  * Hersbach  _et al._ (2020) H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, R. Radu, D. Schepers, _et al._ , Quarterly Journal of the Royal Meteorological Society 146, 1999 (2020).
  * Mohapatra  _et al._ (2012) M. Mohapatra, B. Bandyopadhyay, and A. Tyagi, Natural Hazards 63, 1285 (2012).
  * Kendall (1938) M. G. Kendall, Biometrika 30, 81 (1938).
  * Note (1) A pair of bivariate observations $(x_{1}(t_{i}),x_{2}(t_{i}))$ and $(x_{1}(t_{j}),x_{2}(t_{j}))$ is called concordant (discordant), if the product $(x_{1}(t_{i})-x_{1}(t_{j}))(x_{2}(t_{i})-x_{2}(t_{j}))$ is positive (negative), implying that this particular pair of observations shows an increasing (decreasing) dependence between $x_{1}$ and $x_{2}$.
  * Amerise and Tarsitano (2015) I. L. Amerise and A. Tarsitano, Journal of Applied Statistics 42, 2584 (2015).
  * Adler (1957) L. M. Adler, Journal of the American Statistical Association 52, 33 (1957).
  * Kruskal (1958) W. H. Kruskal, Journal of the American Statistical Association 53, 814 (1958).
  * Gibbons and Chakraborti (2003) J. Gibbons and S. Chakraborti,  _Nonparametric Statistical Inference: Fourth Edition, Revised and Expanded. Statistics: A Series of Textbooks and Monographs_ (New York, NY: Marcel Dekker, Inc., 2003).
  * Goodman and Kruskal (1954) L. A. Goodman and W. H. Kruskal, Journal of the American Statistical Association 49, 732 (1954).
  * Xu  _et al._ (2013) W. Xu, Y. Hou, Y. Hung, and Y. Zou, Signal Processing 93, 261 (2013).
  * Watts and Strogatz (1998) D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998).
  * Saramäki  _et al._ (2007) J. Saramäki, M. Kivelä, J.-P. Onnela, K. Kaski, and J. Kertesz, Physical Review E 75, 027105 (2007).
  * Rubinov and Sporns (2010) M. Rubinov and O. Sporns, Neuroimage 52, 1059 (2010).
  * Rubinov and Sporns (2011) M. Rubinov and O. Sporns, Neuroimage 56, 2068 (2011).
  * Wang  _et al._ (2017) Y. Wang, E. Ghumare, R. Vandenberghe, and P. Dupont, Neural computation 29, 313 (2017).
  * Whittaker and Whittaker (1990) J. Whittaker and J. Whittaker,  _Graphical models in applied multivariate statistics_ , Vol. 19 (Wiley Chichester, 1990).
  * Brier  _et al._ (2015) M. R. Brier, A. Mitra, J. E. McCarthy, B. M. Ances, and A. Z. Snyder, NeuroImage 121, 29 (2015).