./images/

# Multiqubit noise deconvolution and characterization

Simone Roncallo  0000-0003-3506-9027<EMAIL_ADDRESS>Dipartimento di Fisica, Università degli Studi di Pavia, Via Agostino Bassi 6,
I-27100, Pavia, Italy  INFN Sezione di Pavia, Via Agostino Bassi 6, I-27100,
Pavia, Italy  Lorenzo Maccone  0000-0002-6729-5312<EMAIL_ADDRESS>Dipartimento di Fisica, Università degli Studi di Pavia, Via Agostino Bassi 6,
I-27100, Pavia, Italy  INFN Sezione di Pavia, Via Agostino Bassi 6, I-27100,
Pavia, Italy  Chiara Macchiavello  0000-0002-2955-8759
<EMAIL_ADDRESS>Dipartimento di Fisica, Università degli Studi di
Pavia, Via Agostino Bassi 6, I-27100, Pavia, Italy  INFN Sezione di Pavia, Via
Agostino Bassi 6, I-27100, Pavia, Italy

###### Abstract

We present a noise deconvolution technique for obtaining noiseless expectation
values of noisy observables at the output of multiqubit quantum channels. For
any number of qubits or in the presence of correlations, our protocol applies
to any mathematically invertible noise model, even when its inverse map is not
physically implementable, i.e. when it is neither completely-positive nor
trace-preserving. For a generic observable affected by Pauli noise it provides
a quadratic speedup, always producing a rescaling of its Pauli basis
components. We show that it is still possible to achieve the deconvolution
while experimentally estimating the noise parameters, whenever these are
unknown (bypassing resource-heavy techniques such as quantum process
tomography). We provide a simulation, with examples for both Pauli and non-
Pauli channels.

Noise deconvolution; Noise mitigation; Multiqubit quantum channel; Multiqubit
correlated noise; Pauli channel; Pauli transfer matrix; Noise
characterization;

##  I Introduction

Noise in quantum systems can affect the measurement outcome of any observable,
modifying the results of measurement-based protocols or procedures, such as
state and process tomography [1, 2, 3, 4] or quantum simulator experiments
[5]. This has also important consequences in quantum computation: noise
sensitivity remains one of the main drawbacks that prevents quantum computers
from outperforming their classical counterparts. Several noise mitigation
techniques have been considered in the literature, aimed to reduce errors and
potential loss of data in the computation process [6, 7, 8]. Recently, a noise
deconvolution technique was illustrated for observables of single-qubit
systems [9], by means of a tomographic reconstruction formula that acts like a
postprocessing operation on the noisy data, without introducing modifications
to the system.

In this paper we discuss a deconvolution technique that applies to any
multiqubit (possibly correlated) noise model [10, 11, 12, 13, 14], provided
its inverse map exists, even if not physically implementable [15]. We modify
the point of view of [9], adopting an operational framework that is more
suitable for multiqubit implementations. Unlike other approaches in the
literature [6, 7], our protocol is passively implemented at the data
processing stage: it does not require further experimental (or circuital)
configurations, nor active modifications of the original system. For this
reason, its range of applicability is not limited to quantum computing: it
applies to generic quantum measurements on noisy states, e.g. it can be used
to reverse open quantum dynamics [16, 17].

Figure 1: (a) Expectation value of a multiqubit observable $O$, affected by
Pauli noise $\mathcal{N}$. (b) The same scheme is applied to the non-zero
components of $O$. Deconvolution is achieved by rescaling the noisy data with
the corresponding diagonal entry of the Pauli transfer matrix
$\Gamma_{\mathcal{N}}$.

Our method works in general, e.g. it does not require Markovian correlations.
In the specific but important case of Pauli channels [18], it provides a
quadratic speedup over other reconstructions. Given an observable $O$, we show
that its noiseless expectation value can be tomographically reconstructed by
performing local measurements on those elements of the Pauli basis for which
$O$ takes non-zero components and by rescaling them in terms of the
corresponding entries of the Pauli transfer matrix (PTM) of the channel [19,
20, 21]. This procedure works for any number of qubits and does not require
the complete inversion of the noise map, nor the complete calculation of its
PTM. For $n$-qubit Pauli channels the computational complexity of our
procedure, i.e. the number of factors required to complete the deconvolution,
scales with the number $r=1,2,...,d^{2}$ of non-zero components of $O$ in the
Pauli basis, with $d=2^{n}$ the Hilbert space dimension ($r=1$ represents an
observable made by exactly one basis element, while $r=d^{2}$ represents the
worst cases in which $O$ has all non-zero components). This scenario is
summarized in Fig. 1. Then, we discuss a characterization of the noise map
that provides the necessary PTM entries in term of a few measurements on the
Pauli basis, without running a full process tomography of the system. This
represents an efficient alternative whenever the theoretical computation of
the PTM is not doable, e.g. when the noise parameters are unknown. Finally, we
show that our protocols apply also to non-Pauli channels, with reduced
efficiency. In this case, we provide the deconvolution from the complete
inversion of the PTM, and we show that it scales as $d^{4}$ in the worst case.
So, in essence, our procedure experiences a quadratic gain ($d^{2}$ vs
$d^{4}$) for Pauli channels.

In Section II we review the vectorization of operators and channels. In
Section III we present noise deconvolution for multiqubit Pauli channels,
addressing their characterization whenever the noise parameters are unknown.
Then, we generalize our discussion to the non-Pauli case. In Section IV we
consider the explicit deconvolution of $n$-qubit bit-flip and depolarizing
correlated noises, simulating the latter for $n=3$. As a non-Pauli example, we
consider a two-qubit amplitude damping correlated channel, which models qubits
losses inside devices.

##  II Vectorization

We consider the Hilbert space of an $n$-qubit system. The basis for the set of
operators is

$\big{\\{}\sigma_{\alpha_{1}}\otimes\sigma_{\alpha_{2}}\otimes...\otimes\sigma_{\alpha_{n}}\
|\ \alpha_{1},\alpha_{2},...,\alpha_{n}=0,1,2,3\big{\\}}\ ,$ (1)

with $\sigma_{0}=\mathds{1}_{2}$, $\sigma_{1}=\sigma_{x}$,
$\sigma_{2}=\sigma_{y}$ and $\sigma_{3}=\sigma_{z}$. We write the Pauli basis
in the notation

$\big{\\{}\mathcal{P}_{k}\ |\ k=0,1,2,3,...,d^{2}-1\big{\\}}\ ,$ (2)

with $d=2^{n}$ and $\mathcal{P}_{k}$ denoting the generic element of Eq. 1 in
lexicographic order. We introduce the vectorized representation [19], in which
each element of the basis $\mathcal{P}_{k}$ is mapped to a vector
$\ket{k}\rangle$. In this space, an operator $A$ is represented as a $1\times
d^{2}$ column vector

$\ket{A}\rangle=\sum_{k=0}^{d^{2}-1}A_{k}\ket{k}\rangle\ ,$ (3)

with $A_{k}=\langle\bra{k}\ket{A}\rangle$ given by the Hilbert-Schmidt inner
product

$\langle\bra{A}\ket{B}\rangle:=\frac{1}{d}\Tr[A^{\dagger}B]\ ,$ (4)

with $B$ any operator on this Hilbert space.

Consider a quantum system in the state $\rho$. A quantum channel is a linear
completely-positive and trace-preserving (CPTP) map that modifies the system
state as $\rho\to\Phi(\rho)$. In the vectorized framework, a channel $\Phi$ is
represented by a $d^{2}\times d^{2}$ matrix

$\Gamma_{\Phi}=\sum_{j,q=0}^{d^{2}-1}\Gamma_{jq}\ket{j}\rangle\langle\bra{q}\
,$ (5)

with components given by the Hilbert-Schmidt inner product

$\Gamma_{jq}=\langle\bra{j}\Gamma_{\Phi}\ket{q}\rangle=\frac{1}{d}\Tr[\mathcal{P}_{j}\Phi(\mathcal{P}_{q})]\
.$ (6)

This is called the Pauli transfer matrix of the channel [19]. In this
representation, the action of the channel $\Phi(A)$ is given by a matrix-
vector multiplication

$\ket{\Phi(A)}\rangle=\Gamma_{\Phi}\ket{A}\rangle=\sum_{j,q=0}^{d^{2}-1}\Gamma_{jq}A_{q}\ket{j}\rangle\
.$ (7)

The CPTP condition guarantees that $\Gamma_{0q}=\delta_{0q}$, with
$\delta_{jq}$ the Kronecker delta. For unital channels, i.e. when
$\Phi(\mathds{1})=\mathds{1}$, it holds also that $\Gamma_{j0}=\delta_{j0}$.

Consider a channel $\Phi$, its adjoint $\Phi^{*}$ is the map satisfying

$\langle\bra{A}\ket{\Phi^{*}(B)}\rangle=\langle\bra{\Phi(A)}\ket{B}\rangle\ .$
(8)

This implies that the adjoint PTM, here denoted by $\Gamma_{\Phi}^{*}$, is
precisely the Hermitian conjugate of $\Gamma_{\Phi}$.

##  III Noise deconvolution

Noise in open quantum systems is modeled in terms of quantum channels [1],
namely, linear CPTP operations $\mathcal{N}$ that map the ideal, i.e.
noiseless, state $\rho$ into a noisy one $\rho^{\prime}=\mathcal{N}(\rho)$.
Different choices of $\mathcal{N}$ correspond to different noise models, e.g.
the bit-flip, the dephasing, the depolarizing, or the amplitude damping noises
[1]. With the state modified by $\mathcal{N}$, any measurement performed on an
observable $O$ becomes noisy, namely its expectation value $\langle
O\rangle_{\rho}$ is mapped to $\langle O\rangle_{\rho^{\prime}}$. In this
section we introduce noise deconvolution as a technique that provides the
ideal expectation value $\langle O\rangle_{\rho}$ of arbitrary operators,
using the noisy data obtained on $\rho^{\prime}$. Our protocol applies to any
mathematically invertible noise model, i.e. one for which
$\Gamma_{\mathcal{N}}^{-1}$ exists, even when the inverse channel is not
physically implementable, i.e. when $\mathcal{N}^{-1}$ is not CPTP.

We consider the noise deconvolution equation derived in [9]

$\left\langle
O\right\rangle_{\rho}=\left\langle\mathcal{N}^{*^{-1}}(O)\right\rangle_{\rho^{\prime}}\
,$ (9)

which yields the ideal, i.e. noiseless, expectation value of an arbitrary
observable $\left\langle O\right\rangle_{\rho}$ (or even of a non-observable
operator) by evaluating, instead, the inverted adjoint map
$\mathcal{N}^{*^{-1}}(O)$ over the noisy state
$\rho^{\prime}=\mathcal{N}(\rho)$. In other proposals [22], similar inversions
are implemented physically (often only approximately) by introducing suitable
modifications to the channel. Here, instead, we just use the noisy measured
data to calculate the noiseless value. In other words and in contrast to
previous proposals, our noise-inversion reconstruction is implemented entirely
and solely at the data processing stage.

In the vectorized framework, Eq. 9 reads

$\langle\bra{\rho}\ket{O}\rangle=\langle\bra{\rho^{\prime}}\Gamma_{\mathcal{N}}^{*^{-1}}\ket{O}\rangle\
,$ (10)

where $\Gamma_{\mathcal{N}}^{*^{-1}}$ is the inverse adjoint PTM.

We start with the deconvolution of multiqubit Pauli channels, for which the
vectorization guarantees a quadratic speedup in efficiency over the general
case (treated below). The current method generalizes and supersedes the
single-qubit analysis presented in [9]. In this case we derive a
reconstruction formula without completely inverting the noise channel, instead
using only some components of the PTM. We first apply this procedure to an
observable made by one of the possible $n$-fold tensor products of the Pauli
matrices, i.e. one of the elements of the basis, then we extend our
considerations to the expectation value of a generic observable that takes all
the contributions from the Pauli basis. At the end of the section, we
generalize our discussion to the non-Pauli case.

In the Pauli basis, the Kraus representation [1] of an $n$-qubit Pauli channel
is [18]

$\mathcal{N}(\rho)=\sum_{j=0}^{d^{2}-1}\beta_{j}\mathcal{P}_{j}\rho\mathcal{P}_{j}\
,$ (11)

with $\sum_{j}\beta_{j}=1$ and $\beta_{j}\geq 0\ \forall j$. Such channels
represent an example of random unitary maps [23], where each unitary
$\mathcal{P}_{j}$ is applied with probability $\beta_{j}$. In the case of Eq.
11, this produces a depolarizing contraction of the Bloch hypersphere, whose
intensity and direction depend on the choice of $\beta_{j}$ [24].

A Pauli channel yields a diagonal PTM

$\Gamma_{\mathcal{N}}=\sum_{j=0}^{d^{2}-1}\lambda_{j}\ket{j}\rangle\langle\bra{j}\
,$ (12)

