Vardan Martikyan

Camille Marin

Steffen J Glaser

Marc-André Delsuc

Dominique Sugny
email: dominique.sugny@u-bourgogne.fr

Application of Optimal Control Theory to Fourier Transform Ion Cyclotron Resonance

Keywords: Optimal control, Robust protocol, Ion Cyclotron Resonance

We study the application of Optimal Control Theory to Ion Cyclotron Resonance. We test the validity and the efficiency of this approach for the robust excitation of an ensemble of ions with a wide range of cyclotron frequencies. Optimal analytical solutions are derived in the case without any pulse constraint. A gradient-based numerical optimization algorithm is proposed to take into account limitation in the control intensity. The efficiency of optimal pulses is investigated as a function of control time, maximum amplitude and range of excited frequencies. A comparison with adiabatic and SWIFT pulses is done. On the basis of recent results in Nuclear Magnetic Resonance, this study highlights the potential usefulness of optimal control in Ion Cyclotron Resonance.

Introduction

Performing efficient and robust state control by means of external time-dependent system parameter is a fundamental challenge in many technological developments at macroscopic or microscopic scale [START_REF] Bryson | Applied optimal control[END_REF][START_REF] Bressan | Introduction to the Mathematical Theory of Control[END_REF][START_REF] Glaser | Training Schrödinger's cat: quantum optimal control[END_REF][START_REF] Schättler | Geometric optimal control: Theory, Methods and Examples[END_REF]. In this latter case, open-loop control protocol, i.e. without any real time feedback from the experiment during the control process, is generally used for practical and technical reasons. The controls are only designed from a modeling of the system dynamics and the efficiency of the control scenario may suffer from the accuracy of the theoretical description. The robustness of a control process with respect to experimental imperfections is therefore a key parameter in view of experimental implementation. Different techniques extending from adiabatic pulses to optimal control theory (OCT) have been developed in this open-loop framework to find the pulse parameters [START_REF] Bryson | Applied optimal control[END_REF][START_REF] Glaser | Training Schrödinger's cat: quantum optimal control[END_REF][START_REF] Brif | Control of quantum phenomena: past, present and future[END_REF][START_REF] Daems | Robust Quantum Control by a Single-Shot Shaped Pulse[END_REF][START_REF] Van Damme | Robust NOT-gate by single-shot shaped pulses: Demonstration by rephasing of atomic coherences[END_REF].

Optimal control tackles the question of bringing a dynamical system from one state to another with minimum expenditure of time and resources [START_REF] Bryson | Applied optimal control[END_REF][START_REF] Bressan | Introduction to the Mathematical Theory of Control[END_REF][START_REF] Glaser | Training Schrödinger's cat: quantum optimal control[END_REF][START_REF] Schättler | Geometric optimal control: Theory, Methods and Examples[END_REF]. The modern version of OCT was born in the sixties with the Pontryagin Maximum Principle (PMP), which provides a general and rigorous mathematical framework for optimal control techniques [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF][START_REF] Liberzon | Calculus of variations and Optimal control theory[END_REF][START_REF] Bonnard | Optimal Control in Space and Quantum Dynamics[END_REF][START_REF] Jurdjevic | Geometric Control Theory[END_REF][START_REF] Boscain | Introduction to the Foundations of Quantum Optimal Control Theory[END_REF]. OCT has become nowadays a key tool in many different domains extending from space dynamics to robotics or quantum mechanics [START_REF] Bryson | Applied optimal control[END_REF][START_REF] Glaser | Training Schrödinger's cat: quantum optimal control[END_REF][START_REF] Bonnard | Optimal Control in Space and Quantum Dynamics[END_REF]. Optimal process is defined from a cost functional (to minimize) which can depend on the state of the system and the control field. For systems with complex dynamics and optimization targets which are difficult to reach, it is necessary to use optimal control algorithms converging iteratively towards the optimal Submitted to Journal Not Specified, pages 1 -21 www.mdpi.com/journal/notspecified solution. The flexibility of this approach makes it possible to adapt this tool to any experimental situation. Generally, it is possible to include constraints in the algorithms to account for requirements related to a specific material or device [START_REF] Bryson | Applied optimal control[END_REF][START_REF] Glaser | Training Schrödinger's cat: quantum optimal control[END_REF]. The only relative limitation concerns the accuracy of the modeling, even if robustness can be improved by controlling simultaneously an ensemble of systems which differ by the values of one or several constant parameters [START_REF] Li | Ensemble control of Bloch equations[END_REF][START_REF] Kozbar | Exploring the limits of broadband 90 and 180 universal rotation pulses[END_REF][START_REF] Kozbar | Exploring the limits of broadband excitation and inversion: II. Rf-power optimized pulses[END_REF][START_REF] Kozbar | Exploring the limits of broadband excitation and inversion pulses[END_REF]. OCT has been applied

to quantum systems first in the context of physical chemistry to steer chemical reactions or control specific degrees of freedom [START_REF] Brif | Control of quantum phenomena: past, present and future[END_REF][START_REF] Koch | Quantum control of molecular rotation[END_REF], followed by control of spin dynamics [START_REF] Levitt | Spin Dynamics: Basics of Nuclear Magnetic Resonance[END_REF][START_REF] Ernst | Principles of Nuclear Magnetic Resonance in one and two dimensions[END_REF] for applications in Nuclear Magnetic Resonance (NMR) [START_REF] Kozbar | Exploring the limits of broadband 90 and 180 universal rotation pulses[END_REF][START_REF] Kozbar | Exploring the limits of broadband excitation and inversion pulses[END_REF][START_REF] Mao | Selective inversion radiofrequency pulses by optimal control[END_REF][START_REF] Conolly | Optimal Control Solutions to the Magnetic Resonance Selective Excitation Problem[END_REF][START_REF] Rosenfeld | Design of adiabatic selective pulses using optimal control theory[END_REF][START_REF] Nielsen | Optimal control methods in NMR spectroscopy[END_REF][START_REF] Lapert | Singular Extremals for the Time-Optimal Control of Dissipative Spin 1/2 Particles[END_REF] and Magnetic Resonance Imaging [START_REF] Bernstein | Handbook of MRI Pulse Sequences[END_REF][START_REF] Lapert | Exploring the Physical Limits of Saturation Contrast in Magnetic Resonance Imaging[END_REF][START_REF] Vinding | Fast numerical design of spatial-selective RF pulses in MRI using Krotov and quasi-Newton based optimal control methods[END_REF][START_REF] Maximov | Real-time 2D spatially selective MRI experiments: comparative analysis of optimal control design methods[END_REF]. It has become a key tool in this domain to improve the efficiency and the sensitivity of standard experimental setups [START_REF] Glaser | Training Schrödinger's cat: quantum optimal control[END_REF]. In NMR and quantum physics, a well-known optimization method is GRAPE [START_REF] Khaneja | Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms[END_REF], which is a gradient-based algorithm [START_REF] Bryson | Applied optimal control[END_REF]. This approach has demonstrated its efficiency in many different contexts.

Fourier-Transform Ion Cyclotron Resonance (ICR) Mass Spectrometry [START_REF] Comisarow | Fourier transform ion cyclotron resonance spectroscopy[END_REF][START_REF] Comisarow | Frequency-sweep Fourier transform ion cyclotron resonance spectroscopy[END_REF] is a type of mass spectrometry based on cyclotron frequency of ions in a fixed magnetic field [START_REF] Marshall | Fourier Transform Ion Cyclotron Resonance Mass Spectrometry: A primer[END_REF][START_REF] Nikolaev | Realistic modeling of ion cloud motion in a Fourier transform ion cyclotron resonance cell by use of a particle-in-cell approach[END_REF][START_REF] Sehgal | Theory for spiralling ions for 2D FT-ICR and comparison with precessing magnetization vectors in 2D NMR[END_REF]. Ions are excited at their resonant cyclotron frequencies to larger cyclotron radii by an electric field orthogonal to the magnetic field. After the excitation pulse, the ions rotate freely with a frequency characteristic of their mass. The image current induced by the ions on a pair of electrodes is detected. The Fourier transform of the resulting transient signal leads to the mass spectrum after a proper calibration. In a homogeneous magnetic field, ICR allows to access the highest resolution available in mass spectrometry, while leading to extreme sensitivities. This spectrometry has experienced a recent renewal based on several methodological improvements and the search for very high resolutions, which are required to study complex biological or environmental mixtures. Several techniques developed and used in ICR has been inspired by equivalent approaches in NMR. An example is given by two-dimensional ICR [START_REF] Pfandler | Two-dimensional ion cyclotron resonance mass spectroscopy[END_REF][START_REF] Van Agthoven | Towards analytically useful two-dimensional Fourier transform ion cyclotron resonance mass spectrometry[END_REF][START_REF] Van Agthoven | Optimization of the discrete pulse sequence for two-dimensional FT-ICR mass spectrometry using infrared multiphoton dissociation[END_REF][START_REF] Van Agthoven | Two-dimensional FT-ICR/MS with IRMPD as fragmentation mode[END_REF][START_REF] Bray | Nonuniform Sampling Acquisition of Two-Dimensional Fourier Transform Ion Cyclotron Resonance Mass Spectrometry for Increased Mass Resolution of Tandem Mass Spectrometry Precursor Ions[END_REF][START_REF] Van Agthoven | Two dimensional mass spectrometry: new perspectives for tandem mass spectrometry[END_REF] which was proposed in analogy to two-dimensional NMR spectroscopy [START_REF] Levitt | Spin Dynamics: Basics of Nuclear Magnetic Resonance[END_REF][START_REF] Ernst | Principles of Nuclear Magnetic Resonance in one and two dimensions[END_REF]. Following this fruitful approach and given the success and efficiency of optimal control techniques in NMR, a question which naturally arises is the application of this method in ICR. This paper aims at taking a step toward the answer to this open issue. ICR Mass Spectrometry can provide very high resolution mass spectra over a large range of mass to charge ratio. In the ICR experiment, ions are initially at rest in the centre of the trap, and have to be excited in order to generate a resonant signal which can cover, in broad band experiments, frequencies from a few tens of kHz for high m/z up to 1 MHz or higher for the fastest species. However, this implies that all ions have to be excited over this frequency range in an even and controlled manner.

We explore in this study how optimal control can be used to design efficient and robust excitation pulses in ICR. To the best of our knowledge, this has never been studied yet. Due to the wide bandwidth of ICR signal, excitation pulses are usually simple chirped adiabatic pulses with a frequency sweep. Some variations have been proposed such as off-resonance monochromatic pulses for selective excitation of given ions. Based on the linearity of ion dynamics, it has also been proposed to generate pulses by Fourier synthesis from a given excitation profile, in an approach called SWIFT (for Stored-waveform Inverse Fourier Transform) [START_REF] Marshall | Tailored excitation for Fourier transform ion cyclotron resonance mass spectrometry[END_REF][START_REF] Guan | Linear response theory of ion excitation for Fourier transform mass spectrometry[END_REF][START_REF] Guan | Stored waveform inverse Fourier transform ion excitation in trapped-ion mass spectrometry: theory and applications[END_REF]. Optimal control is expected to allow a much wider range of possibilities such as the control of trajectory for given initial and final positions of the ion packet and for a given range of frequencies. In order to evaluate the contribution of OCT in ICR, we consider in this study the simplest modelling which is experimentally relevant. The experiment is considered in a simplified environment, with a uniform magnetic field and a time-dependent homogeneous electric field oriented along a single axis orthogonal to the magnetic field and with no static component. This geometry is unrealistic, as there is no trapping potential, but allows to consider the dynamics of the ions restricted to a plane with a pure cyclotron trajectory and a zero magnetron component. The time-dependent electric field aims at exciting in a robust manner an ensemble of different ions from the centre of the cell to a final position which depends in a controlled way of the ion frequency. The linearity of ion dynamics simplifies drastically the derivation of the optimal control law [START_REF] Liberzon | Calculus of variations and Optimal control theory[END_REF][START_REF] Brockett | Finite Dimensional Linear Systems[END_REF][START_REF] Bonnans | Commande et optimisation de systèmes dynamiques[END_REF]. If there is no constraint on the intensity of the electric field, linear quadratic optimal control theory (LQOCT) can be applied. Many mathematical results have been established in this case [START_REF] Liberzon | Calculus of variations and Optimal control theory[END_REF][START_REF] Li | Ensemble control of finite-dimensional time-varying linear systems[END_REF] and the optimal solution can be derived analytically. When constraints are accounted for, a numerical algorithm has to be used to solve the optimal equations. Note that very few studies have solved optimal control problems of linear systems at the microscopic scale [START_REF] Li | Control of inhomogeneous ensemble[END_REF][START_REF] Li | Ensemble control of finite-dimensional time-varying linear systems[END_REF][START_REF] Martikyan | Comparison between optimal control and shortcut to adiabaticity protocols in a linear control system[END_REF][START_REF] Martikyan | Robust control of an ensemble of springs: Application to ion cyclotron resonance and two-level quantum systems[END_REF][START_REF] Li | Exact broadband excitation of two-level systems by mapping spins to springs[END_REF]. ICR is an interesting example, relevant experimentally, to stimulate further work in this direction.