where $\lambda_{j}=\langle\bra{j}\Gamma_{\mathcal{N}}\ket{j}\rangle$. As the
only requirement, we ask $\Gamma_{\mathcal{N}}$ to be mathematically
invertible, i.e. that $\lambda_{j}\neq 0\ \forall j$.111An example of a non-
invertible map is represented by a depolarizing channel with $p=1$, for which
the Bloch hypersphere collapses to a single point. For now, we also assume
that the noise parameters are known so that a theoretical computation of the
PTM is always possible. This last assumption simplifies our analysis, but it
is not necessary. We discuss the case of unknown noise in the following
paragraphs.

Since all the Kraus operators in Eq. 11 are Hermitian, the adjoint of the
Pauli channel is the channel itself [9], yielding
$\Gamma_{\mathcal{N}}^{*^{-1}}=\Gamma_{\mathcal{N}}^{-1}$. Then the inverse
PTM is diagonal too and its components read $1/\lambda_{j}$.

First, consider an observable given by the $k$th element of the basis, which
in vectorized notation corresponds to a column vector with only one non-zero
component $\ket{O}\rangle=\ket{k}\rangle$. The action of the inverse PTM
yields

$\Gamma^{-1}_{\mathcal{N}}\ket{O}\rangle=\Gamma_{kk}^{-1}\ket{k}\rangle=\frac{1}{\lambda_{k}}\ket{k}\rangle\
.$ (13)

Using Eq. 10, the expectation value follows as
$\langle\bra{\rho}\ket{k}\rangle=\lambda_{k}^{-1}\langle\bra{\rho^{\prime}}\ket{k}\rangle$,
which in standard non-vectorized notation reads

$\langle\mathcal{P}_{k}\rangle_{\rho}=\frac{d}{\Tr[\mathcal{P}_{k}\mathcal{N}(\mathcal{P}_{k})]}\langle\mathcal{P}_{k}\rangle_{\rho^{\prime}}\
.$ (14)

This shows that for an $n$-qubit Pauli channel, the deconvolution of the
expectation value of the $k$th element of the basis is always obtained as a
rescaling of the noisy outcome, which depends on the $k$th element on the
diagonal of the inverse PTM.

We now discuss the deconvolution of a generic observable $O$ subject to an
arbitrary $n$-qubit Pauli channel. We expand $O$ in terms of the basis vectors
$\ket{O}\rangle=\sum_{k}O_{k}\ket{k}\rangle$ to compute the right-hand side of
Eq. 10. Applying the same strategy, each component $O_{k}$ must be rescaled by
the corresponding element on the diagonal of the inverse PTM. Then, the
reconstructed measurement outcome reads

$\langle
O\rangle_{\rho}=\sum_{k=0}^{d^{2}-1}\frac{d}{\Tr[\mathcal{P}_{k}\mathcal{N}(\mathcal{P}_{k})]}O_{k}\langle\mathcal{P}_{k}\rangle_{\rho^{\prime}}\
.$ (15)

Namely, we obtain the ideal noiseless outcome (over $\rho$) by processing the
noisy expectation values (over $\rho^{\prime}=\mathcal{N}(\rho)$) of those
Pauli basis elements $\mathcal{P}_{k}$ that contribute to the expansion of
$O$.

This procedure works for any number $n$ of (even correlated) qubits, and it
involves only local measurements on each element of the basis. The
computational complexity, i.e. the number of factors required, scales with the
number $r=1,2,...,d^{2}$ of non-zero elements of the expansion of $O$ on the
Pauli basis. In the trivial case in which $O$ is exactly an $n$-qubit Pauli
matrix, the deconvolution always requires a single measurement and one PTM
entry, for any number of qubits. On the other hand, when $O$ is a generic
observable, the deconvolution requires a measurement on the entire basis and
the computation of all the $d^{2}$ diagonal components of the PTM. In any
case, this considerably reduces the number of computations of the PTM entries
over the inversion-based method of [9], which for $n$ qubits always requires
$d^{4}$ operations.

So far we have considered $n$-qubit Pauli channels whose parameters are known
_a priori_. Although it simplifies the analysis, this assumption is not
necessary: we can still estimate these parameters without running a full
process tomography of the channel [1, 3], which is intractable for large $n$.
For a general $n$-qubit Pauli channel $\mathcal{N}$, with unknown Kraus
representation coefficients in Eq. 11, the deconvolution is achieved by means
of the following characterization. Prepare the system in the state

$\rho_{k}=\frac{\mathds{1}+\mathcal{P}_{k}}{d}\quad\text{for }k\neq 0\ ,$ (16)

with $\mathds{1}$ the $n$-qubit identity operator (see Appendix A for a
discussion on positivity). Assuming that $\rho_{k}$ evolves to
$\rho_{k}^{\prime}=\mathcal{N}(\rho_{k})$, i.e. that the channel acts
independently of the preparation scheme of Eq. 16, unitality guarantees that

$\rho_{k}^{\prime}=\frac{\mathds{1}+\mathcal{N}(\mathcal{P}_{k})}{d}\ .$ (17)

The diagonal entries of the PTM then read

$\displaystyle\Gamma_{kk}=\langle\mathcal{P}_{k}\rangle_{\rho_{k}^{\prime}}\quad\text{for
}k\neq 0\ ,$ (18)

with
$\langle\mathcal{P}_{k}\rangle_{\rho_{k}^{\prime}}=\Tr[\mathcal{P}_{k}\mathcal{N}(\rho_{k})]$,
which can be used in the deconvolution formula in Eq. 14 or Eq. 15. This means
that, even without knowing the noise parameters, we can still obtain the
reconstruction factors by measuring the $k$th element of the Pauli basis over
a noisy state initially prepared as in Eq. 16. Again, this requires at most
$r$ operations, with $r=1,2,...,d^{2}$ the number of non-zero components of
$O$ in the Pauli basis. We summarize the entire procedure in Protocol 1.

Protocol 1 Multiqubit Pauli channel deconvolution

Input: Observable $O$ Optional: Pauli channel $\mathcal{N}$
Result: Noiseless expectation value $\langle O\rangle_{\rho}$

1:construction of the basis $\triangleright$ $\\{\mathcal{P}_{k}\\}$ for
$k=0,1,...,d^{2}-1$

2:projection of $O$ on the basis $\triangleright$
$O_{k}=\Tr[O\mathcal{P}_{k}]/d$

3:for $k$ such that $O_{k}\neq 0$ do

4: noisy measurement on the basis $\triangleright$
$\langle\mathcal{P}_{k}\rangle_{\mathcal{N}(\rho)}$

5: if $\mathcal{N}$ is known then

6: get PTM $\triangleright$
$\Gamma_{kk}\leftarrow\Tr[\mathcal{P}_{k}\mathcal{N}(\mathcal{P}_{k})]/d$

7: else

8: state preparation $\triangleright$
$\rho_{k}=(\mathds{1}+\mathcal{P}_{k})/d$

9: PTM characterization $\triangleright$
$\Gamma_{kk}\leftarrow\langle\mathcal{P}_{k}\rangle_{\mathcal{N}(\rho_{k})}$

10:noise deconvolution $\triangleright$ $\langle
O\rangle_{\rho}\leftarrow\sum_{k}\Gamma_{kk}^{-1}O_{k}\langle\mathcal{P}_{k}\rangle_{\mathcal{N}(\rho)}$

Before considering the applications, we discuss the case of noise channels
that cannot be expressed in the form of Eq. 11, i.e. whose Kraus
representation contains at least one operator that is different from a Pauli
$\mathcal{P}_{k}$. These channels describe purely quantum mechanical
processes, e.g. a “spontaneous emission” which is modeled by the amplitude
damping channel [1].

Non-Pauli channels are not self-adjoint and have a non-diagonal PTM, so
$\Gamma_{\mathcal{N}}^{*^{-1}}$ cannot be simply obtained as in Eq. 13. In
this case, the tomographic reconstruction formula of a generic observable $O$
reads

$\langle
O\rangle_{\rho}=\sum_{j,q=0}^{d^{2}-1}\big{(}\Gamma^{*^{-1}}\big{)}_{jq}O_{q}\langle\mathcal{P}_{j}\rangle_{\rho^{\prime}}\
,$ (19)

which, in the case of $O=\mathcal{P}_{k}$ ($O=\ket{k}\rangle$ in the
vectorized notation), reduces to

$\langle\mathcal{P}_{k}\rangle_{\rho}=\sum_{j=0}^{d^{2}-1}\big{(}\Gamma^{*^{-1}}\big{)}_{jk}\langle\mathcal{P}_{j}\rangle_{\rho^{\prime}}\
.$ (20)

In this case, the deconvolution procedure still works, although less
efficiently: it requires more computations and the complete calculation of the
inverse adjoint PTM $\Gamma_{\mathcal{N}}^{*^{-1}}$. For Pauli channels, both
equations consistently reduce to Eq. 14 and Eq. 15. In the next section we
consider the explicit deconvolution of a two-qubit amplitude damping
correlated channel.

When noise is not described by a Pauli channel, as it often occurs in real
devices, the quadratic speedup of Eq. 15 is not achievable. However, generic
noise models are usually split in terms of Pauli and non-Pauli (but simpler)
contributions. See for example [9], where the decoherence of Rigetti Aspen-9
is modeled as a composition of a dephasing (Pauli) and an amplitude damping
(non-Pauli) channel. In the vectorized representation, the composition of
multiple channels reduces to the product of their PTM. Even in this case,
their deconvolution can be separately treated, with still a quadratic
advantage on the Pauli terms, which then act as a rescaling of the non-Pauli
rows or columns. Alternatively, twirling techniques can be employed to remove
the off-diagonal elements of the PTM [18, 25], thus mapping the original
channel to an effective Pauli one.

Similarly to the diagonal case, a straightforward characterization can be
employed whenever the channel is unital and its parameters are unknown.
Prepare the system in the state $\rho_{k}$ of Eq. 16 evolved to Eq. 17 and
then iterate Eq. 18 for all the elements $\mathcal{P}_{j}$ of the Pauli basis,
yielding the PTM as

$\displaystyle\Gamma_{jk}=\langle\mathcal{P}_{j}\rangle_{\rho_{k}^{\prime}}\quad\text{for
}j,k\neq 0\ ,$ (21)

which can be inverted and substituted in Eq. 19 or Eq. 20.

Both Eq. 18 and Eq. 21 provide a direct tomographic reconstruction of the
channel PTM, as long as the channel is unital. This represents an alternative
to the standard approach that recovers the channel Kraus representation
through standard quantum process tomography (a generalization of this
procedure to non-unital channels and a comparison with quantum process
tomography will be discussed in [26]).

(a)  (b)

Figure 2: (a) Simulated noise effect on $\langle\sigma_{z}^{\otimes 3}\rangle$
for a three-qubit depolarizing correlated channel, where the degree of
correlation is parametrized by $\mu\in[0,1]$ and with probability $q=0.00052$.
The system is initially prepared in the noiseless state
$\rho=\ket{000}\\!\bra{000}$. The ideal expectation value
$\langle\sigma_{z}^{\otimes 3}\rangle_{\rho}=1$ is represented by a dashed
horizontal line. The simulation is performed using Qiskit Aer for $8192$
shots, with respect to the number $m$ of subsequent applications of the noise
map $\mathcal{N}^{m}(\rho)$. The noise effect is plotted for different values
of $\mu$. Note that when the value of $\mu$ decreases, the correlation of
qubits increases and the effect of the noise decreases. (b) Simulated
deconvolution of the noise channel for $\mu=0.25$. We compare the noisy output
with the noiseless one, recovered with the noise deconvolution procedure
presented here.

##  IV Applications

In this section we discuss noise deconvolution of several examples of Pauli
and non-Pauli noise models. To analyze correlations we demonstrate our general
method on a specific class of channels, in which a parameter $\mu$ measures
the amount of correlations [10, 11, 12, 27, 28].

We start with the class of Pauli correlated channels defined in [10, 11, 12],
whose Kraus representation reads

$\mathcal{N}(\rho)=\sum_{\\{\alpha_{i}\\}=0}^{3}p_{\alpha_{1}\alpha_{2}...\alpha_{n}}A_{\alpha_{1}\alpha_{2}...\alpha_{n}}\rho
A^{\dagger}_{\alpha_{1}\alpha_{2}...\alpha_{n}}\ ,$ (22)

with $i\in\\{1,2,...,n\\}$ and Kraus operators

$A_{\alpha_{1}\alpha_{2}...\alpha_{n}}=\sigma_{\alpha_{1}}\otimes\sigma_{\alpha_{2}}\otimes...\otimes\sigma_{\alpha_{n}}\
,$ (23)

with
$p_{\alpha_{1}\alpha_{2}...\alpha_{n}}=p_{\alpha_{1}}p_{\alpha_{2}|\alpha_{1}}...p_{\alpha_{n}|\alpha_{n-1}}$
given by the Markov chain [12, 29]

$\displaystyle
p_{\alpha_{j}|\alpha_{i}}=(1-\mu)p_{\alpha_{j}}+\mu\delta_{\alpha_{i}\alpha_{j}}\
,$ (24) $\displaystyle\vec{p}=[1-p,p_{x},p_{y},p_{z}]^{T}\ ,$ (25)

with $p=p_{x}+p_{y}+p_{z}$. This kind of correlation is not a requirement of
our deconvolution technique, it provides only an example to test our protocol.
The parameter $\mu\in[0,1]$ represents the degree of correlation between
couple of qubits. On the one hand, $\mu=0$ represents a memoryless channel,
which is when the qubits are completely uncorrelated. On the other hand,
$\mu=1$ describes a full-memory channel, i.e. when the qubits are completely
correlated. The Pauli channels generalize noise models such as the bit-flip,
the bit-phase-flip, the dephasing and the depolarizing which can be re-
obtained from Eq. 22 by a proper choice of $\vec{p}$. For example, we obtain
the bit-flip channel when $\vec{p}=[1-p,p,0,0]^{T}$ or the depolarizing
channel when $\vec{p}=[1-p,p/3,p/3,p/3]^{T}$.

As the first example, consider the $n$-qubit observable $\sigma_{z}^{\otimes
n}$, whose expectation value is affected by a bit-flip correlated noise [1,
12], i.e. with $\vec{p}=[1-p,p,0,0]^{T}$. In this case, the deconvolution
formula yields

$\langle\sigma_{z}^{\otimes
n}\rangle_{\rho}=f_{n}(p,\mu)\langle\sigma_{z}^{\otimes
n}\rangle_{\rho^{\prime}}\ ,$ (26)

where the reconstruction factor is given by

$f_{n}(p,\mu)=\frac{2^{n}}{\Tr[\sigma_{z}^{\otimes
n}\mathcal{N}(\sigma_{z}^{\otimes n})]}\ .$ (27)

By direct computation, we obtain

$\displaystyle f_{1}(p,\mu)=\frac{1}{1-2p}\ ,$ (28) $\displaystyle
f_{2}(p,\mu)=\frac{1}{1+4(\mu-1)(1-p)p}\ ,$ (29) $\displaystyle
f_{3}(p,\mu)=\frac{1}{(1-2p)[1+4(\mu-1)^{2}(p-1)p]}\ .$ (30)

For a depolarizing correlated channel [10, 12], i.e. when
$\vec{p}=[1-3q/4,q/4,q/4,q/4]^{T}$, the same computation gives

$\displaystyle f_{1}(q,\mu)=\frac{1}{1-q}\ ,$ (31) $\displaystyle
f_{2}(q,\mu)=\frac{1}{1+(\mu-1)(2-q)q}\ ,$ (32) $\displaystyle
f_{3}(q,\mu)=\frac{1}{(1-q)[1+(\mu-1)^{2}(q-2)q]}\ .$ (33)

When the qubits are completely uncorrelated ($\mu=0$), the noise deconvolution
factorizes in terms of the single-qubit contributions. For qubits that are
completely correlated ($\mu=1$), the noise has no effect when $n$ is even,
while it is corrected by a single-qubit contribution when $n$ is odd. In Fig.
2a we plot the effect of a three-qubit depolarizing correlated channel on
$\langle\sigma_{z}^{\otimes 3}\rangle$ with respect to the number of
subsequent applications of the noise map, parametrized in terms of $\mu$. In
Fig. 2b we perform a simulation of the deconvolution, which successfully
reproduces the noiseless expectation value.

As a non-Pauli example, we consider a two-qubit amplitude damping correlated
channel [27, 28]. For a single-qubit system, the amplitude damping channel is
defined in terms of two Kraus operators

$E_{0}=\begin{pmatrix}1&0\\\ 0&\sqrt{\eta}\end{pmatrix}\
E_{1}=\begin{pmatrix}0&\sqrt{1-\eta}\\\ 0&0\end{pmatrix},$ (34)

where $1-\eta$ represents the probability of the system losing a qubit, e.g.
by emitting a photon [1], and it plays the role of channel transmissivity,
e.g. in the case of optical fibers [28].

The two-qubit amplitude damping correlated channel can be obtained as a convex
combination of a memoryless amplitude damping channel $\mathcal{N}_{0}$ with a
memoryful one $\mathcal{N}_{1}$

$\mathcal{N}(\rho)=(1-\mu)\mathcal{N}_{0}(\rho)+\mu\mathcal{N}_{1}(\rho)\ ,$
(35)

where

$\mathcal{N}_{0}(\rho)=\sum_{j=0}^{3}A_{j}\rho A_{j}^{\dagger}\
,\quad\mathcal{N}_{1}(\rho)=\sum_{j=0}^{1}B_{j}\rho B_{j}^{\dagger}\ ,$ (36)

with Kraus operators $A_{0}=E_{0}\otimes E_{0}$, $A_{1}=E_{0}\otimes E_{1}$,
$A_{2}=E_{1}\otimes E_{0}$, $A_{3}=E_{1}\otimes E_{1}$ and

$B_{0}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&\sqrt{\eta}\end{pmatrix}\ B_{1}=\begin{pmatrix}0&0&0&\sqrt{1-\eta}\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\end{pmatrix}.$ (37)

Consider the two-qubit observables $\sigma_{x}^{\otimes 2}$,
$\sigma_{y}^{\otimes 2}$ and $\sigma_{z}^{\otimes 2}$ with noisy state
$\rho^{\prime}=\mathcal{N}(\rho)$. From the deconvolution formula and the
computation of the inverse adjoint PTM, it follows that

$\langle\sigma_{x}^{\otimes
2}\rangle_{\rho}=f(\eta,\mu)\big{\\{}\big{[}2\eta(1-\mu)+\mu(\sqrt{\eta}+1)\big{]}\langle\sigma_{x}^{\otimes
2}\rangle_{\rho^{\prime}}\\\ +\mu(\sqrt{\eta}-1)\langle\sigma_{y}^{\otimes
2}\rangle_{\rho^{\prime}}\big{\\}}\ ,$ (38)

while $\langle\sigma_{y}^{\otimes 2}\rangle_{\rho}$ is obtained from Eq. 38
through the index exchange $x\leftrightarrow y$, and that

$\langle\sigma_{z}^{\otimes
2}\rangle_{\rho}=g(\eta,\mu)\big{[}(\mu-1)^{2}(\eta-1)^{2}+\langle\sigma_{z}^{\otimes
2}\rangle_{\rho^{\prime}}\\\
-(\mu-1)(\eta-1)\langle\mathds{1}\otimes\sigma_{z}+\sigma_{z}\otimes\mathds{1}\rangle_{\rho^{\prime}}\big{]}\
,$ (39)

with

$\displaystyle
f(\eta,\mu)=\frac{1}{2\left[\mu\left(\eta-\sqrt{\eta}\right)-\eta\right]\left[\mu(\eta-1)-\eta\right]}\
,$ (40) $\displaystyle
g(\eta,\mu)=\frac{1}{\left[\eta+\mu(1-\eta)\right]^{2}}\ .$ (41)

Note that when the qubits are completely uncorrelated ($\mu=0$) the noise
deconvolution factorizes in terms of the single-qubit contributions.222See [9]
for the deconvolution of the single-qubit amplitude damping channel, with
$\eta=1-\gamma$. When the two qubits are completely correlated ($\mu=1$) the
noise has no effect on $\langle\sigma_{z}^{\otimes 2}\rangle$.

##  V Conclusions

We illustrated a procedure for the deconvolution of multiqubit noise described
by mathematically invertible quantum channels: it returns the ideal
expectation value of arbitrary observables from the noisy data.

In the case of Pauli channels, our prescription bypasses the inversion of the
noise map, providing the deconvolution from a set of Pauli measurements
rescaled by a few components of the PTM (where the number of factors is
quadratically reduced with respect to the general scenario). As discussed and
shown in the simulation, this analysis can be applied to any example of
multiqubit Pauli noise, e.g. the bit-flip or the depolarizing correlated
channels. Then, we presented a characterization technique that provides the
necessary PTM entries as a set of Pauli measurements on a specific class of
input states, guaranteeing the deconvolution whenever a theoretical approach
is not possible, e.g. when the parameters of the channel are unknown.

Finally, we discussed the deconvolution of noise that does not belong to the
class of multiqubit Pauli channels. Our procedure successfully applies also to
these cases, while less efficiently and requiring the complete computation
(and inversion) of the adjoint PTM.

## Acknowledgments

This work received support from MIUR Dipartimenti di Eccellenza 2018-2022
Project No. F11I18000680001, from EU H2020 QuantERA ERA-NET Cofund in Quantum
Technologies, Quantum Information and Communication with High-dimensional
Encoding (QuICHE) under Grant Agreements No. 731473 and No. 101017733, from
the U.S. Department of Energy, Office of Science, National Quantum Information
Science Research Centers, Superconducting Quantum Materials and Systems Center
(SQMS) under Contract No. DE-AC02-07CH11359, and from the National Research
Centre for HPC, Big Data and Quantum Computing (ICSC). We thank S. Mangini for
helpful discussions.

*

##  APPENDIX A POSITIVITY OF THE CHARACTERIZATION STATE

In this appendix we show that the operator used in the noise characterization
procedure fulfills all the requirements of a density operator, i.e.
Hermiticity, unit trace and positive semidefiniteness [1], namely, that it
represents a state of the $n$-qubit system.

Consider the operators of Eq. 16, with $k=1,2,...,d^{2}-1$. While Hermiticity
is straightforward, unit trace follows from $\Tr[\mathcal{P}_{k}]=0$ for
$k\neq 0$.

To prove positive semidefiniteness, we refer to the characterization of
density operators in terms of the coherence vector representation discussed in
[30]. Let $S_{m}$ be the coefficients of the characteristic polynomial of the
$d\times d$ matrix representation of $\rho$,333In this appendix we use the
standard operator framework rather than the vectorized representation. given
by

$S_{m}=\frac{1}{m}\sum_{j=1}^{m}(-1)^{j-1}\Tr[\rho^{j}]S_{m-j}\ ,$ (42)

with $S_{0}=1$ and $m=0,1,...,d$.444Here Eq. 42 is obtained by writing Eq.
(24) from [30] as a series, with $S_{0}=1$ the coefficient of the highest-
order term in the characteristic polynomial of $A$.

A necessary and sufficient condition for $\rho$ being positive semidefinite is
that $S_{m}\geq 0$ $\forall m$.

######  Theorem 1.

The operator $\rho=(1+\mathcal{P}_{k})/d$ is positive semidefinite.

###### Proof.

Consider the operator $A=(1+\mathcal{P}_{k})/2$, which is positive
semidefinite if and only if $\rho$ does. A direct computation yields $A^{j}=A$
and $\Tr[A^{j}]=d/2$. For $A$, Eq. 42 reads

$S_{m}=\frac{d}{2m}\sum_{j=1}^{m}(-1)^{j-1}S_{m-j}\ .$ (43)

Extracting the first term of the series and collecting a minus sign, we get

$S_{m}=\frac{d}{2m}\left(S_{m-1}-\sum_{j=2}^{m}(-1)^{j-2}S_{m-j}\right)\ .$
(44)

Translating $j\to j+1$ (with the sum now running from $1$ to $m-1$) and
writing Eq. 43 for $S_{m-1}$, we obtain a recursive relation

$S_{m}=\frac{d}{2m}\left(1-\frac{2(m-1)}{d}\right)S_{m-1}\ .$ (45)

With $\delta=1+d/2$, this yields

$\displaystyle S_{m}>0\quad\text{for}\ 0\leq m<\delta\ ,$ (46) $\displaystyle
S_{m}=0\quad\text{for}\ \delta\leq m\leq d\ ,$ (47)

which implies that $S_{m}\geq 0$ $\forall m$. This result does not depend on
$k$, nor do the coefficients $S_{m}$ and the roots of the characteristic
polynomial, i.e. the eigenvalue of $A$. Hence, $A$ is positive semidefinite
for any choice of $\mathcal{P}_{k}$ and consequently so is $\rho$. ∎