The remainder of this paper is organized as follows. The formulation of the control problem and the description of the model system are outlined in Sec. 2. After a brief introduction to the principles of OCT, we apply OCT to ICR in Sec. 3. We describe the optimal control algorithm which allows to take into account experimental constraints on the control field. Numerical results in different experimental situations are given in Sec. 4. A comparison is made with the adiabatic and the SWIFT approaches.

We conclude in Sec. 5 with an outlook and future perspectives. The Rotating Wave Approximation is discussed in Appendix A. Technical details about the adiabatic and SWIFT techniques are presented respectively in Appendices B and C. The application of LQOCT is described in Appendix D.

Formulation of the control problem

The model system

We consider the simplest modeling of ion trajectories in ICR. The different ions in the experimental cell are confined in the (x, y)plane and are subjected to a constant magnetic field B and a time-dependent electric field E respectively along the zand xaxes of the laboratory frame. Note that optimal control techniques can also be used if two control fields along the xand ydirections are available. The dynamics are governed by the Lorentz's equation:

m k ˙ v k = q k E + q k ( v k × B), (1)
where m k , q k and v k are the mass, charge and speed of the ion k. ˙ v k denotes the time derivative of v k . Equation (1) can be expressed as:

             ẋk = v x k ẏk = v y k vx k = ω k (e x + v y k ) vy k = -ω k v x k . ( 2
)
with the cyclotron frequency ω k = q k B m k and e = E/B. The coordinates (x k , y k ) and (v x k , v y k ) describe respectively the position and the speed of the ion k in the (x, y)plane. We assume that the frequency ω k belongs to the interval [ω min , ω max ] which is defined by the ion packet under study. As described below, the aim of the control process is to excite the different ions in a robust way with respect to the parameter ω.

The control problem can be defined as follows. Starting from the center of the cell (x k = 0, y k = 0) with a zero speed (v x k = 0, v y k = 0), the goal is to reach at a fixed control time t f a given radius r f and phase ϕ f . As an illustrative example, we will force the phase to vary linearly with ω, contrary to the standard result obtained with chirp pulses, where a quadratic phase dependence is observed (see Appendix B for details). We denote respectively by r k (t) and ϕ k (t) the radius and the phase of ion k at time t. We assume in a first step that there is no constraint on the electric field. A limitation on the maximum pulse intensity is accounted for in Sec. 3.2.

To simplify the notations, we omit below the index k. Using Eq. ( 2), it is straightforward to show that Ω = ωx + v y is a constant of motion. At t = 0, since x(0) = 0 and v y (0) = 0, we deduce that Ω = 0 so v y (t) = -ωx(t). One of the two coordinates v y (t) or x(t) can be eliminated. This also means that we cannot control simultaneously the position and the speed of the ion with only one control. We arrive at:

       ẏ = v y vy = -ω 2 V x Vx = v y + e x
where V x = v x /ω. We introduce the vector X = (y, v y , V x ) whose dynamics are governed by:

Ẋ = AX + Ce x , (3)
with

A =    0 1 0 0 0 -ω 2 0 1 0    , C =    0 0 1   
The dynamics of this linear system can be explicitly integrated as follows. The eigenvalues of A are (0, iω, -iω) and the corresponding eigenvectors can be written as:

X 0 =    1 0 0    , X + =    1 iω 1    , X -=    1 -iω 1   
At time t f , the state of the system is given by:

X(t f ) = t f 0 e A(t f -s) Ce x (s)ds.
We have: e At = Pe Dt P -1 , where D = diag(1, iω, -iω) and

P =    1 1 1 0 iω -iω 0 1 1    , P -1 =    1 0 -1 0 -0.5i/ω 0.5 0 0.5i/ω 0.5   
We deduce that:

e At =    1 sin(ωt)/ω -1 + cos(ωt) 0 cos(ωt) -ω sin(ωt) 0 sin(ωt)/ω cos(ωt)    and
X(t f ) = t f 0 dse x (s)    -1 + cos[ω(t f -s)] -ω sin[ω(t f -s)] cos[ω(t f -s)]    (4)

The Rotating Wave Approximation

The oscillating excitation field e x applied only along the xaxis can be expressed as the sum of two rotating fields, one in the same direction as the ions and the other in the opposite direction. We introduce the Rotating Wave Approximation (RWA) which assumes that the field rotating in opposite direction to the ions has a negligible effect on their trajectories. This approximation is verified if the range of frequencies around the central frequency ω o is not too large, as discussed in Appendix A.

Note that RWA is a standard tool in NMR [START_REF] Ernst | Principles of Nuclear Magnetic Resonance in one and two dimensions[END_REF][START_REF] Mccoy | Nonresonant effects of frequency-selectives pulses[END_REF][START_REF] Emsley | Phase shifts induced by transient Bloch-Siegert effects in NMR[END_REF] where it is derived in a similar but different way due to the non-linearity of the system [START_REF] Shahriar | Driver-phase-correlated fluctuations in the rotation of a strongly driven quantum bit[END_REF]. In particular for ICR, this approximation does not depend on the amplitude of the excitation. Using RWA, we show below that the control of ions is equivalent to the control of an ensemble of springs of different frequencies [START_REF] Li | Ensemble control of finite-dimensional time-varying linear systems[END_REF][START_REF] Li | Exact broadband excitation of two-level systems by mapping spins to springs[END_REF].

The derivation starts with the control of speeds which fulfill:

vxk = ω k v yk + ω k e x vyk = -ω k v xk
In complex coordinates, we have:

vk = -iω k v k + ω k e x (t), (5)
where

v k = v xk + iv yk . We consider that ω k ∈ [ω 0 -δω, ω 0 + δω]
where ω 0 is the carrier frequency of the electric field, e x (t) = e 0 (t) cos(ω 0 t + φ(t)), and δω is small compared to ω 0 . We also assume that the amplitude e 0 (t) and the phase φ(t) vary slowly in time with respect to the frequency ω 0 . We express the speed as: v k = ṽk e -iω 0 t , where ṽk is the complex speed in the frame rotating at frequency ω 0 . We deduce that: vk = -i∆ω k ṽk + ω k e 0 2 (e -iφ + e 2iω 0 t+iφ ),

where ∆ω k = ω kω 0 is the detuning term. In the RWA, we neglect the rapidly oscillating term exp(2iω 0 t) and we arrive at:

vk -i∆ω k ṽk + ω k e 0 2 e -iφ . (6)
It is worth noting here that in the rotating frame, the dynamics are driven by two control parameters, e 0 cos φ and e 0 sin φ.

Note that we recover here the control of an ensemble of springs. An additional step can be done for the position of the ion k,

x k = x k + iy k . We set x k = xk e -iω 0 t . It is then straightforward to show that: ẋk -iω 0 xk = ṽk (t)
Since xk varies slowly with respect to e iω 0 t , we can neglect the time derivative ẋk , which gives:

xk = i ω 0 ṽk (t).
If the RWA is valid, we deduce that the speed control leads also to the control of the position of ions.

Optimal Control Theory

A short introduction to optimal control theory

We briefly introduce in this section the tools of optimal control theory used in this study. In order to keep the introduction as accessible as possible, some mathematical details are not specified. We refer the interested reader to the specialized literature on the subject [START_REF] Bryson | Applied optimal control[END_REF][START_REF] Bressan | Introduction to the Mathematical Theory of Control[END_REF][START_REF] Liberzon | Calculus of variations and Optimal control theory[END_REF][START_REF] Boscain | Introduction to the Foundations of Quantum Optimal Control Theory[END_REF]. We consider a control system described by the following differential equation:

q(t) = f (q(t), u(t)),
where q(t) ∈ R n is the state of the system at time t, f a smooth function and u(t) ∈ R the control law. We assume here that there is no constraint on the control amplitude. The goal of a control problem is to bring the state of the system from the initial state q(0) = q 0 as close as possible to a target state q f in a time t f while minimizing a given cost functional J . For a distance to the target state defined by ||q(t f )q f ||, a standard functional is:

J = 1 2 ||q(t f ) -q f || 2 + λ t f 0 u(t) 2 dt,
where λ is a positive constant which expresses the relative weight between the distance to the target state and the second term, which can be interpreted as the energy of the control. We formulate the optimal control from the Pontryagin Maximum Principle (PMP) which gives necessary conditions for a control solution to be optimal [START_REF] Bressan | Introduction to the Mathematical Theory of Control[END_REF][START_REF] Schättler | Geometric optimal control: Theory, Methods and Examples[END_REF][START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF][START_REF] Bonnard | Optimal Control in Space and Quantum Dynamics[END_REF]. We introduce the Pontryagin Hamiltonian H P (the index P corresponds to Pontryagin) as:

H P = p(t) • f (q(t), u(t)) - λu(t) 2 2 ,
where p(t) ∈ R n is the adjoint state. This state plays qualitatively the role of a Lagrange multiplier for the optimization problem [START_REF] Liberzon | Calculus of variations and Optimal control theory[END_REF][START_REF] Boscain | Introduction to the Foundations of Quantum Optimal Control Theory[END_REF]. The state and the adjoint states of the dynamics fulfill the Hamilton's equation:

   q = ∂H P ∂p = f (q, u) ṗ = -∂H P ∂q = -p • ∂ f (q,u)
∂q with the following initial and final conditions q(0) = 0 and p(t f ) = -∂J ∂q(t f ) = q fq(t f ), while the optimal control u * is given by

∂H p ∂u = 0, i.e. u * = p λ • ∂ f (q, u * ) ∂u .
In the non-linear case, these conditions can be solved only for simple low-dimensional systems [START_REF] Schättler | Geometric optimal control: Theory, Methods and Examples[END_REF][START_REF] Bonnard | Optimal Control in Space and Quantum Dynamics[END_REF][START_REF] Lapert | Singular Extremals for the Time-Optimal Control of Dissipative Spin 1/2 Particles[END_REF] and numerical algorithms are used for more complex dynamics [START_REF] Bryson | Applied optimal control[END_REF][START_REF] Khaneja | Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms[END_REF][START_REF] Werschnik | Quantum optimal control theory[END_REF]. For linear systems, the optimal solutions can be derived explicitly if there is no additional constraint on the control field. This approach is known in the literature as Linear Quadratic Optimal Control [START_REF] Bressan | Introduction to the Mathematical Theory of Control[END_REF][START_REF] Liberzon | Calculus of variations and Optimal control theory[END_REF][START_REF] Li | Ensemble control of finite-dimensional time-varying linear systems[END_REF] and is applied to ICR in Appendix D. When experimental limitations such as maximum pulse intensity are accounted for in the numerical optimization process, the control problem becomes non-linear and the optimal control law is derived numerically from iterative algorithms, which are described in Sec. 3.2.

Optimal gradient-based algorithm