## References

  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_, 10th ed. (Cambridge University Press, 2010).
  * Ariano  _et al._ [2003] G. M. D. Ariano, L. Maccone, and M. Paini, Spin tomography, J. Opt. B 5, 77 (2003).
  * Mohseni  _et al._ [2008] M. Mohseni, A. T. Rezakhani, and D. A. Lidar, Quantum-process tomography: Resource analysis of different strategies, Phys. Rev. A 77, 032322 (2008).
  * Bongioanni  _et al._ [2010] I. Bongioanni, L. Sansoni, F. Sciarrino, G. Vallone, and P. Mataloni, Experimental quantum process tomography of non-trace-preserving maps, Phys. Rev. A 82, 042307 (2010).
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Rev. Mod. Phys. 86, 153 (2014).
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Error mitigation for short-depth quantum circuits, Phys. Rev. Lett. 119, 180509 (2017).
  * Endo  _et al._ [2018] S. Endo, S. C. Benjamin, and Y. Li, Practical quantum error mitigation for near-future applications, Phys. Rev. X 8, 031027 (2018).
  * Endo  _et al._ [2021] S. Endo, Z. Cai, S. C. Benjamin, and X. Yuan, Hybrid quantum-classical algorithms and quantum error mitigation, J. Phys. Soc. Jpn. 90, 032001 (2021).
  * Mangini  _et al._ [2022] S. Mangini, L. Maccone, and C. Macchiavello, Qubit noise deconvolution, EPJ Quantum Technol. 9, 29 (2022).
  * Macchiavello and Palma [2002] C. Macchiavello and G. M. Palma, Entanglement-enhanced information transmission over a quantum channel with correlated noise, Phys. Rev. A 65, 050301(R) (2002).
  * Macchiavello  _et al._ [2004] C. Macchiavello, G. M. Palma, and S. Virmani, Transition behavior in the channel capacity of two-quibit channels with memory, Phys. Rev. A 69, 010303(R) (2004).
  * Macchiavello and Sacchi [2016] C. Macchiavello and M. F. Sacchi, Witnessing quantum capacities of correlated channels, Phys. Rev. A 94, 052333 (2016).
  * Daems [2007] D. Daems, Entanglement-enhanced transmission of classical information in Pauli channels with memory: Exact solution, Phys. Rev. A 76, 012310 (2007).
  * Huang  _et al._ [2011] P. Huang, G. He, Y. Lu, and G. Zeng, Quantum capacity of Pauli channels with memory, Phys. Scr. 83, 015005 (2011).
  * Jiang  _et al._ [2021] J. Jiang, K. Wang, and X. Wang, Physical implementability of linear maps and its application in error mitigation, Quantum 5, 600 (2021).
  * Lautenbacher  _et al._ [2022] L. Lautenbacher, F. de Melo, and N. K. Bernardes, Approximating invertible maps by recovery channels: Optimality and an application to non-Markovian dynamics, Phys. Rev. A 105, 042421 (2022).
  * Gilyén  _et al._ [2022] A. Gilyén, S. Lloyd, I. Marvian, Y. Quek, and M. M. Wilde, Quantum algorithm for Petz recovery channels and pretty good measurements, Phys. Rev. Lett. 128, 220502 (2022).
  * Flammia and Wallman [2020] S. T. Flammia and J. J. Wallman, Efficient estimation of Pauli channels, ACM Trans. Quantum Comput. 1, 1 (2020).
  * Greenbaum [2015] D. Greenbaum, Introduction to quantum gate set tomography (2015), arXiv:1509.02921 [quant-ph] .
  * Nielsen  _et al._ [2021] E. Nielsen, J. K. Gamble, K. Rudinger, T. Scholten, K. Young, and R. Blume-Kohout, Gate set tomography, Quantum 5, 557 (2021).
  * Wood  _et al._ [2015] C. J. Wood, J. D. Biamonte, and D. G. Cory, Tensor networks and graphical calculus for open quantum systems, Quantum Inf. Comput. 15, 759 (2015).
  * Bény [2017] C. Bény, Quantum deconvolution, Quantum Inf. Process. 17, 26 (2017).
  * Audenaert and Scheel [2008] K. M. R. Audenaert and S. Scheel, On random unitary channels, New J. Phys. 10, 023011 (2008).
  * Siudzińska [2020] K. Siudzińska, Classical capacity of generalized Pauli channels, J. Phys. A Math. Theor. 53, 445301 (2020).
  * Cai  _et al._ [2020] Z. Cai, X. Xu, and S. C. Benjamin, Mitigating coherent noise using Pauli conjugation, NPJ Quantum. Inf. 6, 17 (2020).
  * Roncallo  _et al._ [2022] S. Roncallo, L. Maccone, and C. Macchiavello, Pauli transfer matrix direct reconstruction (2022), arXiv:2212.11968 [quant-ph] .
  * D’Arrigo  _et al._ [2013] A. D’Arrigo, G. Benenti, G. Falci, and C. Macchiavello, Classical and quantum capacities of a fully correlated amplitude damping channel, Phys. Rev. A 88, 042337 (2013).
  * D’Arrigo  _et al._ [2015] A. D’Arrigo, G. Benenti, G. Falci, and C. Macchiavello, Information transmission over an amplitude damping channel with an arbitrary degree of memory, Phys. Rev. A 92, 062342 (2015).
  * Hamada [2002] M. Hamada, A lower bound on the quantum capacity of channels with correlated errors, J. Math. Phys. 43, 4382 (2002).
  * Byrd and Khaneja [2003] M. S. Byrd and N. Khaneja, Characterization of the positivity of the density matrix in terms of the coherence vector representation, Phys. Rev. A 68, 062322 (2003).


# An analytical approach to symmetry breaking in multipole RF-traps

M. Marchenay and J. Pedregosa-Gutierrez and M. Knoop and M. Houssin and C.
Champenois  Aix Marseille Univ, CNRS, PIIM, Marseille, France
<EMAIL_ADDRESS>
###### Abstract

Radio-frequency linear multipole traps have been shown to be very sensitive to
mis-positioning of their electrodes, which results in a symmetry breaking and
leads to extra local minima in the trapping potential [1] disturbing the
operation of the trap. In this work, we analytically describe the RF-potential
of a realistic octupole trap by including lower order terms to the well-
established equation for a perfectly symmetric octupole trap. We describe the
geometry by a combination of identified defects, characterised by simple
analytical expressions. A complete equation is proposed for a trap with any
electrode deviation relying on a combination of the simple cases where the
defects are taken individually. Our approach is validated by comparison
between analytical and numerical results for defect sizes up to 4% of the trap
radius. As described in [2], an independent fine-tuning of the amplitude of
the RF voltage applied on each electrode can be used to mitigate the
geometrical defects of a realistic trap. In a different way than in [2], the
knowledge of an analytical equation for the potential allows to design the set
of RF-voltages required for this compensation, based on the experimental
measurement of the ion position in the trap, without information concerning
the exact position of each electrode, and with a small number of iterations.
The requirements, performances and limitations of this protocol are discussed
via comparison of numerical simulations and analytical results.

  * 10 december 2020

##  1 Introduction

The trapping of charged particles by radio-frequency (RF) electric fields has
proven to be both a powerful and versatile tool for experimental exploration
in physics and chemistry. Even if the linear quadrupole trap is mostly used
for mass spectroscopy [3], optical metrology [4], quantum computing [5] and
high precision measurements [6, 7] for its harmonic potential shape, linear
multipole traps happen to be candidate of choice in the design of micro-wave
atomic clocks for spacecraft navigation [8, 9] or for cold collisionand
chemical reaction studies [10, 11, 12, 13]. Both these applications take
advantage of the lower RF driven motion amplitude which results in a lower
Doppler effect and a lower kinetic energy in the collision frame, compared to
the same sample in a quadrupole trap. When laser cooled in a multipole trap,
the ions are expected to organise in a hollow core structure [14], forming
structures like rings [15] and tubes [16] that cannot be observed in
conventional quadrupole traps. Recent reported experimental observations in a
linear octupole trap reveal no such structure but three individual ion clouds.
This is the consequence of built in imperfection and mechanical electrodes
misalignment that result in a symmetry breaking in the RF-potential
arrangement of the trap [1]. A similar phenomenon has been observed in a
22-pole with the apparition of 10 zones where the ions cooled to 170 K
accumulate [17]. We aim at correcting the mechanical defects in the octupole
trap by adjusting each electrode potential as demonstrated in [2] in a
simulated experimental set-up. In this former work, a random geometry with a
mean displacement of each electrodes of 1.3% of the internal radius of the
trap could be corrected on a ring large as a tenth of this internal radius
after 576 iterations of trapping potential comparison and adjustment. The
simulations assume that laser-cooled ions are trapped and observed in the
three local minima of the pseudo-potential and the potential comparisons are
based on their relative positions. The compensation leads to a trapping
potential variation along this ring corresponding to 3.6 mK. To reduce
drastically the number of required steps to compensate for geometrical
defects, we assumed an accurate understanding of the correlation between trap
defects and experimental observations. To that purpose, we have worked on
establishing an analytic equation for the instantaneous RF potential of an
octupole trap with small defects. The resulting equation is entirely defined
by the real position of the electrodes in the radial plan, and the electrode
size. The validation of this approach is based on the comparison between the
analytical potential and the one computed numerically. This comparison focuses
on information that can be accessed by experiments, like the one that
demonstrated the three local minima in the trapping potential [1].

The method used as an analytical description for the rf potential is
introduced in the next section. Section 3 details how to convert the electrode
position into an equation, using a basis of five identified deformations that
can be superposed to describe any configuration. Section 4 explains the
correction protocol that we developed, based on the previous analytical
equations and a set of compensation voltages. The strategy first tested for
defect compensation can then be used to create on-demand configurations of
three parallel trapping zones.

##  2 A 2D-octupole potential as a combination of main and perturbation terms

The objective of this first part is to build an analytical equation for the
time dependent electric potential in the radial plane of an asymmetric linear
octupole trap by adding perturbations to the potential equation of an ideal
octupole trap. Any lack of parallelism between the electrodes is not taken
into account and we only consider a mismatch of electrode positions in a
transverse plane. Using a perturbative approach is motivated by the need to
have a description of the trapping pseudo-potential whereas the defects are
accurately described in the time dependent electric field. The starting point
is the equation for time-dependent potentials in a perfectly symmetric
multipole trap, written as in Eq. 1 where $V_{RF}$ is the amplitude of the RF
potential applied to the electrodes, $\Omega_{RF}/2\pi$ is the radio-frequency
and $K_{n}(x,y)$ is a 2D surface defined by the order $n$ of the multipole,
$2n$ being the number of electrodes in the trap.

$\psi_{n}(x,y,t)=V_{RF}K_{n}(x,y)\cos(\Omega_{RF}t)$ (1)

For an octupole ($n=4$), with an inner radius $r_{0}$ and the orthogonal
$(x,y)$ frame crossing the center of 4 electrodes out of 8, $K_{n}(x,y)$ is

$U_{4}(x,y)=\left(x^{4}-6x^{2}y^{2}+y^{4}\right)/r_{0}^{4}$ (2)

This equation assumes that the electrode shape fits the iso-potential lines,
and imposes a hyperbolic section for the electrodes [18]. However for
practical reasons, the use of circular electrodes is favoured with little
modification of the potential in the center of the trap if the ratio between
the electrode radius $r_{d}$ and the inner radius of the trap $r_{0}$ matches
a specific value that depends on the multipole trap order (0.333 for octupole
traps and 1.1451 for quadrupole traps)[19]. Cylindrical electrodes only are
within the scope of this article and the impact of the radius ratio is
simplified through this section by the use of coefficients $h_{0,p,l,c,h}$
assumed constant for a given radius ratio. The calculation of these
coefficients for different radius ratio is addressed in C.

With the analytical or numerical description of the RF-potential
$\psi_{n}(x,y,t)$, we compute the pseudo-potential using the method first
demonstrated for a quadrupole trap in [20] and extrapolated to a multipole
trap [21, 14]. In this approach, the particles are trapped in a static
potential which corresponds to the time-average of the squared local RF
electric field. One can show that the static pseudo-potential can be computed
by

$V^{*}(\mathbf{r})=\frac{q^{2}\mathbf{E_{0}}^{2}(x,y)}{4m\Omega_{RF}^{2}}.$
(3)

with $\mathbf{E_{0}}(x,y)$ the local amplitude of the RF-electric field,
defined as $-\mathbf{\nabla}(V_{RF}K_{n}(x,y))$. The experimental signature of
a symmetry breaking in a $2n$-pole trap, is the accumulation of trapped ions
into $(n-1)$ local minima, as was observed for an octupole trap [1] and a
22-pole trap [17]. This signature in the pseudo-potential approximation means
that the RF-electric field $E_{0}(x,y)$ has $(n-1)$ non-degenerated $(x,y)$
zeroes. This can happen only if the $n$-order $(x,y)$ polynomial RF-potential
$\psi_{n}$ includes terms of lower orders than $n$. These lower order terms
are cancelled when the $n$ mirror symmetries are obeyed.

Based on the previous analysis, we build the expression of the perturbative
term as a sum of lower order terms. Insight on this decomposition method can
be gained by starting with the simpler case of a symmetric quadrupole
potential with an inner radius $r_{0}$. We call $U_{2}$ the $K_{2}$ polynomial
if the axis frame is chosen crossing the center of the electrodes, and $V_{2}$
if it is equidistant to their center. As shown by Eq. 4 and 5 they can be
expressed as product and combination of the two dipole potential basis terms
$U_{1}(x,y)=x/r_{0}$ and $V_{1}(x,y)=y/r_{0}$.

$\displaystyle U_{2}(x,y)$ $\displaystyle=$
$\displaystyle\left(x^{2}-y^{2}\right)/r^{2}_{0}=U_{1}^{2}(x,y)-V_{1}^{2}(x,y)$
(4) $\displaystyle V_{2}(x,y)$ $\displaystyle=$
$\displaystyle\left(2xy\right)/r^{2}_{0}=2U_{1}(x,y)V_{1}(x,y)$ (5)

In the same manner, and within the same convention as for $(U_{2},V_{2})$ for
the relation of the frame orientation, the equation for the octupolar
potential can be seen as a combination of quadrupole terms and thus, dipole
terms (Eq. 6).

$\displaystyle U_{4}(x,y)$ $\displaystyle=$ $\displaystyle
U_{2}^{2}(x,y)-V_{2}^{2}(x,y)$ (6) $\displaystyle V_{4}(x,y)$ $\displaystyle=$
$\displaystyle 2U_{2}(x,y)V_{2}(x,y)$ (7)

We assume that the lower order terms responsible for extra zeroes in the
octupole potential equation can be mapped onto quadrupole and dipole
contributions, as it is the case for the symmetric geometry. We therefore
propose as an expression for an asymmetric octupole potential $U_{R}(x,y)$,
with a frame axis set crossing the electrode center when it is ”perfect”

$U_{R}(x,y)=h_{0}[U_{4}(x,y)-W(x,y)]$ (8)

with the perturbation potential

$W(x,y)=a_{1}U_{2}(x,y)+a_{2}V_{2}(x,y)+a_{3}U_{1}(x,y)+a_{4}V_{1}(x,y).$ (9)

The global $h_{0}$ scaling parameter accounts for the impact of radius ratio
$r_{d}/r_{0}$. Symmetry breaking could also be responsible for higher order
terms that modify the profile of the potential for large enough distance to
the trap center. We neglect these contributions for this work.

A single perturbation has a clear signature on the pseudo-potential, like
shown on Fig 1. In the first example, a purely quadrupole perturbation, coded
by ($a_{1}=0.1,\ a_{2}=a_{3}=a_{4}=0$) creates three equally spaced potential
minima aligned along the $x$-axis (for $a_{1}<0$ the minima are aligned along
the $y$-axis). Using the equations, one can show that the distance between the
minima is $r_{0}\sqrt{|a_{1}/2|}$ and this perturbation was already proposed
in [22] to trap three parallel ion strings and control their separation. The
second example adds a dipole term perturbation coded by ($a_{3}=-0.1,\
a_{1}=a_{2}=a_{4}=0)$, which creates a near-equilateral triangle organisation
of the minima in the potential. Even if the zeroes are not easily accessible
to an analytic solution, one can show that one of the minima sits on the
$x$-axis and the two other ones are equidistant from this $x$-axis. These
individual cases allow us to associate a linear configuration of potential
minima to a quadrupole perturbation and a triangle configuration to a dipole
perturbation.

Figure 1:  Pseudo-potential shape for two simple cases of added perturbation
(a) Quadrupole term with $a_{1}=0.1,\ a_{2}=a_{3}=a_{4}=0$, and (b) Dipole
term with $a_{3}=-0.1,\ a_{1}=a_{2}=a_{4}=0$, for a $r_{0}=4$ mm trap

##  3 From electrode position to perturbation terms for the potential

The next step of our demonstration is to quantify the weights $a_{i}$ on each
quadrupole and dipole term in accordance with the positions of the trap
electrodes. We aim to define a parameter set independent of an arbitrary
definition frame to establish an unambiguous characterisation of the trap
deformations. In order to build a self-consistent description of the trap
geometry, the electrode configuration is described by the composition of five
deformation classes, based on the superposition of two quadrupole-like
electrode sets. We call $S$ the quadrupole set aligned with the $(x;y)$ frame
and $T$ the quadrupole set turned by a $\pi/4$ angle from this frame. The five
deformations that we review in the following are Compression, Sliding and
Splitting which keep the orthogonality of the sets and Rotation and Shearing
which are angular deformations. The formers are straightforward to turn into
equations, but angular deformations make the problem more complex by coupling
the axis used for the analytical expression.

For each deformation class, we review its physical origin, its relevant
parameters and give as a first approach a simplified equation of the potential
when the trap is deformed by this defect alone. The full equation accounting
for several coupled deformations is detailed in the next part. The validation
of the proposed defect description is based on experimentally accessible data.
In an experiment, the trapping potential can be observed by the fluorescence
emitted by the ions located in the potential minima. By analogy, with our
numerical experiment, we compare the pseudo-potential created by the
analytical RF potential to the one computed by CPO [23], a code that solves
the Laplace equation based on a boundary condition imposed by the surface of
the electrodes. The dimensions of the trial trap are $r_{0}=4$ mm and
$r_{d}=1.5$ mm and this geometry is discretised with pixel size equal to
$r_{0}/1000$ (4 $\mu$m). To be relevant, this comparison must not rely on an
artificial frame center. To avoid this issue, each arrangement of the
potential minima is referenced to the center of the triangle or line they
form.

###  3.1 Compression

A compression describes a situation where there is a modification, in at least
one of the quadrupole sets, of the distance between two electrodes facing each
other, without displacement of the center of each line joining them (see Fig.
10 for a pattern). The CPO calculations show that the minima formed by such a
defect arrange on a line and remain centred on the center of the trap, exactly
like when a quadrupole contribution is added to the instantaneous potential
(Eq. 4 and Fig. 1). To write the perturbation $W_{c}$ induced on the
instantaneous potential by such a defect, we define $L_{S}$ (resp $L_{T}$),
the length difference between the two directions within the $S$-quadrupole set
(resp $T$-set), normalised by the averaged inner radius $\overline{r}_{0}$,
such that

$W_{c}=h_{c}\left[L_{S}U_{2}(x,y)+L_{T}V_{2}(x,y)\right].$ (10)

Here, $h_{c}$ is adjusted to match the minima positions between the analytic
expression and the CPO calculations and which only depends on the radius ratio
$r_{d}/r_{0}$, like discussed in C. For our trap dimensions $h_{c}=0.820$ and
the perturbation code is ($a_{1}=h_{c}L_{S}$, $a_{2}=h_{c}L_{T}$,
$a_{3}=a_{4}=0$). The sign assigned to the length modifications must be
consistent with the convention chosen for the $U_{2}$ and $V_{2}$
contributions. To that end, $L_{S}$ is positive if the length along the
$x$-axis is larger than the one along the $y$-axis and $L_{T}$ is positive if
the length along the $(x-y)$-axis is larger than the one along the
$(x+y)$-axis. The comparison between analytical and numerical calculations of
the pseudo-potential minima, shown on Fig. 2, confirms that this analytical
representation is relevant for a surface grid of 4 $\mu$m pixel size, for a
compression reaching 8% of the averaged inner radius $r_{0}$.

Figure 2: Distance to the trap centre of the two outside minima vs the
compression reduced parameter $L_{S}$ increased from 0 to 0.11 by step of
0.011. The black points are the results computed with the CPO software and the
red points are calculated from the analytical equation Eq. (10). The size of
the error bar corresponds to the pixel size for the representation, which is
3.64 $\mu$m. Trap inner radius $r_{0}=4$ mm and electrode radius $r_{d}=1.5$
mm

For the particular case of an equal compression along the two directions
defining the same quadrupole set, there is no extra minima to be observed.
Indeed, in that case $L_{S}$ and $L_{T}$ remain null but the averaged inner
radius is modified, changing the strength of the trapping potential without
changing the number of local minima. In the basic equations, the normalisation
factor is rescaled by changing $r_{0}$ to $\bar{r_{0}}$

###  3.2 Sliding

Like the compression, a sliding is a deformation described within one
quadrupole set. There, and contrary to compression, the distance between
facing electrodes is conserved but the lines joining them do not cross in
their center (noted $M^{S}$ in the $S$-set and $M^{T}$ in the $T$-set, see
also Fig. 10). The characteristic parameters $(x^{S}_{l},y^{S}_{l})$ and/or
$(x^{T}_{l},y^{T}_{l})$ are defined by the positions of these lines centre
relatively to $M^{S}(X^{S},Y^{S}$) and $M^{T}(X^{T},Y^{T}$), the crossing
points of these two lines in the $S$ and $T$-quadrupole sets. The definition
of $M^{S}$ and $M^{T}$ is necessary in a general case, but for a sliding alone
the two points are not separated. For a sliding in the $S$-set, $x^{S}_{l}$
(resp $y^{S}_{l}$) is the sum of the $x$ (resp $y$) coordinates relatively to
$M^{S}$ of the two lines centres of this set. The same definition is
transposed in the $T$-set. These coordinates are normalised by the distance
between the centres of the facing electrodes. This kind of deformation is
responsible for a triangular organisation of the potential minima, like
induced by a dipole term. With our previous convention, this perturbation is
then coded as $a_{1}=a_{2}=0$, $a_{3}=h_{l}(x^{S}_{l}+x^{T}_{l})$ and
$a_{4}=h_{l}(y^{S}_{l}+y^{T}_{l})$. The comparison between the pseudo-
potential computed by CPO and the one deduced from the instantaneous potential
equation

$W_{l}=h_{l}\left[\left(x^{S}_{l}+x^{T}_{l}\right)U_{1}(x,y)+\left(y^{S}_{l}+y^{T}_{l}\right)V_{1}(x,y)\right]$
(11)

is shown on Fig. 3 with the adjusted value $h_{l}=2.566$ and for the
particular case of a sliding up to 4% of $r_{0}$ in the vertical pair of
electrodes of the $S$-quadrupole set.

Figure 3: Positions of two out of the three minima along the line over which
they organise for small enough sliding of the electrodes set over the $y$-axis
in the $S$ quadrupole set. The reduced sliding parameter $y^{T}_{l}$ is
increased from 0 to $4\times 10^{-2}$ by step of $4\times 10^{-3}$. The black
points are the results computed with the CPO software and the red points are
calculated from the analytical equation (11). The error bars correspond to the
pixel size (4 $\mu$m).Trap inner radius $r_{0}=4$ mm and electrode radius
$r_{d}=1.5$ mm.

This comparison shows that a sliding up to 4% of $r_{0}$ can be completely
taken into account by an extra dipole term.

###  3.3 Splitting

The Splitting defect class corresponds to the separation of the centre of the
$S$ and $T$ quadrupole sets, without any deformation of each set. The relevant
parameters ($x_{0},y_{0}$) are the relative position of the $T$ set centre
relatively to the $S$ one, normalised to $r_{0}$. The pattern formed by the
three minima is similar to the one showed on Fig. 1.b and for small enough
splittings, each minimum settles along a direction that we use as a reference
to plot their positions on Fig. 4, for two out of three minima, because of the
symmetry of the figure. This geometric configuration is similar to the one
induced by a dipole perturbation and can be reproduced by adding a dipole term
$W_{p}^{1}$ to $U_{4}(x,y)$ such that

$W_{p}^{1}=h_{p}\left[x_{0}U_{1}(x,y)+y_{0}V_{1}(x,y)\right].$ (12)

With the adjusted coefficient $h_{p}=1.586$ and within the convention of Eq.
9, this perturbation is coded by $a_{1}=a_{2}=0$, $a_{3}=h_{p}x_{0}$ and
$a_{4}=h_{p}y_{0}$. The comparison between the pseudo-potential minima
resulting from this analytical description and from the CPO calculations are
shown on Fig. 4 for a splitting between the $T$ and $S$ center along the $y$
axis. It shows a slight mismatch, which increases with the size of the defect.
A better fit for larger defects can be reached by adding a correcting term of
the quadrupole kind and the total equation for a splitting perturbation that
we use in the following is

$W_{p}=W_{p}^{1}+\frac{h_{p}^{{}^{\prime}}}{\sqrt{x_{0}^{2}+y_{0}^{2}}}\left[(y_{0}^{2}-x_{0}^{2})U_{2}(x,y)+2x_{0}y_{0}V_{2}(x,y)\right].$
(13)

Figure 4:  Position of two out of the three minima along the line over which
they organise for small enough splitting between the $T$ and $S$ quadrupole
sets. The $y$ oriented splitting is measured by the reduced parameter $y_{0}$,
increased from 0 to $5.5\times 10^{-2}$ by step of $5.5\times 10^{-3}$. The
black crosses are the results computed with the CPO software and the red dots
are calculated from the analytical Eq. (12). The blue stars give the minimum
position calculated with the corrected equation Eq. (13). The error bars
correspond to the pixel size (4 $\mu$m) . Trap inner radius $r_{0}=4$ mm and
electrode radius $r_{d}=1.5$ mm.

The agreement between both calculations is not as good as for the other
defect, with a one pixel mismatch, but we keep this analytical description of
the splitting with a value of $h_{p}^{{}^{\prime}}=0.1$ for any tested
$r_{d}/r_{0}$ ratio.

###  3.4 Angular deformation