The goal of this section is to develop a first-order gradient-based algorithm suited to this control problem [START_REF] Bryson | Applied optimal control[END_REF]. We use a numerical optimization algorithm in order to take into account field amplitude constraint of the form |e x (t)| ≤ e max . This limitation makes the control problem non-linear and this latter cannot be solved analytically any more. Note that this algorithm can be seen as the counterpart of the GRAPE algorithm in NMR [START_REF] Khaneja | Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms[END_REF] and that other limitations such as spectral constraints could be added [START_REF] Lapert | Monotonically convergent optimal control theory of quantum systems with spectral constraints on the control field[END_REF]. For question of numerical stability, we apply the algorithm in the system with the RWA and then we use the derived control law in the original dynamical system.

We start from the differential system (6) written in the rotating frame for the ion k as:

v(k) x = ∆ω k ṽ(k) y + u x v(k) y = -∆ω k ṽ(k) x + u y
where u x = ω 0 2 e 0 cos φ and u y = -ω 0 2 e 0 sin φ. The two controls satisfy the limitation u x (t) 2 + u y (t) 2 ≤ u 2 max with u max = ω 0 2 e max . The corresponding target state is ( ṽ(k)

x f , ṽ(k) y f ). We consider a cost functional J with no penalty on the control field defined as:

J = 1 2 ∑ k [( ṽ(k) x f - ṽ(k) x (t f )) 2 + ( ṽ(k) y f - ṽ(k) y (t f )) 2 ]. (7)
The Pontryagin Hamiltonian can be expressed as:

H P = ∑ k [∆ω k (p (k) y ṽ(k) x -p (k) x ṽ(k) y ) + u x p (k) x + u y p (k) y ].
The adjoint states fulfill the following relations:

ṗ(k) x = -∆ω k p (k) y ṗ(k) y = ∆ω k p (k)
x .

(

) 8
and the gradients are given by:

∂H P ∂u x = ∑ k p (k)
x ,

∂H P ∂u y = ∑ k p (k) y
The correction to the control fields δu x (t) and δu y (t) at each step of the algorithm is proportional to these gradients [START_REF] Bryson | Applied optimal control[END_REF]. The final adjoint states can be expressed as:

   p (k) x (t f ) = ( ṽ(k) x f - ṽ(k) x (t f )) p (k) y (t f ) = ( ṽ(k) y f - ṽ(k) y (t f )).
and Eq. ( 8) can be directly integrated backward in time. We thus consider the following gradient-based algorithm.

1. Choose guess fields u x (t) and u y (t).

2. Propagate forward the state of every ion k and compute (v

(k) x (t f ), v (k)
y (t f )).

3. Propagate backward the adjoint state of the system from Eq. (8).

4. Compute the corrections δu x (t) and δu y (t) to the control fields, δu

x (t) = ∑ k p (k)
x , δu y (t) =

∑ k p (k)
y where is a small positive constant.

5. Define the new control fields u x (t) → u x (t) + δu x (t) u y (t) → u y (t) + δu y (t).

6. Truncate the new control fields u x (t) and u y (t) to satisfy the constraint u x (t) 2 + u y (t) 2 ≤ u max :

u x (t) → u x (t)u max u x (t) 2 + u y (t) 2 , u y (t) → u y (t)u max u x (t) 2 + u y (t) 2 .
7. Go to step 2 until a given accuracy is reached.

Similar algorithms are used in NMR for taking into account pulse constraints [START_REF] Kozbar | Exploring the limits of broadband 90 and 180 universal rotation pulses[END_REF][START_REF] Kozbar | Exploring the limits of broadband excitation and inversion: II. Rf-power optimized pulses[END_REF][START_REF] Kozbar | Exploring the limits of broadband excitation and inversion pulses[END_REF]. The efficiency of such algorithms has been shown in different domains. The use of a gradient causes this type of algorithm to converge towards a local maximum of the optimization problem. Numerical simulations with different guess fields allow to partly overcome this limitation, even if the global maximum is not reached with certainty. The efficiency of this algorithm in ICR is illustrated numerically in Sec. 4.

Numerical results

We present numerical results obtained either with LQOCT (see Appendix D for details) or with the gradient-based algorithm. A comparison with the SWIFT approach described in Appendix C is also done. Different experimental constraints have to be satisfied by the control pulse. The objective is to excite ions in a wide range of frequencies around a central frequency of the order of 500 kHz. The excitation has to be as uniform as possible in radius and phase in the range [ f min , f max ] and close to zero outside. As a benchmark example, we choose in this section to consider the interval [400,600] kHz. Using the linearity of the dynamics, these results can be transposed to another range of frequencies by a scaling of the excitation pulse duration and of the pulse amplitude. For instance, if the total process time is increased by a factor α, then the range of frequencies and the amplitude of the electric field are divided by the same parameter α. The description of the optimal control of this infinite dimensional dynamical system is mathematically quite intricate, even if some results can be established [START_REF] Li | Ensemble control of finite-dimensional time-varying linear systems[END_REF][START_REF] Martikyan | Robust control of an ensemble of springs: Application to ion cyclotron resonance and two-level quantum systems[END_REF]. For practical and numerical reasons, it is more convenient to consider a finite set of N systems by discretizing the frequency interval. We thus consider the simultaneous control of each element of this set. In this paper, we consider a regular discretization, but other choices could be possible, and the frequency step is chosen small enough to avoid the discretization effect. Note that the same approach is used in NMR to control a spin ensemble [START_REF] Kozbar | Exploring the limits of broadband 90 and 180 universal rotation pulses[END_REF][START_REF] Kozbar | Exploring the limits of broadband excitation and inversion: II. Rf-power optimized pulses[END_REF][START_REF] Kozbar | Exploring the limits of broadband excitation and inversion pulses[END_REF]. The required sharp excitation profile is modeled by the following function:

r( f ) = r 0 2 [tanh(µ( f -f min )) + tanh(µ( f max -f ))],
where µ is a free parameter allowing to adjust the slope of the excitation gate. In the numerical simulations, µ is fixed to 0.1. We impose that the final radius is r f = 5 cm and a final phase varying linearly with the frequency, with a maximum variation of the order of 1 o /1 Hz, which corresponds to 17.45 rad/kHz. The final phase ϕ

(k) f of the ion k is expressed as ϕ (k) f = -aω k t f
, where a is a parameter characterizing the slope of the angular variation. The magnetic field is set to 7 T and the maximum electric field amplitude that can be generated is of the order of 10 3 V.m -1 . The control time can be very long, of the order of few hundred ms, but more stable numerical results were achieved for duration of the order of few ms.

We first present in Fig. 1 a series of simulations without limitation on the field intensity. The excitation pulse duration is chosen to be equal to 1 ms. The electric field is computed from a set of the order of 500 frequencies regularly spaced in the interval under study. Since the derived solutions are very sharp, this duration can be modified to some extent without changing the control pulse. Figure 1 compares the results achieved by LQOCT and by the SWIFT approach. The optimal solutions can be computed by using or not the RWA. Note that the pulse computed in the RWA is then applied without any approximation to the original system. In the case displayed in Fig. 1, very similar efficiencies are obtained for the two optimal excitations. The optimal pulses are very similar to the SWIFT solution, even if the analytical expressions of the pulses are different. This point is not so surprising since for a linear system and a continuous range of frequencies, the control field is expected to be unique.

This statement can be rigorously shown in the case of an ensemble of springs under some specific mathematical assumptions [START_REF] Li | Ensemble control of finite-dimensional time-varying linear systems[END_REF][START_REF] Martikyan | Robust control of an ensemble of springs: Application to ion cyclotron resonance and two-level quantum systems[END_REF]. However, the optimal control method offers greater flexibility since one can play with different parameters such as the cost functional or the number of discrete frequencies to adjust the final result. We study in Fig. 2 the role of the phase slope of the excitation profile on the structure of the pulse. Figure 2 shows that this slope changes the position of the peak of the pulse. This position can be deduced from a Fourier transform of the profile. Very good results have been obtained

for slopes in the range [0.05,0.95] with a maximum pulse amplitude almost constant. Pulse distortion appears when the slope parameter a is close to 0 or 1. For a = 0, it becomes very difficult to control all the ions which have to reach a fixed target state in a space-fixed frame, independently of their own frequency. Note that similar results were achieved in NMR [START_REF] Gershenzon | Linear phase slope in pulse design: Application to coherence transfer[END_REF][START_REF] Gershenzon | Optimal control design of excitation pulses that accommodate relaxation[END_REF], which highlights the similarities between the control of the two dynamics.

We now focus on ion control with amplitude constraint. The numerical simulations have been carried out by assuming the RWA. The same set of discretized frequencies is chosen. We optimize piecewise constant functions with a time step lower than 1 µs to avoid discretization effect. The dynamics are integrated numerically through the formulas given in Sec. 2. More than 1000 iterations are usually needed to converge to an efficient solution. In a first step, we apply the gradient-based algorithm described in Sec. 3.2 with only one control field, namely E 0 (t) = e 0 (t)B and the phase φ(t) of the electric field is set to 0. We consider the same control problem as before and the optimal solutions derived above are used as guess field for the optimization algorithm. Comparison between the optimal and the SWIFT approaches for the robust control of an ensemble of ions in the frequency range [400, 600] kHz. The small insert is a zoom of the profile around the frequency f = 600 kHz. Panels (a) and (e) display the evolution of the final radius and phase as a function of f . The black, blue (dark gray) and red (light gray) solid lines depict respectively the optimal solutions computed without and with the RWA and the SWIFT pulse. The SWIFT and optimal control laws are plotted in panels (b), (c) (optimal without RWA) and (d) (optimal with RWA). The number of discretized frequency points is set to 601 in the optimization process in the range [350,650] kHz. 100 V.m -1 , while maintaining an almost ideal excitation profile. This reduction was made possible by a distribution of the energy along the control interval. Outside of t 0.5 ms, the amplitude of the optimized field is much larger than the one of the guess pulse. As a comparison, Fig. 3 also presents the profile obtained from the optimal pulse whose amplitude has been arbitrarily limited to 100 V.m -1 , showing the non-trivial transformation made by the algorithm.

The optimization algorithm fails to converge towards a very good excitation profile, when the maximum amplitude is much smaller than 100 V.m -1 . This obstacle can be partly overcome by considering in a second step two control fields (in the rotating frame) denoted E 0x = e 0 B cos φ and E 0y = e 0 B sin φ. An example is displayed in Fig. 4 for a maximum amplitude of 100 and 50 V.m -1 . An almost perfect excitation profile is achieved in these two cases. Note the different structures of the fields along the xand ydirections, namely even and odd functions. This observation was also made in some optimal control problems in NMR [START_REF] Kozbar | Exploring the limits of broadband excitation and inversion pulses[END_REF].

A systematic analysis of the efficiency of the optimized control fields with respect to the maximum pulse amplitude and to the control duration is provided in Fig. 5. The efficiency of the control process is measured from the cost functional J given in Eq. ( 7). As could be expected, we observe that better results are achieved for larger maximum amplitude and control time. However, the final fidelity varies in a quite complex way with the control time. A saturation is observed for times of the order of few ms. It is not clear if this point is due to an intrinsic limitation of the control protocol or to convergence problems of the algorithm. Further investigations are needed to clarify this issue.

Conclusions

We have applied optimal control techniques to the robust excitation of ions in ICR. We have considered the simplified but realistic conditions of a two-dimensional trajectory and of a homogeneous magnetic field. In this model system, we propose different ways to solve the optimal control problems. Such methods are directly inspired from NMR in which OCT is a standard and efficient tool. In the case without pulse limitation, the linearity of the dynamical equations allows to use LQOCT, which has the advantage to lead to an analytical formula of the control law. Very good results have been obtained both for the final radii and phases of the ions. A specific range of frequencies has been considered in this paper, but the same approach can be extended to broadband excitation from 100 to 900 kHz. However, this solution is both in shape and in amplitude very similar to the SWIFT pulse. The two solutions are expected to be equal for a continuous range of frequencies. More original control laws are derived when the pulse intensity is limited. Due to this constraint, optimal iterative algorithms have to be used and we adapt to ICR the standard GRAPE algorithm, well-known in NMR. Even if this algorithm has some limitations, it allows to reduce the pulse intensity, by a factor larger than three in the examples under study. On the basis of NMR results, this algorithm is expected to be very efficient in the case of other excitation profiles. The very encouraging and promising results obtained in this investigation must now be confirmed by experimental implementation. Numerical simulations of this study are not fully realistic. Effects such as the magnetron motion, field geometry, field inhomogeneities or ion collisions are neglected. However, the model system we consider describes quite faithfully the main cyclotronic behavior and permits to grasp rapidly the main features of ion trajectories. Numerical codes have been developed to account for such experimental details. The relative simplicity of the application of numerical optimal algorithms makes it possible to adapt it straightforwardly to a new class of control problems. They could thus be combined with such codes.

We are therefore quite confident about the extension of optimization procedures to these additional experimental constraints and limitations. Work is in progress on these different issues.

where φ is a smooth function, which is assumed to be rapidly varying with respect to h. A stationary point t 0 is defined by φ (1) (t 0 ) = 0, where φ (n) denotes the nth time derivative of φ. A Taylor expansion around t = t 0 leads to:

φ(t) = φ(t 0 ) + (t -t 0 )φ (1) (t 0 ) + (t -t 0 ) 2 2 φ (2) (t 0 ) + • • •
We arrive at: ĥ(ω)

h(t 0 )e iφ(t 0 ) +∞ -∞ e i ξ 2 2 φ (2) (t 0 ) dξ 2π φ (2) (t 0 ) h(t 0 )e i(φ(t 0 )+ π 4 ) .
For a chirp excitation, the phase φ(t) is defined by φ(t) = ω i t + st 2 2 . The instantaneous frequency ω(t) can be expressed as:

ω(t) = φ(t) = ω i + st,
where s = ω(t). In the example under study, the rate s is given by s = (ω fω i )/t f . We assume that s > 0 and we deduce that the Fourier transform of the control field is given by: êx (ω) =

t f 0 e x (t)e -iωt dt = e 0 2 t f 0 [e i(ω i t+ st 2 2 -ωt) + e -i(ω i t+ st 2 2 +ωt) ]dt.
We denote by φ 1 and φ 2 the arguments of the two exponential terms. It is straightforward to verify that φ1 (t) = 0 for t = t (ω) 1

= ω-ω i s and that φ2 (t) = 0 for t = t < 0 and assuming that t

(ω) 1

is not too close to 0 and t f , we can consider that the integral is defined from -∞ to +∞. We finally get: êx (ω) = e 0 π 2s e i( π

4 +φ 1 (t (ω) 1
)) .

The phase spectrum φ(ω) = π 4 + φ 1 (t (ω)

1 ) can be written as:

φ(ω) = π 4 - (ω -ω i ) 2 2s .
Coming back to the original control problem, we obtain:

v k (t f ) e 0 π 2s exp[i( 3π 4 -ωt f + (ω k -ω i ) 2 2s
)].

In the range of validity of this approximation, we observe that the final radius of ions at time t = t f is a constant, while the phase varies quadratically with the frequency ω. A numerical example is given in Fig. 5. The frequency of the chirped pulse goes from 400 kHz to 600 kHz.

Appendix C. Excitation by the SWIFT approach

In this paragraph, we describe the application of the SWIFT method to the model system. We consider a specific approach in which the control law and the corresponding dynamics can be expressed analytically. The dynamics are governed by the differential system (2). In the RWA described in Sec. 2.2, the dynamics can be approximated as:

vk = -i∆ω k ṽk + ω 0 2 e 0 e -iφ ,
where ṽk = ṽxk + i ṽyk and the control field is expressed as e x (t) = e 0 (t) cos(ω 0 t + φ(t)). The differential equation can be integrated and leads to:

ṽk (t f ) = t f 0 e -i∆ω k (t f -t) ω 0 2 e 0 e -iφ dt
We deduce that:

ṽ * k (t f )e -i∆ω k t f = t f 0 e -i∆ω k t ω 0 2
e 0 e iφ dt.

Introducing u(t) = e 0 e iφ and assuming that u is different from zero only in the interval [0, t f ], we obtain:

√ 2π ω 0 2 û(∆ω k ) = ṽ * k (t f )e -i∆ω k t f .
where we use the following definition for the Fourier transform:

f (t) = 1 √ 2π +∞ -∞ f (ω)e iωt dω; f (ω) = 1 √ 2π +∞ -∞ f (t)e -iωt dt.
The target states are defined as:

r ∆ω k = r 0 Π( ∆ω k δω ) φ ∆ω k = a∆ω k + φ 0
where Π is the gate function, with Π(x) = 1 if |x| ≤ 1 2 and 0 otherwise. The parameter δω is the width of the distribution and φ 0 is an arbitrary constant. We have:

x k = r ∆ω k e i(a∆ω k +φ 0 ) .
In the RWA, starting from ṽk = -iω 0 xk , we arrive at: ṽk (t) = -iω 0 r ∆ω k e i(a∆ω k +φ 0 ) e iω 0 t and û(∆ω) = 2i √ 2π r ∆ω k e -i∆ω k (t f +a) e -iφ 0 -iω 0 t f , which gives

u(t) = FT -1 [ 2i √ 2π r ∆ω k e -i∆ω k (t f +a) e -iφ 0 -iω 0 t f ]. Since 1 √ 2π +∞ -∞ Π( ω δω )e iωt dω = δω √ 2π sinc( δωt 2 ),
we obtain:

u(t) = r 0 δωe iφ 1 π sinc[ δω 2 (t -t 0 )],
with t 0 = t f + a and φ 1 an arbitrary phase. The original control field e x (t) = e 0 (t) cos(ω 0 t + φ(t)) is then given by:

e x (t) = r 0 δω π sinc( δω 2 (t -t 0 )) cos(ω 0 t + φ 1 ).
Since the choice of the initial phase φ 1 is arbitrary, we finally get:

e x (t) = r 0 δω π sinc( δω 2 (t -t 0 )) cos(ω 0 (t f -t)).
The next step consists in integrating exactly the system dynamics using the original system and Eq. ( 4).

We have to compute terms of the form:

I c (t, ω s , ω) = sinc(ω s (t -t 0 )) cos(ω(t f -t))dt I s (t, ω s , ω) = sinc(ω s (t -t 0 )) sin(ω(t f -t))dt.
For that purpose, we use the Sine and the Cosine Integral functions Si and Ci which are defined by: Si

(x) = x 0 sinc(t)dt, Ci(x) = - ∞ x cos t t dt, x > 0.
We have the following results:

261 I c (t, ω s , ω) = sin[(t 0 -t f )ω] 2ω s Ci[(t 0 -t)(ω + ω s )] -Ci[(t 0 -t)(ω -ω s )] + cos[(t 0 -t f )ω] 2ω s Si[(t 0 -t)(ω -ω s )] -Si[(t 0 -t)(ω + ω s )] 262 I s (t, ω s , ω) = cos[(t 0 -t f )ω] 2ω s Ci[(t 0 -t)(ω + ω s )] -Ci[(t 0 -t)(ω -ω s )] + sin[(t 0 -t f )ω] 2ω s Si[(t 0 -t)(ω + ω s )] -Si[(t 0 -t)(ω -ω s )]
The final state of the dynamics is given by the following expressions:

x ω (t f ) = t f 0 dte x (t) sin[ω(t f -t)] y ω (t f ) = t f 0 dte x (t)(-1 + cos[ω(t f -t)]).
We then deduce:

x ω (t f ) = r 0 δω 2π [I s (t f , ω s , ω 0 + ω) -I s (0, ω s , ω 0 + ω) -I s (t f , ω s , ω 0 -ω) + I s (0, ω s , ω 0 -ω)]
and

y ω (t f ) = r 0 δω 2π [-2I c (t f , ω s , ω 0 ) + 2I c (0, ω s , ω 0 ) +I c (t f , ω s , ω 0 + ω) -I c (0, ω s , ω 0 + ω) + I c (t f , ω s , ω 0 -ω) -I c (0, ω s , ω 0 -ω)]
The results achieved with this approach are described in Sec. 4.

Appendix D. Application of LQOCT to ICR

We apply in this section the PMP to ICR processes in the case without any amplitude constraint. We denote by X k the state associated with the frequency ω k as defined in Eq. (3) of Sec. 2.1 and by (X

(k) 1 , X (k) 2 , X (k)
3 ) the coordinates. {ω k } is the set of discrete frequencies used in the numerical optimization. The optimal problem is defined through the cost functional J to minimize:

J = 1 2 ∑ k [(X (k) 1 (t f ) -X (k) 1 f ) 2 + (X (k) 2 (t f ) -X (k) 2 f ) 2 ] + λ 2 t f 0 e 2 x dt.
Since there is no final condition on X 3 (t f ), this term does not appear in the expression of J . The Pontryagin Hamiltonian is given by:

H P = ∑ k [p (k) 1 X (k) 2 -ω 2 k p (k) 2 X (k) 3 + p (k) 3 X (k) 2 + p (k) 3 u] - λ 2 e 2 x .
For the adjoint state, we have:

       ṗ(k) 1 = 0 ṗ(k) 2 = -p (k) 1 -p (k) 3 ṗ(k) 3 = ω 2 k p (k) 2
with the final conditions:

       p (k) 1 (t f ) = X (k) 1 f -X (k) 1 (t f ) p (k) 2 (t f ) = X (k) 2 f -X (k) 2 (t f ) p (k) 3 (t f ) = 0 (A2) Note that p (k)
1 is a constant of the motion. We deduce the dynamics of the adjoint state:

       p (k) 1 (t) = p (k) 1 (t f ) p (k) 2 (t) = A (k) cos(ω k t) + B (k) sin(ω k t) p (k) 3 (t) = -p (k) 1 (t f ) + ω k [A (k) sin(ω k t) -B (k) cos(ω k t)] (A3) with      A (k) = sin(ω k t f ) p (k) 1 (t f ) ω k + p (k) 2 (t f ) cos(ω k t f ) B (k) = sin(ω k t f )p (k) 2 (t f ) - p (k) 1 (t f ) ω k cos(ω k t f )
The optimal control e *

x can be expressed as:

e * x (t) = 1 λ ∑ k p (k) 3 (t) (A4)
which can be transformed into:

e * x (t) = 1 λ ∑ k [-p (k) 1 (t f ) + p (k) 1 (t f ) cos(ω k (t f -t)) -p (k)
The last step consists in computing the trajectory corresponding to this optimal control field. We obtain for an ion of frequency ω:

X 1 (t f ) = 1 λ ∑ k [p (k) 1 (t f )(t f - sin(ωt f ) ω - sin(ω k t f ) ω k ) + p (k) 2 (t f )(1 -cos(ω k t f ) + p (k) 1 (t f ) 2 [ sin((ω k + ω)t f ) ω k + ω + sin((ω k -ω)t f ) ω k -ω ] + ω k p (k) 2 (t f ) 2 [ cos((ω k + ω)t f ) -1 ω k + ω + cos((ω k -ω)t f ) -1 ω k -ω ]]
and

X 2 (t f ) = -ω λ ∑ k [p (k) 1 (t f ) cos(ωt f ) -1 ω + p (k) 1 (t f ) 2 [ 1 -cos((ω k + ω)t f ) ω k + ω + cos((ω k -ω)t f ) -1 ω k -ω ] + ω k p (k) 2 (t f ) 2 [ sin((ω k + ω)t f ) ω k + ω - sin((ω k -ω)t f ) ω k -ω ]]
Such results can be written in a compact form as follows:

λX (j) 1 (t f ) = ∑ k [R jk p (k) 1 (t f ) + S jk p (k) 2 (t f )] λX (j) 2 (t f ) = ∑ k [T jk p (k) 1 (t f ) + U jk p (k) 2 (t f )]
where the matrices R, S, T and U are known explicitly and the index j labels the ion of the ensemble. We finally arrive at the following system to fulfill:

   ∑ k R jk X (k) 1 f + S jk X (k) 2 f = λX (j) 1 (t f ) + ∑ k R jk X (k) 1 (t f ) + S jk X (k) 2 (t f ) ∑ k T jk X (k) 1 f + U jk X (k) 2 f = λX (j) 2 (t f ) + ∑ k T jk X (k) 1 (t f ) + U jk X (k)

(t f )

In matrix form, for N = 2, we have:

       ∑ k R 1k X (k) 1 f + S 1k X (k) 2 f ∑ k T 1k X (k) 1 f + U 1k X (k) 2 f ∑ k R 2k X (k) 1 f + S 2k X (k) 2 f ∑ k T 2k X (k) 1 f + U 2k X (k) 2 f        =      λ +
           X (1) 1 (t f ) X (1) 2 (t f ) X (2) 1 (t f ) X (2) 2 (t f )      
This linear system allows to compute the final state of the system X k (t f ), then the adjoint state from Eq. (A2) and (A3) and the optimal control field with Eq. (A4). We observe that the control law is expressed as a linear combination of cosine and sine functions of the frequencies ω k of the finite discretized set. Numerical results are presented in Sec. 4. Note that the same method can be applied in the RWA starting from Eq. ( 6) (see Ref. [START_REF] Martikyan | Robust control of an ensemble of springs: Application to ion cyclotron resonance and two-level quantum systems[END_REF] for details).

Figure 3

 3 displays the best result achieved with this limitation. The maximum field amplitude can be reduced from 130 to

Figure 1 .

 1 Figure 1.Comparison between the optimal and the SWIFT approaches for the robust control of an ensemble of ions in the frequency range [400, 600] kHz. The small insert is a zoom of the profile around the frequency f = 600 kHz. Panels (a) and (e) display the evolution of the final radius and phase as a function of f . The black, blue (dark gray) and red (light gray) solid lines depict respectively the optimal solutions computed without and with the RWA and the SWIFT pulse. The SWIFT and optimal control laws are plotted in panels (b), (c) (optimal without RWA) and (d) (optimal with RWA). The number of discretized frequency points is set to 601 in the optimization process in the range[350, 650] kHz.

Figure 2 .

 2 Figure 2. Same as Fig. 1 but for different slopes of the excitation profile. The parameter a is fixed respectively to 0.25 and 0.75 in panels (a) and (b).

Figure 3 .

 3 Figure3. Panel (a) displays the final radius r as a function of the frequency f . The black, blue (dark gray) and red (light gray) curves represent respectively the ideal profile, the one obtained with the optimization algorithm and the one corresponding to a pulse whose amplitude has been abruptly limited (see the text for details). The amplitude E 0 in the rotating frame of the optimal fields with (black curve) and without (red or light gray curve) constraints are depicted in panel (b).

Figure 4 .Figure 5 .

 45 Figure 4. Plot of the optimal amplitudes E 0x (blue or dark gray) and E 0y (red or light gray) for a maximum amplitude of 100 (panel (a)) and 50 V.m -1 (panel (b)). Panel (c) represents the corresponding total amplitude E 0 = E 2 0x + E 2 0y .

Figure A1 .

 A1 Figure A1. Excitation of an ensemble of ions by an adiabatic pulse: Evolution of the final radius as a function of the frequency. The parameters are set to t f = 10 ms, E 0 =3.2 V/m and B 0 =7 T. The red (light gray) solid line represents the stationary phase approximation. The vertical blue (dark gray) solid lines indicate the range of frequency of the pulse.

Funding: This research project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie-Sklodowska-Curie Grant Agreement No. 765267 (QUSCO). S.J.G. acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy, Grant No.EXC-2111-390814868. C. B. and M. A. D. acknowledge the European Project EU FT-ICR MS (H2020-INFRAIA-02-2017.

Author Contributions: V. M. and D. S. carried out the numerical simulations. M. A. D. proposed the initial idea about the application of OCT to ICR. M. A. D., S. G. and D. S. worked on the design of the research project. M. A. D. and C. M. wrote the initial numerical code to simulate the dynamics. The connection between ICR and NMR was analyzed by S. J. G. and D. S. All authors participated in the analysis of the results and in the paper redaction.

Conflicts of Interest:

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

Appendix A. The Rotating Wave Approximation

We discuss in this section the validity of the Rotating Wave Approximation described in Sec. 2.2. We consider the following dynamical system:

which corresponds to Eq. ( 5) of the main text. Equation (A1) describes a spring of frequency ω/(2π) excited by an external field of constant amplitude u and of frequency ω 0 /(2π). Introducing the frame rotating at ω 0 with the transformation z = ze -iω 0 t , we arrive at:

where ∆ω = ωω 0 is the detuning. In the RWA, we neglect the fast oscillating term and we get:

where zr denotes the approximate z-variable. We set δ z = zzr and we obtain:

This differential system can be exactly integrated:

This leads to:

We deduce that the relative error due to the RWA can be expressed as:

A rough approximation gives:

Numerical simulations show that this formula overestimates the error and that RWA can be used in a quite wide interval around the carrier frequency of the excitation pulse.

Appendix B. Adiabatic Excitation of ICR Process

The goal of this paragraph is to compute the final states of the ions in the case of an adiabatic excitation of the form e x = e 0 cos(ω i t + s 2 t 2 ) where ω i is the initial frequency and s the sweep rate. We recall that integrals of the form:

can be computed from the Erfi function. This result allows to compute exactly the dynamics of the system. Starting from Eq. ( 4), the final state of the ICR process can be expressed as follows:

and we finally obtain:

An approximation of the dynamics can be derived by using the stationary phase approximation. For that purpose, we start from Eq. ( 5) and we assume that t f 0 e x (t)dt = 0. We have:

The stationary phase approximation can be stated as follows. We consider the following integral: ĥ(ω) =


Marc-André Delsuc
email: madelsuc@unistra.fr

Kathrin Breuker
email: kathrin.breuker@uibk.ac.at

Maria A Van Agthoven
email: maria.van-agthoven@uibk.ac.at

Phase Correction

Keywords: mass spectrometry, two-dimensional, tandem mass spectrometry, phase correction, absorption mode, data processing

published or not. The documents may come

Introduction

Two-dimensional mass spectrometry (2D MS) is a tandem mass spectrometry method that relies on manipulating ion motions to correlate precursor and fragment ion signals [START_REF] Pfändler | Two-dimensional Fourier transform ion cyclotron resonance mass spectrometry[END_REF][START_REF] Pfändler | Broad-band two-dimensional Fourier transform ion cyclotron resonance[END_REF][START_REF] Bensimon | A method to generate phase continuity in two-dimensional Fourier transform ion cyclotron resonance mass spectrometry[END_REF]. Various 2D MS methods have been developed for quadrupolar ion traps, linear ion traps, and Fourier transform ion cyclotron resonance mass spectrometry (FT-ICR MS) [START_REF] Snyder | Logical MS/MS scans: A new set of operations for tandem mass spectrometry[END_REF][START_REF] Snyder | Two-Dimensional Tandem Mass Spectrometry in a Single Scan on a Linear Quadrupole Ion Trap[END_REF][START_REF] Szalwinski | MS Spectra Recorded in the Time Domain Using Repetitive Frequency Sweeps in Linear Quadrupole Ion Traps[END_REF][START_REF] Szalwinski | Complex mixture analysis by two-dimensional mass spectrometry using a miniature ion trap[END_REF][START_REF] Van Agthoven | Two-dimensional mass spectrometry in a linear ion trap, an in silico model[END_REF]. In addition to precursor-fragment correlations, two-dimensional partial correlation mass spectrometry (2D PC MS) has also been developed to explore correlations between complementary fragment ion signals [START_REF] Driver | Two-Dimensional Partial-Covariance Mass Spectrometry of Large Molecules Based on Fragment Correlations[END_REF]. In 2D FT-ICR MS, the most frequently used pulse sequence is the one proposed by Pfändler et al. to modulate ion radii in the ICR cell before fragmentation [START_REF] Pfändler | Two-dimensional Fourier transform ion cyclotron resonance mass spectrometry[END_REF][START_REF] Pfändler | Broad-band two-dimensional Fourier transform ion cyclotron resonance[END_REF][START_REF] Bensimon | A method to generate phase continuity in two-dimensional Fourier transform ion cyclotron resonance mass spectrometry[END_REF][START_REF] Van Agthoven | Towards analytically useful two-dimensional Fourier transform ion cyclotron resonance mass spectrometry[END_REF][START_REF] Van Agthoven | Two-dimensional mass spectrometry: New perspectives for tandem mass spectrometry[END_REF]. The pulse sequence and data processing have been optimized for analytical use [START_REF] Van Agthoven | Two-dimensional Fourier transform ion cyclotron resonance mass spectrometry: Reduction of scintillation noise using Cadzow data processing[END_REF][START_REF] Chiron | Efficient denoising algorithms for large experimental datasets and their applications in Fourier transform ion cyclotron resonance mass spectrometry[END_REF][START_REF] Van Agthoven | Optimization of the discrete pulse sequence for two-dimensional FT-ICR mass spectrometry using infrared multiphoton dissociation[END_REF][START_REF] Bray | Nonuniform Sampling Acquisition of Two-Dimensional Fourier Transform Ion Cyclotron Resonance Mass Spectrometry for Increased Mass Resolution of Tandem Mass Spectrometry Precursor Ions[END_REF]. 2D FT-ICR MS has been used for applications ranging from small molecules and agrochemicals, to polymers, bottomup and top-down proteomics [START_REF] Van Agthoven | Two-dimensional FT-ICR/MS with IRMPD as fragmentation mode[END_REF][START_REF] Van Agthoven | Two-Dimensional ECD FT-ICR Mass Spectrometry of Peptides and Glycopeptides[END_REF][START_REF] Van Agthoven | Differentiating Fragmentation Pathways of Cholesterol by Two-Dimensional Fourier Transform Ion Cyclotron Resonance Mass Spectrometry[END_REF][START_REF] Floris | 2D FT-ICR MS of Calmodulin: A Top-Down and Bottom-Up Approach[END_REF][START_REF] Van Agthoven | Two-Dimensional Mass Spectrometry for Proteomics, a Comparative Study with Cytochrome c[END_REF][START_REF] Floris | Bottom-Up Two-Dimensional Electron-Capture Dissociation Mass Spectrometry of Calmodulin[END_REF][START_REF] Van Agthoven | Can Two-Dimensional IR-ECD Mass Spectrometry Improve Peptide de Novo Sequencing?[END_REF][START_REF] Paris | Facile Determination of Phosphorylation Sites in Peptides Using Two-Dimensional Mass Spectrometry[END_REF][START_REF] Marzullo | Advantages of Two-Dimensional Electron-Induced Dissociation and Infrared Multiphoton Dissociation Mass Spectrometry for the Analysis of Agrochemicals[END_REF]. 2D MS has also been adapted for high-resolution analysis of precursor ions over narrower mass-to-charge ratio (m/z) ranges [START_REF] Halper | Narrowband Modulation Two-Dimensional Mass Spectrometry and Label-Free Relative Quantification of Histone Peptides[END_REF].

Since 2D MS relies on Fourier transformation in both dimensions, one area in which its performance can be improved is the phase correction for absorption mode mass spectrometry. Until now, 2D mass spectra have been shown in magnitude mode. In onedimensional FT-ICR MS, phase correction for absorption mode was proposed by Xian et al. and by Qi et al. [START_REF] Xian | Automated Broadband Phase Correction of Fourier Transform Ion Cyclotron Resonance Mass Spectra[END_REF][START_REF] Qi | Phase correction of Fourier transform ion cyclotron resonance mass spectra using MatLab[END_REF]. Automation for phase correction in one-dimensional FT-ICR MS was achieved by Kilgour et al. [START_REF] Kilgour | An Algorithm for Automated Generation of Absorption Mode Spectra for FT-ICR MS[END_REF][START_REF] Kilgour | Producing absorption mode Fourier transform ion cyclotron resonance mass spectra with non-quadratic phase correction functions[END_REF][START_REF] Kilgour | Producing Isotopic Distribution Models for Fully Apodized Absorption Mode FT-MS[END_REF]. Processing data sets in absorption mode has improved the performance of one-dimensional FT-ICR MS by a factor of √2≅ 1.41 in signal-to-noise and a factor of 2 in resolving power, resulting in more complete and reliable analytical information for complex samples [START_REF] Qi | Absorption-Mode Fourier Transform Mass Spectrometry: The Effects of Apodization and Phasing on Modified Protein Spectra[END_REF][START_REF] Smith | Absorption mode FTICR mass spectrometry imaging[END_REF][START_REF] Kilgour | Improved optimization of the Fourier transform ion cyclotron resonance mass spectrometry phase correction function using a genetic algorithm[END_REF][START_REF] Qi | Absorption-mode spectra on the dynamically harmonized Fourier transform ion cyclotron resonance cell[END_REF][START_REF] Qi | Absorption-Mode: The Next Generation of Fourier Transform Mass Spectra[END_REF][START_REF] Cho | Application of Phase Correction to Improve the Interpretation of Crude Oil Spectra Obtained Using 7 T Fourier Transform Ion Cyclotron Resonance Mass Spectrometry[END_REF][START_REF] Qi | Application of phase correction to improve the characterization of photooxidation products of lignin using 7 Tesla Fourier-transform ion cyclotron resonance mass spectrometry[END_REF]. In 2D MS, we can therefore expect an improvement by a factor of √2× √2=2 in signal-to-noise and in resolving power in each dimension.

A previous study has shown that the phase correction function is quadratic in the horizontal fragment ion dimension and linear in the vertical precursor ion dimension [START_REF] Van Agthoven | Phase relationships in two-dimensional mass spectrometry[END_REF]. Here, we demonstrate how to perform phase correction on a narrowband 2D mass spectrum and we show how absorption mode 2D MS compares to magnitude mode.

Theory

The pulse sequence for 2D MS is shown in Figure 1. The chirp pulses used for the excitation are needed to excite a large range of frequencies for the fragment ions and they induce a quadratic phase dependence, requiring a quadratic phase correction in the fragment ion dimension [START_REF] Xian | Automated Broadband Phase Correction of Fourier Transform Ion Cyclotron Resonance Mass Spectra[END_REF][START_REF] Qi | Phase correction of Fourier transform ion cyclotron resonance mass spectra using MatLab[END_REF][START_REF] Comisarow | Frequency-sweep Fourier transform ion cyclotron resonance spectroscopy[END_REF]. In Figure 1, the quadratic dependence of the phase of ion signals on their resonant cyclotron frequencies is expressed through the parameters α for fragment ions and β for precursor ions. In Bruker FT-ICR mass spectrometers, the initial phase of each pulse in the pulse sequence rotates with the minimum frequency of the pulses f min . Figure 1 shows the initial phase of each pulse, which leads to a corresponding phase offset for the signals of ions that are excited from the center of the ICR cell in that pulse. During the first pulse, the precursor ions are excited from the center of the ICR cell to higher radius and start accruing phase [START_REF] Van Agthoven | Phase relationships in two-dimensional mass spectrometry[END_REF]. The initial phase of the first pulse is zero. At the end of the first pulse, the precursor ions have been excited to a radius r 0 . The precursor ions rotate during the encoding delay t 1 at their own cyclotron frequency. The second pulse has a starting phase of 2πf min (T 1 + t 1 ). The phase difference between the second pulse and the precursor ion motion is 2πt 1 (f ICR -f min ) where f ICR is the cyclotron frequency of the ion, f min is the minimum frequency in the pulse, and t 1 is the encoding delay. Therefore, depending on the phase difference, precursors are excited to higher radii or de-excited towards the center of the ICR cell. The frequency of the precursor ion radius modulation is f ICR -f min . At the end of the second pulse, the radius of the precursor ions can be expressed as [START_REF] Pfändler | Broad-band two-dimensional Fourier transform ion cyclotron resonance[END_REF][START_REF] Van Agthoven | Two-dimensional mass spectrometry: New perspectives for tandem mass spectrometry[END_REF][START_REF] Guan | A theory for two-dimensional Fourier-transform ion cyclotron resonance mass spectrometry[END_REF]:

r ( t 1 ) =r 0 √ 2 ( 1+ cos 2 π ( f ICR -f min ) ( t 1 -T 1 )) (1)
After the second pulse, a fragmentation period τ m with a method that has a maximum fragmentation efficiency at the center of the ICR cell (e.g., ECD or IRMPD) is applied to the precursor ions [START_REF] Van Agthoven | Optimization of the discrete pulse sequence for two-dimensional FT-ICR mass spectrometry using infrared multiphoton dissociation[END_REF][START_REF] Van Agthoven | Can Two-Dimensional IR-ECD Mass Spectrometry Improve Peptide de Novo Sequencing?[END_REF]. Ions within the fragmentation zone are fragmented. The abundance of the fragment ions is modulated with the same frequency as the radius of the precursors, i.e., f ICR -f min . The third pulse has a starting phase of 2πf min (2T 1 + t 1 + τ m ) and excites all ions to high radius before detection. Fragment ions start accruing phase during the third pulse [START_REF] Van Agthoven | Phase relationships in two-dimensional mass spectrometry[END_REF].

In the horizontal dimension, the phase φ p,h accrued by precursor ions that are not fully de-excited to the center of the ICR cell at the end of the second pulse, can be expressed with a quadratic phase dependence as:

φ p ,h (f 1 )=c 2 ' f 1 2 +c 1 ' f 1 +c 0 ' +2 π f 1 t 1 c 1 ' =c 1 +T 1 + τ m +T 2 +T 3 (2)
in which c 0′ , c 1 , and c 2′ are the parameters for the quadratic phase incurred during excitation, f 1 is the cyclotron frequency of the ion, t 1 is the encoding delay, T 1 is the duration of the second pulse, τ m is the fragmentation period, T 2 is the duration of the third pulse, and T 3 is the delay between excitation and detection [START_REF] Xian | Automated Broadband Phase Correction of Fourier Transform Ion Cyclotron Resonance Mass Spectra[END_REF][START_REF] Qi | Phase correction of Fourier transform ion cyclotron resonance mass spectra using MatLab[END_REF][START_REF] Van Agthoven | Phase relationships in two-dimensional mass spectrometry[END_REF].

The phase φ f,h accrued by fragment ions and unfragmented precursor ions that have been fully de-excited to the center of the ICR cell at the end of the second pulse, can be expressed with a quadratic phase dependence as:

φ f , h ( f 2 )=c 2 ' ' f 2 2 + c 1 '' f 2 +c 0 ' ' c 0 ' ' =c 0 +2 π f min (2 T 1 +t 1 +τ m ) (3)
in which c 0 , c 1′ ', and c 2′ ' are the parameters for the quadratic phase incurred during excitation, f 2 is the cyclotron frequency of the ion, t 1 is the encoding delay, T 1 is the duration of the first and the second pulse, τ m is the fragmentation period. As a rule, in phase-sensitive Fourier transform mass spectrometry, a shift of origin induces a zero order phase shift (c 0′ ' in Equation ( 3)), a delay induces a first order phase rotation (c 1′ ' in Equation ( 3)), and a frequency-dependent excitation induces a second order phase rotation (c 2′ ' in Equation ( 3)).

In the vertical dimension, equation 1 shows that the phase φ f,v for precursor ion signal can be expressed with a linear phase dependence as:

φ p ,v ( f 1 ) =2 π ( f 1 -f min ) ( t 1 -T 1 ) (4)
in which f 1 is the cyclotron frequency of the ion, t 1 is the encoding delay, f min is the lowest frequency in the pulse, and T 1 is the duration of the second pulse.

Multiple studies have shown that, in the vertical dimension, for fragmentation methods that have a maximum efficiency at the center of the ICR cell, the phase φ f,v of the fragment ion signal is shifted by π [START_REF] Marzullo | Advantages of Two-Dimensional Electron-Induced Dissociation and Infrared Multiphoton Dissociation Mass Spectrometry for the Analysis of Agrochemicals[END_REF][START_REF] Van Agthoven | Phase relationships in two-dimensional mass spectrometry[END_REF][START_REF] Paris | Facile Determination of Phosphorylation Sites in Peptides Using Two-Dimensional Mass Spectrometry[END_REF]. The phase of the fragment ion signal can therefore be expressed as:

φ f , v ( f 1 ) =2 π ( f 1 -f min ) ( t 1 -T 1 ) + π (5)
in which f 1 is the cyclotron frequency of the precursor ion, t 1 is the encoding delay, f min is the lowest frequency in the pulse, and T 1 is the duration of the second pulse.

Because signals have a t 1 amplitude modulation, a data set recorded in 2D mass spectrometry experiments requires processing using hypercomplex Fourier transformation [START_REF] Delsuc | Spectral Representation of 2D NMR Spectra by Hypercomplex Numbers[END_REF][START_REF] Schuyler | Formalism for hypercomplex multidimensional NMR employing partial-component subsampling[END_REF][START_REF] Ernst | Principles of Nuclear Magnetic Resonance in One and Two Dimensions[END_REF][START_REF] Le Bihan | Instantaneous frequency and amplitude of orthocomplex modulated signals based on quaternion Fourier transform[END_REF][START_REF] Sehgal | Theory for spiralling ions for 2D FT-ICR and comparison with precessing magnetization vectors in 2D NMR[END_REF]. The resulting spectrum is a 4-quadrant data set:

S ( f 1 , f 2 ) =RR ( f 1 , f 2 ) +iRI ( f 1 , f 2 ) + jIR ( f 1 , f 2 ) + kII ( f 1 , f 2 ) (6)
In which S is the signal in frequency domain, f 1 and f 2 are the frequencies, and i, j, and k are constants following the rules:

i 2 = j 2 =-1 ; k 2 =1 (7)
i . j= j . i=k [START_REF] Van Agthoven | Two-dimensional mass spectrometry in a linear ion trap, an in silico model[END_REF] i .k =k . i=-j [START_REF] Driver | Two-Dimensional Partial-Covariance Mass Spectrometry of Large Molecules Based on Fragment Correlations[END_REF] j . k =k . j=-i [START_REF] Van Agthoven | Towards analytically useful two-dimensional Fourier transform ion cyclotron resonance mass spectrometry[END_REF] With these rules, phase corrections can be applied independently to S(f 1 , f 2 ) along the f 1 and f 1 axes:

S (f 1 , f 2 ) e i φ 1 =S (f 1 , f 2 ) (cos φ 1 +isin φ 1 ) (11) S (f 1 , f 2 )e j φ 2 =S ( f 1 , f 2 ) (cosφ 2 + j sin φ 2 ) (12)
The phase corrections in the horizontal dimension and in the vertical dimensions can be applied independently from each other.

Experimental Methods

Sample Preparation

The sample used was an equimolar mixture of unmodified and modified (mono-, di-and trimethylated) model peptides (C-terminal GK-biotinylated histone H3 sequences, AnaSpec, Fremont, CA, USA) prepared with the same protocol as described in a previous article. [START_REF] Halper | Narrowband Modulation Two-Dimensional Mass Spectrometry and Label-Free Relative Quantification of Histone Peptides[END_REF] The sequence of the histone peptides was ATKAARKSAP ATGGVKKPHR YRPGGK biotin .

Instrument Parameters

Both the standard tandem mass spectrum and the two-dimensional mass spectrum were acquired on a 7 T Apex Ultra FT-ICR mass spectrometer (Bruker Daltonik GmbH, Bremen, Germany) with an electrospray ion source operated in positive mode and direct injection at a flow rate of 70 μL/hour. In both spectra, transients were acquired for 489.34 ms (512k datapoints) with an excitation pulse before detection with a frequency range of 535714.29-74728.13 Hz (corresponding to m/z 202.22-1450). The excitation pulse was made of 739 decrements of 624 Hz with 20 μs per frequency. The amplitude of the excitation was 106 V pp and the delay between excitation and detection was 3 ms.

For the standard tandem mass spectrum, ions were accumulated for 0.1 s in the first hexapole and 0.1 s in the second hexapole. In the quadrupole, ions were isolated at m/z 491 with an isolation window of m/z 30, and in the collision cell, they were activated with a 1.5 V collision energy before transfer to the Infinity ICR cell [START_REF] Caravatti | The infinity cell: A new trapped-ion cell with radiofrequency covered trapping electrodes for Fourier transform ion cyclotron resonance mass spectrometry[END_REF]. The fragmentation method was ECD, using a hollow cathode [START_REF] Tsybin | Combined infrared multiphoton dissociation and electron capture dissociation with a hollow electron beam in Fourier transform ion cyclotron resonance mass spectrometry[END_REF]. The heater was set at 1.3 A, the lens at 20 V, the bias at 1.7 V, and the irradiation lasted for 0.07 s.

The 2D mass spectrum was acquired in narrowband modulation mode [START_REF] Halper | Narrowband Modulation Two-Dimensional Mass Spectrometry and Label-Free Relative Quantification of Histone Peptides[END_REF]. The acquisition lasted 20 min. The ions were accumulated for 0.2 s in the first hexapole and 0.2 s in the second hexapole. In the quadrupole, ions were isolated at m/z 491 with an isolation window of m/z 10. In the collision cell, ions were activated with a 1.5 V collision energy. The ECD cathode was heated at 1.4 A, the lens set at 20 V, the bias at 2.0 V, and the ions were irradiated for 0.06 s. In the encoding sequence, the pulses had a frequency range of 535714.29-74728.13 Hz (corresponding to m/z 202.22-1450), with 739 decrements of 624 Hz with 1.0 μs per frequency. The amplitude of the encoding pulses was 106 V pp . The encoding delay incremented 1024 times by 50 μs, corresponding to a 10 kHz bandwidth, with 1 scan per increment. The signal was folded over 14 times, leading to a frequency range of 214659.79-224659.79 Hz (corresponding to m/z 482.177-504.602 for precursor ion modulation) [START_REF] Halper | Narrowband Modulation Two-Dimensional Mass Spectrometry and Label-Free Relative Quantification of Histone Peptides[END_REF]. The parameters of the excitation-detection sequence were identical to the parameters used for the standard tandem mass spectrum, with transients lasting 489.34 ms (512k datapoints).

Data Processing

The standard one-dimensional tandem mass spectrum and the 2D mass spectrum were processed and visualized using the Spectrometry Processing Innovative Kernel (SPIKE) software developed by the University of Strasbourg (France) and CASC4DE (Illkirch-Graffenstaden, France) in a 64-bit Python 3.7 programming language on an open-source platform distributed by the Python Software Foundation (Beaverton, OR, USA) [START_REF] Chiron | Spike a processing software dedicated to Fourier spectroscopies[END_REF]. Processed data files were saved using the HDF5 file format. No denoising was applied to either dataset [START_REF] Van Agthoven | Two-dimensional Fourier transform ion cyclotron resonance mass spectrometry: Reduction of scintillation noise using Cadzow data processing[END_REF][START_REF] Chiron | Efficient denoising algorithms for large experimental datasets and their applications in Fourier transform ion cyclotron resonance mass spectrometry[END_REF][START_REF] Bray | Nonuniform Sampling Acquisition of Two-Dimensional Fourier Transform Ion Cyclotron Resonance Mass Spectrometry for Increased Mass Resolution of Tandem Mass Spectrometry Precursor Ions[END_REF].

The transient for the one-dimensional tandem mass spectrum was apodised with a sine-bell window that was shifted, with a maximum of the bell at 15% of the transient, zero-filled twice, and Fourier transformed. The resulting spectrum was phase-corrected quadratically using a locally developed plugin to the SPIKE program (available at www.github.com/spike-project).

The program used for the data processing of the 2D mass spectrum in phasecorrected absorption mode is listed in the Appendix. Each transient in the 2D data set was apodised with a sine-bell window that was shifted, with a maximum of the bell at 15% of the transient, zero-filled twice and Fourier transformed. The modulation induced by the frequency generator was then removed by a complex rotation of the opposite frequency in the vertical dimension along the t 1 delay. The demodulation frequency was determined to be 74659.79 Hz. The exact demodulation frequency was offset by 68 Hz from the lowest frequency entered in the Apex Control software . The correct demodulation frequency was measured by processing the 2D mass spectrum with the nominal frequency entered in the software, which causes peak-splitting on the autocorrelation line, and measuring the frequency difference between the split peaks. Digital demodulation with the corrected frequency eliminated the peak-splitting [START_REF] Van Agthoven | Two-Dimensional ECD FT-ICR Mass Spectrometry of Peptides and Glycopeptides[END_REF]. The parameters determined for the one-dimensional tandem mass spectrum for quadratic phase correction were then applied to each row in the spectrum. The imaginary part of each row of the data set was then removed to reduce the computer memory burden. Then, each column of the data set was apodised with a slightly shifted sine-bell window (with a maximum of the bell at 15% of the transient), zero-filled twice and Fourier transformed. The theoretical linear phase correction determined above was applied on each column without modification. The imaginary part of the signal in each column was dropped and the dataset stored on disk.

The same 2D MS data set was processed in magnitude mode with identical apodisation and digital demodulation. However, due to the limitations of computer memory, the data-set was only zero-filled once in each dimension.

Results and Discussion

Figure 2 shows the result of quadratic phase correction with the SPIKE plugin for the one-dimensional ECD tandem mass spectrum of the four histone peptides. Although phase correction for FT-ICR mass spectra has been fully automated, this program requires user optimization [START_REF] Kilgour | An Algorithm for Automated Generation of Absorption Mode Spectra for FT-ICR MS[END_REF][START_REF] Kilgour | Producing absorption mode Fourier transform ion cyclotron resonance mass spectra with non-quadratic phase correction functions[END_REF][START_REF] Smith | Absorption mode FTICR mass spectrometry imaging[END_REF][START_REF] Kilgour | Improved optimization of the Fourier transform ion cyclotron resonance mass spectrometry phase correction function using a genetic algorithm[END_REF][START_REF] Kilgour | Absorption mode Fourier transform mass spectrometry with no baseline correction using a novel asymmetric apodization function[END_REF]. The phase correction method is based on initial phase correction of small regions of the spectrum with linear phase functions, in a similar way to the methods proposed by Xian et al. [START_REF] Xian | Automated Broadband Phase Correction of Fourier Transform Ion Cyclotron Resonance Mass Spectra[END_REF] and Qi et al. [START_REF] Qi | Phase correction of Fourier transform ion cyclotron resonance mass spectra using MatLab[END_REF]. The coefficients of the overall quadratic phase correction function are estimated by using the fact that these linear phase functions are tangents of the overall quadratic phase function. The user can then optimize the coefficients until the whole spectrum is properly phase-corrected.

Figure 2a shows the absorption mode mass spectrum before phase correction and Figure 2b shows the absorption mode mass spectrum after phase correction. Figure 2b shows that the signal for all the ions in the mass spectrum has been properly corrected. An insert shows a zoom-in on the isotopic distribution of the z 24 4+ fragment of the monomethylated histone peptide and compares the phase-corrected absorption mode mass spectrum to the magnitude mode mass spectrum. Both the resolving power and the signal-to-noise ratio are almost double for the phase corrected absorption mode spectrum than for the magnitude mode spectrum, as predicted by theory [START_REF] Xian | Automated Broadband Phase Correction of Fourier Transform Ion Cyclotron Resonance Mass Spectra[END_REF][START_REF] Qi | Phase correction of Fourier transform ion cyclotron resonance mass spectra using MatLab[END_REF]. The one-dimensional tandem mass spectrum was acquired with excitation pulse parameters and an excite-detect delay that was identical to the excitation pulse and the excite-detect delay of the 2D MS pulse sequence. Qi et al. have showed that the phase correction function is very stable with experimental conditions [START_REF] Qi | Variation of the Fourier Transform Mass Spectra Phase Function with Experimental Parameters[END_REF]. The phase correction function determined for the one-dimensional tandem mass spectrum can therefore also be used for the phase correction of the fragment ion signals in the 2D mass spectrum, which are, analytically, the most informative signals [START_REF] Van Agthoven | Phase relationships in two-dimensional mass spectrometry[END_REF].

The only difference between the excitation pulse for the tandem mass spectrum and the 2D mass spectrum is the initial phase of the pulse, which is determined by the continuous phase pulse generator. As is shown in Figure 1, the phase shift of the excitation pulse in the 2D MS experiment is 2πf min (2T 1 + t 1 + τ m ), in which f min is the minimum frequency in the pulse, T 1 is the duration of the first pulse, τ m the irradiation period, and t 1 is the incremental delay. The term 2πf min (2T 1 + τ m ) is constant throughout the experiment (and can be corrected with a constant phase offset in the phase correction function), but 2πf min t 1 is different for each transient.

Digital demodulation eliminates the phase contribution of 2πf min t 1 in each transient. In magnitude mode, digital demodulation is used to reduce the number of harmonics visible in the 2D mass spectrum [START_REF] Van Agthoven | Two-dimensional mass spectrometry: New perspectives for tandem mass spectrometry[END_REF][START_REF] Van Agthoven | Two-Dimensional ECD FT-ICR Mass Spectrometry of Peptides and Glycopeptides[END_REF]. Here, by applying digital demodulation before the quadratic phase correction in the horizontal dimension, we can eliminate the effect of the continuous phase pulse generator on the phase of fragment ion signals.

Table 1 gives an overview of the theoretically and empirically calculated coefficients for both phase correction functions used to calculate the absorption mode 2D mass spectrum (see Appendix A for the python code). The differences between the theoretical and calculated values shows that an operator is still needed to obtain the optimal phase correction function.

Table 1. Comparison between theoretically and empirically calculated coefficients for the phase correction functions in the vertical dimension (linear phase correction function) and in the horizontal dimension (quadratic phase correction function). The theoretical values are obtained from equation 28 of ref. [START_REF] Kilgour | An Algorithm for Automated Generation of Absorption Mode Spectra for FT-ICR MS[END_REF], ignoring image charge effect, and ref [START_REF] Van Agthoven | Phase relationships in two-dimensional mass spectrometry[END_REF]. Here, f min is the lowest frequency in the excitation pulse and ∆f the frequency range in the excitation pulse, T 1 is the duration of the first encoding pulse, T 2 the duration of the third excitation pulse and T 3 is the delay between excitation and detection, f N1 is the Nyquist frequency in the vertical precursor ion dimension and f N2 is the Nyquist frequency in the horizontal fragment ion dimension (see Figure 1).Zero order coefficients are expressed in degrees, first and second order coefficients are expressed in turns over the whole spectrum.

Coefficient Theory Theoretical Value Empirical Value 0 order (vertical)

f min T 1 88.2 59.9 1 st order (vertical) fN1T 1 7.44 7.8 0 order (horizontal) N/A N/A -9

1 st order (horizontal)

T 2 f min +T 3 ∆ f ∆ f f N 2 632 564 2 nd order (horizontal) -T 2 ∆ f f N 2 2 4607 4595
Figure 3. Extracted fragment ion scan from the 2D mass spectrum for the 1 × 13 C isotope of (M + 6H) 6+ ions of the monomethylated histone peptide (K7 1m) at m/z 489.453, zoom-in on fragment ion z 24 4+ , comparison between absorption mode and magnitude mode.

Figure 3 focuses on the extracted fragment ion scan of the 1 × 13 C isotope of (M+6H) 6+ ions of the monomethylated histone peptide (K7 1m). The fragment z 24 4+ contains 24 out of the 26 residues of the peptide. Consequently, if the precursor ion contains exactly one 13 C, then the probability of z 24 4+ retaining the 13 C is very high [START_REF] Van Agthoven | Two-dimensional mass spectrometry: New perspectives for tandem mass spectrometry[END_REF][START_REF] Halper | Narrowband Modulation Two-Dimensional Mass Spectrometry and Label-Free Relative Quantification of Histone Peptides[END_REF][START_REF] O'connor | High-resolution ion isolation with the ion cyclotron resonance capacitively coupled open cell[END_REF][START_REF] O'connor | Isotopic Assignment in Large-Molecule Mass Spectra by Fragmentation of a Selected Isotopic Peak[END_REF][START_REF] Smith | Ultrahigh Resolution Ion Isolation by Stored Waveform Inverse Fourier Transform 21 T Fourier Transform Ion Cyclotron Resonance Mass Spectrometry[END_REF]. In Figure 3, the relative intensity of the 1 × 13 C isotope of z 24 4+ is much higher than the relative intensity of the 12 C isotope. We can also see a peak for the 2 × 13 C isotope of z 24

4+

, which is caused by scintillation noise from the fragment of the 2 × 13 C isotope of the precursor ion. The isotopic distribution shown in Figure 3 is a consequence of the precursor ion signals being isotopically resolved in the vertical dimension. In the horizontal dimension, the phase-corrected absorption mode has doubled the resolving power and signal-to-noise ratio of the magnitude mode, as expected [START_REF] Xian | Automated Broadband Phase Correction of Fourier Transform Ion Cyclotron Resonance Mass Spectra[END_REF][START_REF] Qi | Phase correction of Fourier transform ion cyclotron resonance mass spectra using MatLab[END_REF]. Figure 4 shows the extracted precursor ion scan of the 12 C isotope of the c 6 fragment of the histone peptides. Since the modifications are all on the 7 th residue from the Nterminus, all histone peptides in the sample have identical c 6 fragments. This precursor ion scan is therefore an appropriate column of the data set to test the linear phase correction function. Equations ( 4) and ( 5) predict that the phase correction function has a slope of f N T 1 , in which f N is the Nyquist frequency in the vertical dimension (i.e., 10 kHz in the present experiment) and T 1 is the duration of the first pulse (i.e., 739 μs). Therefore, the theoretical first order phase correction should be of 7.39 turns over the whole spectrum. The resulting spectrum is well-phased, but we observed that 7.8 turns give slightly better results (see Appendix A). Figure 4 (left) shows that the phase correction function yields a well-corrected spectrum in the vertical dimension. A comparison with the extracted precursor ion scan from the magnitude mode spectrum (Figure 4 right) shows that the absorption mode spectrum has doubled the resolving power and signal-to-noise ratio compared to the magnitude mode spectrum, as is predicted [START_REF] Xian | Automated Broadband Phase Correction of Fourier Transform Ion Cyclotron Resonance Mass Spectra[END_REF][START_REF] Qi | Phase correction of Fourier transform ion cyclotron resonance mass spectra using MatLab[END_REF]. However, as this experiment was acquired in narrowband mode vertically, the frequency bandwidth is about twenty times narrower than in broadband 2D mass spectra, which may lead to a different phase correction [START_REF] Halper | Narrowband Modulation Two-Dimensional Mass Spectrometry and Label-Free Relative Quantification of Histone Peptides[END_REF].

The zoom-in on the isotopic distribution of the c 6 fragment of the dimethylated histone peptide (Figure 4 centre) shows that the peaks are baseline-resolved. The fullwidth at half-maximum (FWHM) of the peaks is approximately 35 mDa (R = 14000 at m/z 490), which is sufficient to isotopically resolve the fragments of biomolecules up to 14 kDa. Figure 5a shows a zoom-in on the 2D mass spectrum in phase-corrected absorption mode with only the positive values in the 2D mass spectrum. A similar 2D mass spectrum of the same sample was analyzed in detail in a previous article [START_REF] Halper | Narrowband Modulation Two-Dimensional Mass Spectrometry and Label-Free Relative Quantification of Histone Peptides[END_REF]. We notice that all isotopic distributions are isotopically resolved in both dimensions. An insert with a zoom-in on the isotopic distribution of the c 14 2+ fragment of the monomethylated histone peptide shows that clearly. Figure 5a shows mostly fragment ion peaks with very few artefacts. Figure 5b shows the negative values in the same region, which correspond to the 2ω harmonic of the ion signal in the vertical dimension [START_REF] Van Agthoven | Optimization of the discrete pulse sequence for two-dimensional FT-ICR mass spectrometry using infrared multiphoton dissociation[END_REF]. The modulation frequency of the unmodified and the monomethylated peptides is between 146.34 and 147.34 kHz.

Since the frequency range is 10 kHz, this leads to 14 foldovers for these signals. Because of the amplitude modulation of the signal in 2D MS and the real Fourier transform used during processing along the vertical axis, every time a signal is folded over in the 2D mass spectrum, the phase is shifted by 180° (or π radians). In this experiment, the signal has been folded over 14 times, and the phase is shifted by a multiple of 360°. The frequencies for the 2ω harmonic are between 292.68 and 294.68 kHz, so the signal has been folded over 29 times. Because the signal has been folded over an odd number of times, the 2ω harmonic of each dissociation line is a decreasing hyperbola, but over an m/z 3 range, it may appear linear [START_REF] Van Agthoven | Two-dimensional mass spectrometry: New perspectives for tandem mass spectrometry[END_REF]. The phase of the 2ω harmonic has been shifted by 29 times 180°, and the peaks corresponding to the 2ω harmonic have negative intensities, which simplifies data visualization and interpretation of the 2D mass spectrum.

The insert in Figure 5b shows the 2ω harmonic peaks for the isotopes for the z 24

4+

fragment of the unmodified histone peptide with frequency-to-mass conversion that has been adapted to the 2ω harmonics. An additional phase correction of -33° was applied to the extracted precursor ion scans to correctly phase the 2w harmonic peaks. As has been discussed in a previous article, the harmonic signals in a 2D mass spectrum can be used to increase the resolving power, just like they can in one-dimensional mass spectra [START_REF] Halper | Narrowband Modulation Two-Dimensional Mass Spectrometry and Label-Free Relative Quantification of Histone Peptides[END_REF][START_REF] Thomas | Petroleomic depth profiling of Staten Island salt marsh soil: 2ω detection FTICR MS offers a new solution for the analysis of environmental contaminants[END_REF][START_REF] Grosshans | Can Fourier transform mass spectral resolution be improved by detection at harmonic multiples of the fundamental ion cyclotron orbital frequency?[END_REF][START_REF] Vorobyev | Towards data acquisition throughput increase in Fourier transform mass spectrometry of proteins using double frequency measurements[END_REF][START_REF] Park | Parallel Detection of Fundamental and Sixth Harmonic Signals Using an ICR Cell with Dipole and Sixth Harmonic Detectors[END_REF]. Here, we almost double the resolving power from R = 14,000 to R = 27,000. Figure 6 shows the precursor ion scans for the 12 C isotopes of all four histone peptides on the autocorrelation line of the 2D mass spectrum. A comparison between equations ( 2) and (4) shows that, in the horizontal fragment ion dimension, the phase of precursor ion signals behaves differently from the phase of fragment ion signals, and that the phase correction function for precursor ion signals depends on t 1 . Furthermore, as shown in equations ( 2) and (3), in the vertical precursor ion dimension, the phase of precursor ion signals is shifted by -180° compared to the phase of fragment ion signals [START_REF] Van Agthoven | Phase relationships in two-dimensional mass spectrometry[END_REF]. Figure 6 confirms that, with phase correction functions that aim to correct the phases of fragment ion signals, precursor ion peaks on the autocorrelation line have negative intensities, and that their phase shifts with their cyclotron frequency (and their m/z ratio).

Conclusions

In this study, we have for the first time obtained a phase-corrected absorption mode two-dimensional mass spectrum by performing a quadratic phase correction in the horizontal fragment ion dimension and a linear phase correction in the precursor ion dimension. The resulting 2D mass spectrum shows the expected performance improvements in terms of resolving power and signal-to-noise ratio [START_REF] Xian | Automated Broadband Phase Correction of Fourier Transform Ion Cyclotron Resonance Mass Spectra[END_REF][START_REF] Qi | Phase correction of Fourier transform ion cyclotron resonance mass spectra using MatLab[END_REF][START_REF] Van Agthoven | Phase relationships in two-dimensional mass spectrometry[END_REF]. An added benefit of phase correction is the clear differentiation between "ground state" ion signals and folded-over harmonics, which, in our case, have negative intensities. This result can be expanded to broadband 2D MS, for harmonics that have been folded over. Careful experimental design can lead to good use of this effect.

Phase-corrected absorption mode narrowband 2D MS is shown here to enable isotopically resolved tandem mass spectrometry for biomolecules up to 14 kDa on a 7 T FT-ICR mass spectrometer. This performance opens the way for the identification, location, and relative quantification, of modifications such as citrullination or deamidation in peptides and proteins by two-dimensional mass spectrometry [START_REF] Halper | Narrowband Modulation Two-Dimensional Mass Spectrometry and Label-Free Relative Quantification of Histone Peptides[END_REF][START_REF] Lam | Does deamidation of islet amyloid polypeptide accelerate amyloid fibril formation?[END_REF][START_REF] Stensland | Targeted analysis of protein citrullination using chemical modification and tandem mass spectrometry[END_REF]. In future studies, phase correction for 2D MS needs to be expanded to broadband 2D mass spectra and automated for easy data processing.

Figure 1 .

 1 Figure 1. Pulse sequence for two-dimensional mass spectrometry and evolution of phases for (a) precursor ion radii, (b) fragment ion motion, (c) unfragmented precursor ion motion, (d) excitation pulses. (e) Duration of each segment of the pulse sequence. The factor α (resp. β) is related to the moment that the swept frequency of the chirp reaches the cyclotron frequency of the fragment (resp. unfragmented precursor).

Figure 2 .

 2 Figure 2. Tandem mass spectrum of the [M + 6H] 6+ histone peptides with ECD fragmentation in absorption mode (a) before phase correction and (b) after phase correction. Insert: zoom-in on the z 244+ fragment ion isotopic distribution of the monomethylated histone peptide (K7 1m), comparison between magnitude mode and phase-corrected absorption mode (FWHM: full-width at half-maximum).

Figure 4 .

 4 Figure 4. Precursor ion scan extracted from the absorption mode 2D mass spectrum for the 12 C isotope of the c 6 fragment of the histone peptides (m/z 616.38892), zoom-in on the c 6 fragment ion isotopic distribution in the 2D mass spectrum for the dimethylated histone peptide (K7 2m), comparison between the absorption mode and magnitude mode of the precursor ion scan of the c 6 fragment.

Figure 5 .

 5 Figure 5. Phase-corrected absorption mode 2D mass spectrum of histone peptides, zoomed in between m/z 600-700 horizontally and m/z 485-495 vertically, (a) positive values plotted (insert: zoom-in on the isotopic distribution of the c 14 2+ fragment of the methylated histone peptide) and (b) negative values plotted (insert: vertical precursor ion scans of the isotopic distribution of z 24 4+ of the unmodified histone peptide, with frequency-to-mass conversion for the 2ω harmonic).

Figure 6 .

 6 Figure 6. Extracted precursor ion scans for the 12 C isotopes of the four histone peptides from the phase-corrected absorption mode 2D mass spectrum.

Acknowledgments:

The authors thank Sarah Heel, Michael Palasser, and Christopher A. Wootton for helpful conversations. M.A.v.A. and K.B. thank Der Wissenschaftsfonds (Austrian Science Fund, FWF) for the Lise Meitner Fellowship project M 2757-B.

Funding: This research was funded by Der Wissenschaftsfonds (Austrian Science Fund, FWF) with the Lise Meitner Fellowship project M 2757-B. Open Access Funding by the Austrian Science Fund (FWF).

Data Availability Statement: In this section, please provide details regarding where data supporting reported results can be found, including links to publicly archived datasets analyzed or generated during the study. Please refer to suggested Data Availability Statements in section "MDPI Research Data Policies" at https://www.mdpi.com/ethics. You might choose to exclude this statement if the study did not report any data.

Author Contributions: All authors have contributed to the manuscript equally, read and agreed to the published version of the manuscript.

Conflicts of Interest:

The authors declare no conflict of interest.

Sample Availability