An angular deformation is more complicated to account for as the frame used to
describe the geometry is modified, which has an impact on the analytical
expression of other deformations. An angular deformation can always be seen as
a combination of the rotation of one quadrupole set compared to the other one
and a shearing effect within one or both quadrupole sets.

We start with the shearing effect alone as it does not modify the global
orientation of one quadrupole set relatively to the other one (see Fig. 10 for
a representation). One angle per set ($\beta_{S}$, $\beta_{T}$) is required to
characterise this defect pattern. In the case of a pure shearing, $\beta_{S}$
is the angle of the axis passing by the centres of the two electrodes of the
$S$-set closest to the $x$-axis of the frame to the $x$-direction, the other
trap axis being constrained to have a $-\beta_{S}$ angle to the $y$-direction
of the frame (the angle sign-convention being the trigonometric one). The
three aligned local minima induced by such a deformation is the signature of
an extra quadrupole term, and the minima configuration can be reproduced by
the analytical expression

$W_{h}=h_{h}\left[\beta_{T}U_{2}(x,y)+\beta_{S}V_{2}(x,y)\right]$ (14)

with the adjusted parameter $h_{h}=1.404$ fitting our trap dimensions. The
comparison between analytical and numerical results show an agreement for a
shearing angle up to its maximal size $\pi/15$ up to a the smaller pixel size
of $2.6~{}\mu$m. The coding of this perturbation is then
$a_{1}=h_{h}\beta_{T}$, $a_{2}=h_{h}\beta_{S}$, $a_{3}=a_{4}=0$

As for a rotation between the two quadrupole sets, defined by the smallest
angle between two neighbour electrodes within the octupole trap
$\pi/4-2|\delta|$, it does not induce any extra minimum in the pseudo-
potential. Nevertheless, we see in the next section that it has an impact when
combined with other geometric defects as it rotates the alignment of the
minima. In that case, $\delta$ must be oriented and it is positive for a
counter-clockwise rotation of the $S$-set. We now consider more realistic
situations that combine different basis defects.

###  3.5 A general form for any configuration

In the previous sections, an analytical expression is proposed for five
identified defect classes. We propose now to demonstrate that any geometric
configuration can be represented as a linear combination of the five
identified defects and that their expression can be summed up to analytically
describe this configuration. We first focus on situations where the defects do
not involve any angular deformation and where the three defects Compression,
Sliding and Splitting are combined. In that case, there is no cross-effect
between the three basic defects and the resulting impact on the potential is a
simple linear sum of the contribution of each individual defects. The
comparison between the numerical simulations by CPO and the analytical
expressions resulting from a simple sum of these three perturbation terms show
an agreement for defects resulting in a mis-position as large as 4% of the
inner radius. As an example, we propose to focus on the combination of the
dipole-kind sliding within the $S$-quadrupole set, defined by
$y_{l}^{S}=-0.004$ superposed with a quadrupole-kind compression within the
same set. The compression parameter $L_{S}$ is scanned from 0.108 to -0.111
with step of size -0.011. The positions of the three potential minima
calculated both by CPO and by the analytical description are shown on Fig. 5.
When the compression parameter is scanned, the organisation of the minima
starts from a triangle (labelled 1 minima) through a nearly balanced triangle
(labelled 11 minima) to a line oriented along the compression and sliding axis
(labelled 21 minima). The agreement between the analytical description and the
numerical calculations shows the relevance of our method for mis-positioning
up to 4% of the inner radius.

Figure 5: Positions in the radial plane of the three minima of the pseudo-
potential, computed by CPO (black cross) and by the analytical equation (red
dots). They are superposed on the same pixel (size 4 $\mu$m). The trap defects
are described by a constant sliding effect within the $S$-quadrupole set,
defined by $y_{S}^{Sd}=-0.004$, superposed to a compression within the same
set, defined by $L_{S}$ scanned from 0.108 (minima label 1) to -0.111 (minima
label 21) with step of size -0.011

If we now consider an angular deformation, the total perturbative term is not
a simple sum of the basic ones and the analytical expression of each
perturbative term must be modified to represent accurately the trap
configuration. For example, to account for the rotation $\delta$ of one
quadrupole set compared to the other one (like defined in 3.4), the quadrupole
terms, $U_{2}$ and $V_{2}$ have to be corrected in the expression of a
compression $W_{c}$, and shearing $W_{h}$ perturbation. We call
$\overline{R(\theta)}$ the matrix

$\left(\begin{array}[]{cc}\cos(\theta)&\sin(\theta)\\\
\sin(\theta)&\cos(\theta)\end{array}\right)$ (15)

which is not a conventional rotation matrix because it acts on quadrupole
terms. In case of a compression+rotation configuration, $W_{c}$ becomes
$W_{c}^{R}$ which can be analytically expressed by changing
$(U_{2}(x,y),V_{2}(x,y))$ into $(U_{2}^{c}(x,y),V_{2}^{c}(x,y))$ such that
$(U_{2}^{c};V_{2}^{c})=\overline{R(-2\delta)}(U_{2};V_{2})$. In case of a
shearing+rotation configuration, $W_{h}$ becomes $W_{h}^{R}$ with
$(U_{2}^{h};V_{2}^{h})=\overline{R(\delta)}(U_{2};V_{2})$. For a sliding
defect, the rotation effect is already included in the definition and the
splitting description does not need an axis rotation to remain accurate but we
observe a rescaling of the perturbative contributions of the three non-angular
defects. Altogether, the perturbation terms adding up to the octupole
instantaneous potential write

$\displaystyle W^{R}(x,y)$ $\displaystyle=$
$\displaystyle\left(1+\frac{|\delta|}{\pi}\right)W_{h}^{R}$ (16)
$\displaystyle+$
$\displaystyle\left(1+\frac{2|\delta|}{\pi}\right)\left[W_{c}^{R}(x,y)+W_{l}+W_{p}\right]$

The last step is to combine a shearing effect to the other deformations : in
the case of a shearing+compression defect, the three potential minima still
align around the trap center but the modification in the inner radius has an
impact on the scaling of the perturbation terms. More precisely, it is the
difference $\Delta r$ between the inner radius in the $T$-quadrupole set and
the $S$-quadrupole set that counts and $W_{h}$ must be coded like
$a_{1}=h_{h}\beta_{T}(1-3\Delta r/r_{0})$, $a_{2}=h_{h}\beta_{S}(1+3\Delta
r/r_{0})$, $a_{3}=a_{4}=0$.

In the case of a shearing+sliding defect, the minima organise on a triangle
and the analytical expression depends on the involved quadrupole sets. Indeed,
if these two defects are sported by the same set, the projection procedure
described in 3.2 accounts for the shearing effect and there is no need for an
extra calculation step. In the other case, the definition of the sliding
parameter depends on the shearing of the other set. As an example, if the
$T$-set is sheared by $\beta_{T}$ and the sets are deformed by a sliding
defined by ($x_{l}^{S},y_{l}^{S}$) and ($x_{l}^{T},y_{l}^{T}$), the
perturbation term $W_{l}$ must be coded as $a_{1}=a_{2}=0$,
$a_{3}=h_{l}((1+4\beta_{T}/\pi))x_{l}^{S}+x_{l}^{T})$ and
$a_{4}=h_{l}((1-4\beta_{T}/\pi))y_{l}^{S}+y_{l}^{T})$.

As for a combination of shearing+splitting defect, the minima also organise as
any triangle and the calculated potential matches the numerical calculation if
the perturbation $W_{p}$ is now coded by $a_{1}=a_{2}=0$,
$a_{3}=h_{p}(1+\sin(2\beta_{T})\cos(2\beta_{S})x_{0}-\sin(2\beta_{S})y_{0})$
and
$a_{4}=h_{p}(1-\sin(2\beta_{T})\cos(2\beta_{S})y_{0}-\sin(2\beta_{S})x_{0})$.
The equation summing all the corrected contributions is collected in B.

To illustrate a general case, we propose on Fig. 6 to observe this matching
for electrode positions that can be described by a compression
($L_{T}=0.55$)+splitting ($x_{0}=0.007$)+shearing with a shearing parameter
$\beta_{T}$ scanned from 0.088 to -0.088. Once again, the analytical
description and the numerical calculations agree within the 4 $\mu$m pixel
size.

Figure 6: Positions in the radial plane of the three minima of the pseudo-
potential, computed by CPO (black crosses) and by the analytical equation (red
dots). The electrode position defects are described by a fixed splitting
effect defined by $x_{0}=0.007$, superposed to a compression within the
$T$-quadrupole set, defined by $L_{T}=0.55$ and a shearing within the same set
scanned from $\beta_{T}=0.088$ (minima label 1) to $\beta_{T}=-0.088$ (minima
label 21) with step of size -0.009

###  3.6 Completeness of the defect basis

In the previous sections, agreement between numerical and analytical results
were demonstrated for designed geometric configurations. To confirm the
relevance of the description, we propose to check that any random electrode
configuration can be mapped onto the 5 basic defects. To that purpose, we
decomposed 200 electrode configurations onto the five defects basis. These
configurations were built such that the positions of the centers of the eight
electrodes are randomly chosen at a distance from their ideal positions
ranging 0.5% to 4% of the distance to the trap centre $r_{0}+r_{d}$. The
comparisons between numerical and analytical results are ordered by the
average distance $\overline{d}$ between the local potential minima calculated
with CPO and with the analytical description. The resolution is limited by the
4 $\mu$m pixel size relatively to the distance $r_{0}+r_{d}=5.5$ mm. To
compare very different potential minima configuration, we rescale the average
distance $\overline{d}$ relatively to the mean distance $d_{b}$ between the
three minima and their mass center, and write it as $\overline{d}_{s}$. For a
position mismatch of each electrode equal to 1% of $r_{0}+r_{d}$, 182
(respectively 198) cases out of 200 have a mean relative distance
$\overline{d}_{s}$ lower than 2% (respectively 4%). These ratios do not change
much when the mis-positioning is 2% of $r_{0}+r_{d}$. When it reaches 4%,
these numbers decrease to 132 cases out of 200 that have a mean distance
$\overline{d}_{s}$ lower than 2% and 173 cases with a mean distance lower than
4%. This very good statistics demonstrate that the description in terms of
defect basis is relevant to reproduce the potential minima configuration for
defect size of few % of the trap size. This limit is reachable to machined
macroscopic trap of millimetre size.

##  4 Control of the RF voltages

In this last section, we show how the analytical description that we propose
above allows the calculation of a set of relative RF-amplitude that must be
applied to compensate for the geometric defects. Then, starting from a
compensated potential, on can create designed line or triangle minima
configurations. The description of identified defects as dipole- or
quadrupole-kind shows that several geometrical configurations can be
responsible for the same potential minima configuration. As an example, a
compression described by $L_{s}=0.0554$, which corresponds to a 2% positioning
error relative to $r_{0}+r_{d}$, creates the same pattern as a shearing
described by $\beta_{T}=0.0319$, which corresponds to an angular mis-
positioning of 0.032 rad (with position resolution of $2\times 10^{-3}r_{0}$).
Our interest being in the shape of the potential itself, we developed a
strategy that bypass the identification of the geometric defects and the
positions of the electrodes as there is no one-to-one connection with a
minimum configuration.

###  4.1 Producing controlled patterns in a perfect octupole

The perturbation induced by the structure being a sum of quadrupole and dipole
terms, we deduce a set of applied voltages on the electrodes that generate
quadrupole and dipole perturbations in the potential. The compensation
strategy is based on the identification of the $a_{1,2,3,4}$ coefficients
involved in the analytical equation to apply the appropriate counter-
perturbation $-a_{1,2,3,4}$ through custom applied voltages on the electrodes.
The method we propose is to map a $a_{1,2,3,4}$ set to a set of voltage bias
$\delta V_{k}$ added to the average RF-voltage on each electrode, with $k$ its
index. Applying $-\delta V_{k}$ to the trap electrodes should, in principle,
compensate for the geometric defects. This is only partly verified because the
$-\delta V_{k}$ are applied on a trap showing a perturbed geometry but we
choose this strategy as a starting point.

There are different $\delta V_{k}$ patterns that can be applied to produce the
same perturbation in the potential, and one is given here that presents the
advantage of having the most accurate representation of dipole terms and
sharing the compensation between the largest number of electrodes :

$\displaystyle\delta
V_{0}=V_{rf}\left(+\frac{a_{1}}{q_{cal}}-\frac{a_{4}}{d_{cal}}\right);\ \delta
V_{1}$ $\displaystyle=$ $\displaystyle
V_{rf}\left(-\frac{a_{2}}{q_{cal}}-\frac{1}{\sqrt{2}}\frac{a_{3}+a_{4}}{d_{cal}}\right);$
$\displaystyle\delta
V_{2}=V_{rf}\left(-\frac{a_{1}}{q_{cal}}-\frac{a_{3}}{d_{cal}}\right);\ \delta
V_{3}$ $\displaystyle=$ $\displaystyle
V_{rf}\left(+\frac{a_{2}}{q_{cal}}-\frac{1}{\sqrt{2}}\frac{a_{3}-a_{4}}{d_{cal}}\right);$
(17) $\displaystyle\delta
V_{4}=V_{rf}\left(+\frac{a_{1}}{q_{cal}}+\frac{a_{4}}{d_{cal}}\right);\ \delta
V_{5}$ $\displaystyle=$ $\displaystyle
V_{rf}\left(-\frac{a_{2}}{q_{cal}}+\frac{1}{\sqrt{2}}\frac{a_{3}+a_{4}}{d_{cal}}\right);$
$\displaystyle\delta
V_{6}=V_{rf}\left(-\frac{a_{1}}{q_{cal}}+\frac{a_{3}}{d_{cal}}\right);\ \delta
V_{7}$ $\displaystyle=$ $\displaystyle
V_{rf}\left(+\frac{a_{2}}{q_{cal}}+\frac{1}{\sqrt{2}}\frac{a_{3}-a_{4}}{d_{cal}}\right);$

The electrodes are numbered clockwise, starting by the bottom. The linear
dependency of the $\delta V_{k}$ with the $a_{1,2,3,4}$ was established
through numerical simulations by generating potentials according to Eq. 9 with
varying $a_{1,2,3,4}$ coefficients and collecting the positions of the minima
as reference. The CPO software was then used to generate the potential map of
a perfect octupole trap ($(r_{0},~{}r_{d})$=(4, 1.5) mm) with individually
tuneable potentials on the electrodes. The potentials on the electrodes of the
simulated octupole trap were applied according to Eq. 17, with adjusted values
of the $d_{c}$ and $q_{c}$ calibration coefficient to ensure the superposition
of the positions of the potential minima between the surface equation and the
CPO simulations. With this protocol we find $d_{cal}=0.912\pm 0.008$ and
$q_{cal}=0.796\pm 0.008$, the uncertainty comes from the range of values that
allows for a superposition at the pixel size which is reduced to 1 $\mu$m to
improve the precision. The fit was ensured for values of $a_{i}$ up to $0.1$
which, for the sake of comparison, is a little larger than the amplitude
expected for a single involved deformation class reaching of 4% of the
$(r_{0}+r_{d})$ distance if the potential deformation were of mechanical
origin. This ought to ensure we are not out of range of the validity domain of
this fit when we apply the correction scheme in our trap where the defect
amplitude is expected to be about 2% of the $(r_{0}+r_{d})$ distance at most.
Fig. 7 shows a fit to the pixel size between minima generated from the surface
equation Eq. 9 for varying $a_{1,2,3,4}$ coefficients and the minima generated
in the perfect octupole trap by converting the $a_{1,2,3,4}$ coefficients into
adapted RF voltages according to Eq. 17. The selected tuning patterns allows
the controlled generation of any perturbation in the octupole trap that writes
as a sum of quadrupole and dipole terms.

Figure 7: Red crosses: positions of the minima in the pseudo-potential
generated from Eq. 13 for $a_{1,2,3,4}$ coefficients following the dependency
: $\\{0.4(0.1-j),~{}0.2j,~{}-0.1(0.1-j),~{}0.2j\\}$ with $j$ spanning from 0
to 0.1 with steps of 0.01. Black dots: corresponding positions of the minima
for the simulated perfect octupole trap with the calculated RF voltages
$V_{k}=V_{rf}+\delta V_{k}$ according to Eq. 17. The size of the grid in both
cases is $d_{px}=$1 $\mu$m.

###  4.2 Diagnosing the perturbations from the positions of the minima

In order to compensate the perturbations of the potential by applying the
proper counter-perturbation through a set of adapted voltages, it is necessary
to establish a diagnosis of the perturbations from experimentally accessible
data. Working directly with the perturbations of the potential instead of the
deformation of the structure presents the advantage of reducing the number of
unknown parameters from 16 to five ($a_{1,2,3,4}$ and $\bar{r_{0}}$). In the
context of an experiment, we can bypass the knowledge of the positions of the
electrodes to find the $a_{i}$ coefficients by working directly with the
relative positions of the three potential minima in the octupole frame,
assuming the trapped ions bunch in these minima once laser cooled [2]. The set
of parameter ($a_{1,2,3,4}$, $\bar{r_{0}}$) being unique to a given
perturbation of the potential, we use the position of the minima to evaluate
the parameter set $a_{1,2,3,4}$ by assuming a value for $\bar{r_{0}}$, usually
$\bar{r_{0}}=r_{0}$. Starting from guess initial values for the $\\{a_{i}\\}$,
the pseudo-potential minima positions generated from Eq. 9 are compared to the
observed minima. Like in [2], the figure of merit is the average of the three
distances $\bar{d}$ between the calculated and the measured minima and an
optimisation routine looks for the set of $a_{i}$ coefficients which minimises
this figure.

To validate the possibility to deduce the perturbations from the positions of
the minima we have generated a number of test cases with known perturbations
and applied the search protocol to see if the returned coefficients correspond
to the inputs values, with an expected result $\bar{d}=0$. Three different
pixel sizes were tested ($d_{px}=$ 2 $\mu$m, 4 $\mu$m and 8 $\mu$m) to
evaluate the impact of the resolution on the positioning of the minima,
determined in practice by the collection optics. The positions of the minima
were collected from calculated potentials with randomised $a_{i}$ coefficient
with realistic order of magnitude, which corresponds to $\pm 0.1$ for a 4%
amplitude on the defects and $\pm 0.05$ for a 2% case. In our tests,
coefficients can have any value in the interval $[-0.1,0.1]$. 50 randomised
cases were tested for each pixel sizes. In regard to the sampling of the
surface, a success condition for the search protocol can be set to
$\bar{d}<d_{px}$ or $\bar{d}<2\times d_{px}$. The success rate of the search
protocol is summarised in Table. 1.

| Success rate if : | Success rate if :
---|---|---
$d_{px}$ ($\mu$m) |  $\bar{d}<d_{px}$ (%) |  $\bar{d}<2\times d_{px}$ (%)
$2$ | 78 | 96
$4$ | 84 | 92
$8$ | 66 | 90
Table 1: Success rate of the code predicting the $\\{a_{i}\\}$ coefficients
pondering the perturbations from the positions of the minima in the resulting
pseudo-potential.

To estimate the closeness of the $\\{a_{i}\\}$ coefficients returned by the
code to the input coefficients, we use the standard deviation
$\sigma_{1,2,3,4}$ of their difference over the cases fulfilling the success
condition $\bar{d}<d_{px}$. The results are summarised in Table. 2, and show
that a smaller pixel size leads to a better accuracy on the estimated
coefficients. We validate our search protocol as adequate to the purpose of
identifying the perturbations in the potential from the position of the minima
in the frame of our compensation protocol.

$d_{px}$ ($\mu$m) | $\sigma_{1}$ | $\sigma_{2}$ | $\sigma_{3}$ | $\sigma_{4}$
---|---|---|---|---
$2$ | 0.000293 | 0.000370 | 0.000223 | 0.000249
$4$ | 0.000820 | 0.000878 | 0.000470 | 0.000420
$8$ | 0.001787 | 0.001785 | 0.000973 | 0.000757
Table 2: Standard deviation of the difference between the input coefficients
and the coefficients estimated from the search protocol, for the ‘successful
cases’ fulfilling the condition $\bar{d}<d_{px}$.

###  4.3 Compensation of geometric defects

To test the compensation in an octupole trap affected with structural
deformations, we impose that the generated geometries are free from a general
rotation. From the positions of the minima we can diagnose the weights
$a_{1,2,3,4}$ characterising the structure induced perturbation of the
potential and with the set of equations (17) we can apply the appropriate
counter-perturbation characterised by $-a_{1,2,3,4}$ by tuning the voltages on
the electrodes to compensate for the structure deformations. Simulations have
shown that compensation of the potential is not achieved after one step of
this protocol, which is due to the deformation of the structure modifying
slightly the value of the applied counter-perturbation in the potential. Fig.
8 shows a practical example in an octupole trap with a random displacements of
the electrodes with a norm of 2% of the $(r_{0}+r_{d})$ distance : the
positions labelled 0 correspond to the positions of the minima before
compensation is applied, and their average distance to their barycentre is
$\bar{d}_{b}=$ 830 $\mu$m. After a first correction step, the minima, labelled
1, are brought closer to the centre but not merged since they are still
distant of $\bar{d}_{b}=$ 250 $\mu$m (see Fig. [fig:dist_step_conv]).

Figure 8: Positions of the pseudo-potential minima in the radial plan for ten
step of the iterative correction protocol, one color per step, on a structure
with 2% random error on the position of the electrodes. The right panel is a
zoom on the central region of the left panel. The black dots correspond to the
position of the minima at the end of the tenth iteration of the protocol.
Figure 9: Average distance $\bar{d}$ of the minima to their barycentre
calculated from the position plotted in Fig. 8, for each step of the iterative
correction protocol.

By repeating the correction protocol several times, each iteration ought to
allow the potential to be brought closer to the ideal shape by getting finer
information on the tuning voltages. More precisely, at each step $j$ the
correction protocol goes as follow: the positions of the minima are gathered
from the fluorescence of the ions in the trap. The associated coefficients set
$\\{a_{i}\\}$ characterising the present potential deformations is determined
via the search protocol. A correction set $\\{a_{i}\\}_{j}$ is calculated and
used to update the correction voltages by building up with the previous values
and the total correction voltage on each electrode is cumulative with the
previous steps. To illustrate the step by step functioning of the correction
protocol, ten iterations of the correction protocol are applied on the test
trap. Fig. 8 shows the positions of the minima in the radial plane for each
correction step, and the average distance of the minima to their barycentre is
plotted on Fig. 9. The first three correction steps reduce the average
distance between the minima of 93.5% and the next steps further reduce it of
another 3%. Given the pixel size of 4 $\mu$m, the average distance between the
minima reaches its limit at the seventh iteration, with a value of 24 $\mu$m.
Iterating the protocol further rearranges the positions of the minima but does
not bring them closer to merging.

Two parameters govern the achievable quality of the compensation: they are the
pixel size, and the allowed resolution on the coefficients $a_{i}$. In an
experiment they are limited by the resolution of the optics and the ion motion
for the positioning of the minima, and the available resolution on the tuning
voltages $\delta V_{k}$ (see Eq. 17). The optimum scenario for the
compensation was identified for pixel sizes of 4 $\mu$m, and a resolution of
$10^{-6}$ on the $a_{i}$. In both cases there is a notable deterioration of
the results for bigger pixel sizes and lower resolutions, but no notable gain
is achieved with more stringent requirements. With this pixel size and
resolution, five iterations were sufficient, on average, to see no further
gain in the average distance $\bar{d}_{b}$ over 10 test cases with a 2%
perturbation range. The average distance between the minima keeps around 40
$\mu$m on average, which means that even in the best-case parametrisation the
spatial merging of the minima is not achieved. Nevertheless, regarding their
potential depth, they are reduced to a few microkelvins compared to their
initial depths of the order of a few hundreds of kelvins. The microkelvin
scale is negligible in regard of the temperature of laser cooled ions which is
of the order of 10 mK and we can consider that the three local minima in the
potential are compensated after five iterations.

So far, the contribution of the DC-deconfining term induced by the DC-trapping
along the trap axis [24] has been neglected in the protocol even if its impact
on the trapping potential shape can not be avoided. By turning a flat bottom
well into ring shaped geometry, it smears the three potential minima into two
or even one, setting a limit to the compensation protocol that stops before
the best estimated performance are reached [2]. The efficiency of the proposed
method is then dependent on the maximum RF amplitude and minimal DC-voltage
that the set-up can reach to observe ions in the most favorable configuration
for a relevant compensation. Once the tuning parameters identified, the
operating voltages can be modified as the compensation voltages can be adapted
accordingly. With a DC-voltage added to the RF-pseudo-potential, it was found
in [2] that a relative resolution of $10^{-4}$ on $\delta V_{k}$ was
sufficient to reach the best reachable performances, with residual potential
variation along the ring of the order of few millikelvins.

##  5 Conclusion

The description of a realistic linear octupole trap proposed in this paper is
based on the adjunction of lower order contributions in the RF potential
created by the electrodes. The protocol for diagnosing of the impact of
structural defects on the potential and their compensation relies on the
positioning of the minima in the RF-induced pseudo-potential through the
imaging of the trapped cold ions by their laser induced fluorescence. This
protocol was motivated by the practical realisation of the compensation of the
potential asymmetries introduced by structural defects, to form a trapping
potential with a shape as symmetric as possible. In practice, it is limited by
the impact of the DC-deconfining contribution in the radial plane that
modifies the ion organisation when they are few. In the large number limit,
they organise independently of this contribution [14] and another strategy
needs to be found to go further in the defect compensation. Nevertheless, the
proposed description is a tool that allows for the control design of three ion
chains [22], a relevant configuration to simulate and study 3D frustrated spin
system like demonstrated in [25] for 2D system where the transition from
paramagnetic to ferromagnetic order was observed for 3 ions.

## acknowledgement

MM acknowledges DGA for funding and stimulating scientific environnement. This
experiment has been financially supported by CNES (contract 151084).

##  Appendix A Geometrical representation of the defect basis set

Figure 10: Representation of the parameters defining a defect for a) :
Compression, b) : Sliding, c): Splitting, d): Shearing, and e)Rotation.

##  Appendix B Complete equation for a configuration involving the five
identified defects

This is the full equation like described in Eq. 9 :

$\displaystyle a_{1}$ $\displaystyle=$ $\displaystyle
h_{c}\left(1+\frac{2|\delta|}{\pi}\right)\left(L_{S}\cos(2\delta)-L_{T}\sin(2\delta)\right)$
$\displaystyle+$ $\displaystyle
h_{h}\left(1+\frac{|\delta|}{\pi}\right)\left(\beta_{T}\left(1-3\frac{\Delta
r}{r_{0}}\right)\cos(\delta)+\beta_{S}\left(1+3\frac{\Delta
r}{r_{0}}\right)\sin(\delta)\right)$ $\displaystyle-$ $\displaystyle
0.1\frac{x_{0}^{2}-y_{0}^{2}}{\sqrt{x_{0}^{2}+y_{0}^{2}}}$ $\displaystyle
a_{2}$ $\displaystyle=$ $\displaystyle
h_{c}\left(1+\frac{2|\delta|}{\pi}\right)\left(L_{T}\cos(2\delta)-L_{S}\sin(2\delta)\right)$
$\displaystyle+$ $\displaystyle
h_{h}\left(1+\frac{|\delta|}{\pi}\right)\left(\beta_{T}\left(1-3\frac{\Delta
r}{r_{0}}\right)\sin(\delta)+\beta_{S}\left(1+3\frac{\Delta
r}{r_{0}}\right)\cos(\delta)\right)$ $\displaystyle+$ $\displaystyle
0.1\frac{2x_{0}y_{0}}{\sqrt{x_{0}^{2}+y_{0}^{2}}}$ $\displaystyle a_{3}$
$\displaystyle=$
$\displaystyle[h_{l}\left(\left(1+4\frac{\beta_{T}}{\pi}\right)x^{S}_{l}+x^{T}_{l}-4\frac{\beta_{T}}{\pi}y^{T}_{l}\right)$
(20) $\displaystyle+$ $\displaystyle
h_{p}(1+\sin(2\beta_{T})\cos(2\beta_{S}))x_{0}-\sin(2\beta_{S})y_{0}]\left(1+\frac{2|\delta|}{\pi}\right)$
$\displaystyle a_{4}$ $\displaystyle=$
$\displaystyle[h_{l}\left(\left(1-4\frac{\beta_{T}}{\pi}\right)y^{S}_{l}+y^{T}_{l}-4\frac{\beta_{T}}{\pi}x^{T}_{l}\right)$
(21) $\displaystyle+$ $\displaystyle
h_{p}(1-\sin(2\beta_{T})\cos(2\beta_{S}))y_{0}-\sin(2\beta_{S})x_{0}]\left(1+\frac{2|\delta|}{\pi}\right)$

##  Appendix C Scaling the perturbation terms to any radius sizes

The comparisons between analytical expressions and numerical calculations
shown in the previous section where done for $r_{d}=1.5$ mm and $r_{0}=4$ mm
and we stated that the radius scale effects were all taken into account
through the scaling constants $h_{0,c,l,p,h}$. The correlations between
several numerical calculations run with 1 mm$<r_{0}<$5 mm for $r_{d}=0.5$ mm
and 0.1 mm$<r_{d}<$1.9 mm for $r_{0}=4$ mm show that these constants only
depend on the ratio $r_{d}/r_{0}$. For configurations with only one basic
defect introduced in the electrode geometry, the dependence of constants
$h_{0,c,l,p,h}$ with this ratio is computed by scanning $r_{d}$ from
$0.02r_{0}$ to $0.55r_{0}$ for $r_{0}=4$ mm. An exemple of such a numerical
study is shown on Fig. 11 for two splitting defects. The agreement between the
$h_{p}$ coefficient for two amplitudes of the defect confirms that they scale
for $r_{d}/r_{0}$ only.

Figure 11: Coefficient $h_{p}$ scaling the splitting perturbative term vs the
trap radius ratio $r_{d}/r_{0}$. The blue and red points correspond to a
splitting defect of respectively $y_{0}=0.055$ and $y_{0}=0.0275$ which
corresponds in the case of an octupole with $(r_{d},r_{0})=$(1.5 mm, 4 mm) to
a distance between the S and T-sets of $110~{}\mu$m and $220~{}\mu$m (2% and
4% of the $r_{d}+r_{0}$ total length). The $r_{d}/r_{0}$ error bar is due to
the $4~{}\mu$m pixel size of the simulated CPO surface. The $h_{p}$ error bar
covers the different values that allow us to reproduce the position of the
minima within one pixel. The continuous lines are the best fitted polynomial
curves.

Within the $4~{}\mu$m pixel size, the dependence of each scaling constant with
$r_{d}/r_{0}$ can be fitted by a polynomial equation like

$\displaystyle h_{0}$ $\displaystyle=$
$\displaystyle-0.4557-7.028(r_{d}/r_{0})+53.7(r_{d}/r_{0})^{2}-254(r_{d}/r_{0})^{3}\dots$
$\displaystyle+687(r_{d}/r_{0})^{4}-973(r_{d}/r_{0})^{5}+556(r_{d}/r_{0})^{6}$
$\displaystyle h_{c}$ $\displaystyle=$
$\displaystyle+0.565+1.138(r_{d}/r_{0})-2.073(r_{d}/r_{0})^{2}+3.021(r_{d}/r_{0})^{3}-1.995(r_{d}/r_{0})^{4}$
$\displaystyle h_{l}$ $\displaystyle=$
$\displaystyle+1.141+4.869(r_{d}/r_{0})-5.880(r_{d}/r_{0})^{2}+10.696(r_{d}/r_{0})^{3}-6.975(r_{d}/r_{0})^{4}$
$\displaystyle h_{p}$ $\displaystyle=$
$\displaystyle+0.614+5.639(r_{d}/r_{0})-16.260(r_{d}/r_{0})^{2}+29.980(r_{d}/r_{0})^{3}-22.084(r_{d}/r_{0})^{4}$
$\displaystyle h_{h}$ $\displaystyle=$
$\displaystyle+4.208+4.989(r_{d}/r_{0})-3.753(r_{d}/r_{0})^{2}-0.0025(r_{d}/r_{0})^{3}+2.1279(r_{d}/r_{0})^{4}$

## References

  * [1] J. Pedregosa-Gutierrez, C. Champenois, M. R. Kamsap, G. Hagel, M. Houssin, and M. Knoop.  Symmetry breaking in linear multipole traps.  Journal of Modern Optics, 0(0):1–9, 2017.
  * [2] J. Pedregosa-Gutierrez, C. Champenois, M. Houssin, M. R. Kamsap, and M. Knoop.  Correcting symmetry imperfections in linear multipole traps.  Review of Scientific Instruments, 89(12):123101, 2018.
  * [3] Donald J. Douglas, Aaron J. Frank, and Dunmin Mao.  Linear ion traps in mass spectrometry.  Mass Spectrometry Reviews, 24(1):1–29, 2005.
  * [4] S. M. Brewer, J.-S. Chen, A. M. Hankin, E. R. Clements, C. W. Chou, D. J. Wineland, D. B. Hume, and D. R. Leibrandt.  ${}^{27}{\mathrm{al}}^{+}$ quantum-logic clock with a systematic uncertainty below ${10}^{-18}$.  Phys. Rev. Lett., 123:033201, Jul 2019.
  * [5] Thomas Monz, Philipp Schindler, Julio T. Barreiro, Michael Chwalla, Daniel Nigg, William A. Coish, Maximilian Harlander, Wolfgang Hänsel, Markus Hennrich, and Rainer Blatt.  14-qubit entanglement: Creation and coherence.  Phys. Rev. Lett., 106:130506, Mar 2011.
  * [6] S. Alighanbari, G. Giri, F. Constantin, Vladimir Korobov, and Stephan Schiller.  Precise test of quantum electrodynamics and determination of fundamental constants with hd+ ions.  Nature, 581:1–7, 05 2020.
  * [7] Sayan Patra, M. Germann, J.-Ph. Karr, M. Haidar, L. Hilico, V. I. Korobov, F. M. J. Cozijn, K. S. E. Eikema, W. Ubachs, and J. C. J. Koelemeij.  Proton-electron mass ratio from laser spectroscopy of hd+ at the part-per-trillion level.  science, 369:1238, 2020.
  * [8] J.D. Prestage and G.L. Weaver.  Atomic clocks and oscillators for deep-space navigation and radio science.  Proceeding of the IEEE, 95(11):2235, 2007.
  * [9] E.A. Burt, W.A. Diener, and R.L. Tjoelker.  A compensated multi-pole linear ion trap mercury frequency standard for ultra-stable timekeeping.  IEEE Tr. on Ultrasonics, ferroelectrics ans frequency control, 55, 2008.
  * [10] Dieter Gerlich.  Ion-neutral collisions in a 22-pole trap at very low energies.  Physica Scripta, 1995:256, 01 2007.
  * [11] Roland Wester.  Radiofrequency multipole traps: tools for spectroscopy and dynamics of cold molecular ions.  J. Phys. B, 42(15):154001, 2009.
  * [12] Alan Günther, Pablo Nieto, David Müller, Alexander Sheldrick, Dieter Gerlich, and Otto Dopfer.  Berlintrap: A new cryogenic 22-pole ion trap spectrometer.  Journal of Molecular Spectroscopy, 332:8 – 15, 2017.  Molecular Spectroscopy in Traps.
  * [13] Malcolm Simpson, Markus Nötzold, Alice Schmidt-May, Tim Michaelsen, Björn Bastian, Jennifer Meyer, Robert Wild, Franco A. Gianturco, Milan Milovanović, Viatcheslav Kokoouline, and Roland Wester.  Threshold photodetachment spectroscopy of the astrochemical anion cn-.  The Journal of Chemical Physics, 153(18):184309, 2020.
  * [14] C. Champenois.  About the dynamics and thermodynamics of trapped ions.  J. Phys. B, 42:154002, 2009.
  * [15] C. Champenois, M. Marciante, J. Pedregosa-Gutierrez, M. Houssin, M. Knoop, and M. Kajita.  Ion ring in a linear multipole trap for optical frequency metrology.  Phys. Rev. A, 81(4):043410, 2010.
  * [16] F. Calvo, C. Champenois, and E. Yurtsever.  Crystallization of ion clouds in octupole traps: Structural transitions, core melting, and scaling laws.  Phys. Rev. A, 80(6):063401, 2009.
  * [17] R. Otto, P. Hlavenka, S. Trippel, J. Mikosch, K. Singer, M. Weidemueller, and R. Wester.  How can a 22-pole ion trap exhibit 10 local minima in the effective potential?  J. Phys. B, 42:154007, 2009.
  * [18] M H Friedman, A L Yergey, and J E Campana.  Fundamentals of ion motion in electric radio-frequency multipole fields.  Journal of Physics E: Scientific Instruments, 15(1):53–56, jan 1982\.
  * [19] A. J. Reuben, G. B. Smith, P. Moses, A. V. Vagov, M. D. Woods, D. B. Gordon, and R. W. Munn.  Ion trajectories in exactly determined quadrupole fields.  Int. J. Mass Spectrometry and Ion Processes, 154(1-2):43 – 59, 1996\.
  * [20] H.G. Dehmelt.  Radiofrequency spectroscopy of stored ions I: storage.  Advances in Atomic and Molecular Physics, 3:53–72, 1967.
  * [21] D. Gerlich.  Inhomogeneous rf fields: a versatile tool for the study of processes with slow ions.  In Cheuk-Yiu Ng and Michael Baer, editors, State-selected and state-to-state ion-molecule reaction dynamics, Part I, volume 82 of Advances in Chemical Physics Series. John Wiley and Sons, 1992.
  * [22] M. Marciante, C. Champenois, J. Pedregosa-Gutierrez, A. Calisti, and M. Knoop.  Parallel ion strings in linear multipole traps.  Phys. Rev. A, 83:021404, Feb 2011.
  * [23] F.H. Read and Inc. N.J. Bowring, Scientific Instrument Services.  Cpo-charged particle optics software.
  * [24] M. Drewsen and A. Brøner.  Harmonic linear Paul trap: Stability diagram and effective potentials.  Phys. Rev. A, 62(4):045401, 2000.
  * [25] A. Friedenauer, H. Schmitz, J. T. Glueckert, D. Porras, and T. Schaetz.  Simulating a quantum magnet with trapped ions.  nature physics, 4:757, 2008.