INTRODUCTION H + 2 REMPI ion source

The general aim of the "Trapped Ions" team at Laboratoire Kastler Brossel (LKB) is to perform high-resolution spectroscopy of the H + 2 molecular ion in order to test the Quantum Electrodynamics (QED) theory and improve the determination of fundamental constants, particularly the proton-to-electron mass ratio mp me . This constant is presently determined with a relative uncertainty of 9.5 × 10 -11 [START_REF]Codata mp/me page[END_REF][START_REF] Mohr | Codata recommended values of the fundamental physical constants: 2014[END_REF] from separate measurements of the proton [START_REF] Van Dyck | High precision penning trap mass spectroscopy and a new measurement of the protons atomic mass[END_REF] and electron [START_REF] Sturm | High-precision measurement of the atomic mass of the electron[END_REF] masses in Penning traps, respectively by mass spectrometry and by g-factor measurements. Very recently, the accuracy of the proton mass measurement was improved to 3.2 × 10 -11 [START_REF] Heiße | High-precision measurement of the proton's atomic mass[END_REF], resulting in a new value of mp me with 4.3 × 10 -11 accuracy. The new value of m p differs from the CODATA one by about 3 standard deviations.

The idea of using spectroscopy of rovibrational transitions in hydrogen molecular ions (H + 2 or its isotopes), combined with precise theoretical predictions, to determine nucleus-to-electron mass ratios was first formulated about 40 years ago [START_REF] Wing | Observation of the infrared spectrum of the hydrogen molecular ion HD +[END_REF]. On the experimental side a few one-photon ro-vibrational transitions in HD + have been measured with accuracies of 1-2 ppb in the last decade, limited by Doppler broadening [START_REF] Koelemeij | Vibrational spectroscopy of HD + with 2-ppb accuracy[END_REF][START_REF] Bressel | Manipulation of individual hyperfine states in cold trapped molecular ions and application to HD + frequency metrology[END_REF][START_REF] Biesheuvel | Probing QED and fundamental constants through laser spectroscopy of vibrational transitions in HD +[END_REF]. On the theory side, the precision has been improved through a systematic calculation of QED corrections and now reaches 7.6 × 10 -12 [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF].

The project of our team consists in probing a Doppler-free two-photon transition between the 1sσ g |ν = 0, L = 2 and 1sσ g |ν = 1, L = 2 ro-vibrational states of H + 2 at 9.17 µm. Doppler-free spectroscopy of hydrogen molecular ions is also in progress or planned in the groups of J. Koelemeij and W. Ubachs in Amsterdam [START_REF] Karr | Hydrogen molecular ions: new schemes for metrology and fundamental physics tests[END_REF] and S. Schiller in Düsseldorf [START_REF] Schiller | Quantum state preparation of homonuclear molecular ions enabled via a cold buffer gas: An ab initio study for the h 2 + and the d 2 + case[END_REF]. The expected accuracy of our measurement is a few 10 -11 , which would represent the most stringent test of QED in a molecular system and an important cross-check of the mp me value obtained from mass spectrometry and g-factor measurements. The precision may then be further improved to 10 -12 or better; it was shown [START_REF] Karr | Hydrogen molecular ions for improved determination of fundamental constants[END_REF] that by combining a few transition measurements of such accuracy in H + 2 and HD + it would be possible not only to improve the determination of the mp me and m d me mass ratios, but also to cross-check the values of the nuclear radii r p , r d and of the Rydberg constant, thus shedding light on the current discrepancy between different determinations of these constants [START_REF] Pohl | The size of the proton[END_REF][START_REF] Antognini | Proton structure from the measurement of 2s-2p transition frequencies of muonic hydrogen[END_REF][START_REF] Pohl | Laser spectroscopy of muonic deuterium[END_REF][START_REF] Beyer | The rydberg constant and proton size from atomic hydrogen[END_REF].

In order to perform our high-resolution spectroscopy experiment on H + 2 we must trap it. To do this we use a linear Paul trap which is the standard trapping technique in optical spectroscopy experiments on ion clouds, variants of hyperbolic Paul trap being preferred for single ions. One could also trap ions with a Penning trap [START_REF] Ghosh | Ion Traps[END_REF][START_REF] Major | Charged Particle Traps: Physics and Techniques of Charged Particle Field Confinement[END_REF] but the magnetic field in the Penning trap would mix the hyperfine level structure and change the frequency of the transition due to the Zeeman effect. In chapter 1 I introduce the theory of the linear Paul trap. Additionally to trapping the ions we need to cool them to avoid the shift and broadening of the transition we want to measure caused by the second order Doppler effect. Unfortunately, the level structure of H + 2 , as is the case for most molecules, is not suitable for laser cooling [START_REF] Champenois | Laser Cooling Techniques Applicable to Trapped Ions[END_REF] in particular because it is a homonuclear diatomic molecule, where all the ro-vibrational transitions are dipole forbidden so we plan to use sympathetic cooling. Sympathetic cooling is a technique to cool any ionic species close to the Doppler limit using an auxiliary ionic species which can be laser cooled. In this case the plan is to employ a Coulomb crystal of Be + laser cooled by a 313 nm laser to sympathetically cool H + 2 . The mechanism by which the cooling is transferred from one species to the other is Coulomb repulsion, the ions are maintained together by the trapping force and collide with each other which causes them to thermalise.

Due to the extremely long lifetime of the 1sσ g |ν = 1, L = 2 rovibrational states, the transition cannot be detected by fluorescence. In order to detect the 1sσ g |ν = 0, L = 2 → 1sσ g |ν = 1, L = 2 transition we will use a second laser at around 213 nm selectively photodissociate the ions excited to the ν = 1 level (more precisely, the photodissociation cross-section of the excited state is larger by almost 2 orders of magnitude than that of the ground state ν = 0). Maximal loss of H + 2 ions will therefore occur when the quantum cascade laser (QCL) at 9.17 µm is precisely tuned to the 1sσ g |ν = 0, L = 2 → 1sσ g |ν = 1, L = 2 transition. Fig. 0.2 illustrates the principle of the experiment.

In order to probe this transition we need to create H + 2 ions in the 1sσ g |ν = 0, L = 2 state. If one creates H + 2 ions by electron-impact ionisation of H 2 the ions are distributed among many states in a distribution which is known theoretically [START_REF] Busch | Photodissociation of H + 2 and D + 2 : Experiment[END_REF] and empirically [START_REF] Busch | Photodissociation of H + 2 and D + 2 : Experiment[END_REF][START_REF] Weijun | Vibrational population of H + 2 after electroionization of thermal H 2[END_REF]. Our team tried to carry out the experiment with such a non state-selective source but the signal was too weak [START_REF] Karr | Photodissociation of trapped H + 2 ions for REMPD spectroscopy[END_REF] because most of the ions are in states other than 1sσ g |ν = 0, L = 2 . Part of my work was to design and build a state-selective H + 2 ion source using a 3+1 Resonance Enhanced Photoionisation (REMPI) transition at 303 nm. In chapter 2 I discuss the design and show experimental results obtained with this ion source. The photodissociation of trapped H + 2 ions by a UV laser, which is another ingredient of the projected spectroscopy experiment, is studied in chapter 3.

Sympathetic cooling Simulations for GBAR

The GBAR (Gravitational Behaviour of Antihydrogen at Rest) international collaboration's goal is to measure the effect of Earth's gravitational field on H antimatter atoms. This experiment would represent the first test of the Weak Equivalence Principle with antimatter. This issue raises a lot of interest because any difference between matter and antimatter could, for example, shed light on the observed asymmetry between matter and antimatter in the universe. So far, the only available experimental data was obtained by analysing the trajectories of a few cold (T ∼ 1 K). H atoms released from a trap, resulting in very loose bounds on the Earth's gravitational acceleration for antimatter ḡ: -65 < ḡ g < 110 [START_REF] Collaboration | Description and first application of a new technique to measure the gravitational mass of antihydrogen[END_REF]. Three projects aiming for a more precise measurement of ḡ are underway at CERN: ALPHA-g [START_REF]Alpha collaboration website[END_REF], AEGIS [START_REF]Aegis collaboration website[END_REF] and GBAR [START_REF]Gbar collaboration website[END_REF].

The approach chosen by GBAR is to study the free-fall of H atoms. For this purpose the atoms have to be cooled to a very low temperature (T∼µK) in order to discriminate the ion's motion due to gravity from thermal agitation. Since the envisaged height of the free-fall is about 40 cm the goal is to have H with a root mean square velocity of less than 0.4 m s -1 . H could be Doppler cooled with a 121 nm laser but the Doppler and recoil limits (v recoil = 3.3 m s -1 ) would be too high for our purposes. The idea of Walz and Hänsch [START_REF] Walz | A proposal to measure antimatter gravity using ultracold antihydrogen atoms[END_REF] is to form an antimatter ion H+ (composed of an anti-proton and two positrons, the antimatter equivalent of H -), cool it by sympathetic cooling techniques and finally photodetach the excess positron very close to threshold to obtain a slow H atom. The H+ ions will be produced by the collaboration at the Antiproton Decelerator(AD) at CERN, using two successive charge exchange reactions between antiprotons and positronium atoms. This production of a three-body antimatter atom will, in itself, be a world first.

The antiprotons used to create H+ are provided to GBAR by CERN at high energies so a significant part of the collaboration's effort is devoted to reaching the 10-20µK goal temperature. One of the steps will consist in capturing and cooling the H+ ions from ≈1 eV to ≈1 mK using laser Doppler cooling. Unfortunately H+ cannot be laser cooled because it only has one bound state. It is therefore planned to sympathetically cool the H+ with a crystal of lasercooled Be + in a linear Paul trap, just like in our H spectroscopy experiment. This is what led the "Trapped Ions" team to become involved with the GBAR project as responsible of the design and operation of a "capture trap" where this cooling step will take place. In Fig. 0.3 we show a schematic of the GBAR experiment.

Achieving this sympathetic cooling is challenging because of the unfavourable mass ratio of 1 to 9 between Be + and H+ and because the H+ will be injected at relatively high energies. Indeed, sympathetic cooling of low energy ions such as those created directly in the trap is known to work experimentally but sympathetic cooling of ions created and loaded from outside the trap has only very recently been achieved [START_REF] Schmöger | Coulomb crystallization of highly charged ions[END_REF] on a quite different system : 40 Ar 13+ ions stopped by Be + . This is why I simulated this process during my thesis. The goals were to determine if the sympathetic cooling in the capture trap could work, and if so, under what conditions. Understanding sympathetic cooling better could also have applications for any other experiment which uses this technique, such as studies of cold highly charged ions [START_REF] Schmöger | Coulomb crystallization of highly charged ions[END_REF] or heavy molecular ions of biological interest [START_REF] Ostendorf | Sympathetic cooling of complex molecular ions to millikelvin temperatures[END_REF]. In order to carry out these simulations we needed the performance of Graphics Processing Units (GPU) and therefore I had to learn the CUDA programming language. In annex B I give an introduction to CUDA targeted towards a physicist having some general purpose programming knowledge. In chapter 4 I describe the simulation code and some tests performed in order to check its validity. Finally, in chapter 5 the results obtained from these simulations are presented.

THE PAUL TRAP

In this chapter we present the theory of the hyperbolic and linear Paul trap. The linear Paul trap will be used both in our H + 2 spectroscopy experiment and in GBAR. For this reason we chose this geometry for our simulations described in Chapters 4 and 5. In Fig. 1.1 one can see a photo of the trap constructed for the H + 2 experiment. Meanwhile, the hyperbolic Paul trap is also important to us because we used a trap of this type (which was the first trap operated in our team, see Fig. 1.2), to test the ion source discussed in Chapter 2 and perform the study of H + 2 photodissociation described in Chapter 3.

Hyperbolic Paul Trap 1.Theory

The hyperbolic Paul trap [START_REF] Ghosh | Ion Traps[END_REF][START_REF] Major | Charged Particle Traps: Physics and Techniques of Charged Particle Field Confinement[END_REF] was invented by Wolfgang Paul and earned him a Nobel prize in 1989. Trapping cannot be achieved with constant potentials due to the Maxwell equations which, in free space, forbid having a positive curvature of the potential in all three dimensions. Because of • V = 0 if the curvature is positive, and therefore trapping, in one or two dimensions then it must be negative in at least one other dimension. Paul found a way to trap ions using alternating voltages.

The hyperbolic Paul trap, as the name would suggest, has electrodes shaped like hyperboles represented in Fig. 1.3. If the cap electrodes at the top and bottom are at 0 potential and the ring electrode is at a potential V then the electrostatic potential Φ in the trap is

Φ = V d 2 x 2 + y 2 -2z 2 + 2z 2 0 (1.1)
with r 0 the shortest distance between the centre of the trap and the ring electrode, z 0 the shortest distance between the centre of the trap and the cap electrodes and d = r 2 0 + 2z 2 0 . The potential described in Eq. 1.1 shows equipotential surfaces with a hyperbolic shape. Because of the unicity of solutions of Laplace's equation's for given boundary conditions the electrodes must have this hyperbolic shape.

As we mentioned the Paul trap does not use a static voltage as in Eq. 1.1 but an alternating one, U (t) = U 0 + V 0 cos (Ωt), giving rise to the potential Φ = U 0 + V 0 cos (Ωt) 2d 2

x 2 + y 2 -2z 2 + 2z 2 0 .

(1.2)

Given Eq. 1.2 we can compute the electric field E in the trap

E = U 0 + V 0 cos (Ωt) d 2 (2zẑ -xx -yŷ) .
(1.3) This electric field can be separated into an RF and a DC part

E Ω = V 0 cos (Ωt) d 2 (2zẑ -xx -yŷ) (1.4) E DC = U 0 d 2 (2zẑ -xx -yŷ) (1.5)
From the expression of the electric field in Eq. 1.3 we can write the 3 independent equations of motion of an ion of charge q and mass m in the trap as

d 2 x dt 2 + q md 2 (U 0 + V 0 cos (Ωt)) x = 0 d 2 y dt 2 + q md 2 (U 0 + V 0 cos (Ωt)) y = 0 d 2 z dt 2 - 2q md 2 (U 0 + V 0 cos (Ωt)) z = 0. (1.6)
By introducing the stability parameters

a x = a y = - 4qU 0 md 2 Ω 2 , a z = 8qU 0 md 2 Ω 2 , q x = q y = 2qV 0 md 2 Ω 2 , q z = - 4qV 0 md 2 Ω 2 , (1.7)
and making a change of variable τ = 1 2 Ωt we can rewrite Eq. 1.6 as

d 2
x dτ 2 + (a x + 2q x cos (2τ )) x = 0 d 2 y dτ 2 + (a y + 2q y cos (2τ )) y = 0 d 2 z dτ 2 + (a z + 2q z cos (2τ )) z = 0.

(1.8)

Equations 1.8 are known as the Mathieu equations and they have oscillating or divergent solutions [START_REF] Landau | Mechanics[END_REF], depending on the values of the stability parameters a x ,a y ,a z ,q x ,q y and q z . It can be shown [START_REF] Major | Charged Particle Traps: Physics and Techniques of Charged Particle Field Confinement[END_REF] that there are so-called stability regions in the a x ,q x plane. In Fig. 1.4 on the left we show the first stability region of the linear Paul trap. Under the adiabatic approximation where a x , a y , a z , q x , q y , q z 1 it can be shown [START_REF] Major | Charged Particle Traps: Physics and Techniques of Charged Particle Field Confinement[END_REF] that the motion of the ion separates into a low amplitude forced motion at the frequency Ω called "micromotion" and a lower frequency oscillation called "secular motion" or "macromotion", which is a harmonic motion in the average trapping potential. For example, the x coordinate is of the form x (t) ≈ x 0 1 + q x 2 cos (Ωt) cos (ω x t), (1.9)

with ω x the secular frequency of the motion along x given by

ω x = 1 2 β x Ω (1.10)
with

β x = a x + q 2 x 2 .
(1.11) Here the red and blue areas correspond to the x and y directions respectively. By courtesy of J. Heinrich [START_REF] Heinrich | A Be + Ion Trap for H + 2 Spectroscopy[END_REF].

The macromotion is described by the so-called pseudopotential

U ef f = q 2 E 2 Ω (r) 4mΩ 2 = q 2 V 2 0 4mΩ 2 d 2 4z 2 + x 2 + y 2
(1.12)

and the depth of the trapping potential can be estimated by evaluating the pseudopotential at the position of the electrodes:

D z = qV 2 0 4mz 2 0 Ω 2 + U 0 2 (1.13) D r = qV 2 0 4mr 2 0 Ω 2 - U 0 2 (1.14)
Note that the actual potential in a real trap is only approximately given by Eq. 1.2. Indeed, to be perfect the Paul trap would need electrodes that extend to infinity. In practice they are finite and holes are drilled for ion loading, extraction, and optical access e.g. for laser cooling and fluorescence detection (see Fig. 1.2). Nevertheless, the above expressions are valid near the centre of the trap with r 0 and z 0 representing effective sizes, as opposed to real ones.

Orders of Magnitude

In this section we give orders of magnitudes of trapping parameters of a hyperbolic Paul trap. The parameters we give are those we typically use when operating our trap [START_REF] Bielsa | H + 2 two photons vibrational spectroscopy: 9.166 µm laser source development[END_REF]. In Tab. 1.1 we give stability parameters and trapping depths for H + 2 .

a z q z D z D r 0.002 0.191 3.02 eV 2.08 eV Tab. 1.1: Stability parameters a z and q z , trapping depths D z and D r for H + 2 in a hyperbolic Paul trap with parameters z 0 = 3 mm, r 0 = √ 2z 0 = 4.24 mm, Ω = 2π • 14.24 rad s -1 , U 0 = 0.76 V and V 0 = 143 V. 1.15r 0 because this geometry was found to minimize anharmonicities of the potential [START_REF] Reuben | Ion trajectories in exactly determined quadrupole fields[END_REF]. By courtesy of J. Heinrich [START_REF] Heinrich | A Be + Ion Trap for H + 2 Spectroscopy[END_REF].

Linear Paul Trap 1.2.1 Theory

The linear Paul Trap [START_REF] Prestage | New ion trap for frequency standard applications[END_REF][START_REF] Ghosh | Ion Traps[END_REF][START_REF] Major | Charged Particle Traps: Physics and Techniques of Charged Particle Field Confinement[END_REF] is a variation of the original hyperbolic Paul trap. The linear Paul trap can take many shapes but in Fig. 1.5 we show a representation of the geometry we use, with electrodes formed by four segmented cylindrical rods. The confinement in the x,y plane is ensured by applying a RF potential to a pair of diametrically opposed rods according to the same principle as for the 3D trap. DC potentials applied to the end electrodes provide confinement along the direction of the trap axis z. We will first describe the transverse motion of the ions and then the axial one.

While the original Paul trap has hyperbolic electrodes in order to generate an almost pure quadrupolar field, the linear Paul trap always approximates this original hyperbolic design (the rods are most often circular instead of hyperbolic), which gives rise to larger anharmonicities. In the trap shown in Fig. 1.5, the ratio between the electrode radius (r e = 4 mm) and the trap characteristic size (r 0 = 3.5 mm) is chosen in order to minimize the anharmonicities.

Near the axis, where x, y r 0 , the electrostatic potential Φ r is of the form

Φ r = U 0 + V 0 cos (Ωt) 2r 2 0 x 2 -y 2 + r 2 0 (1.15)
where r 0 is an effective size which slightly differs from the physical size R due to the deviation from a perfect quadrupolar potential. For the trap of Fig. 1.5 one finds r 0 = 0.9966R (see J. Heinrich's thesis [START_REF] Heinrich | A Be + Ion Trap for H + 2 Spectroscopy[END_REF]). The radial equations of motion, along x and y, are therefore:

d 2 x dt 2 + q mr 2 0 (U 0 + V 0 cos (Ωt)) x = 0 d 2 y dt 2 - q mr 2 0 (U 0 + V 0 cos (Ωt)) y = 0 (1.16)
By defining the stability parameters (note the sign on a y and q y compared to Eq. 1.7)

a x = -a y = 4qU 0 mr 2 0 Ω 2 q x = -q y = 2qV 0 mr 2 0 Ω 2 τ = Ωt 2 (1.17)
we find, just as we did with the hyperbolic Paul trap, the so-called Mathieu equations:

d 2 x dτ 2 + (a x + 2q x cos (2τ )) x = 0 d 2 y dτ 2 + (a y + 2q y cos (2τ )) y = 0 (1.18)
This leads us to the same conclusion that for certain values of a x and q x the trajectories are stable. In x(t) = x 0 cos (ω x,r t) 1 + q x 2 cos (Ωt)

y(t) = y 0 cos (ω y,r t) 1 + q y 2 cos (Ωt) (1.19)
with:

ω x,r = Ω 2 q 2 x 2 + a x and ω y,r = Ω 2 q 2 y 2 + a y . (1.20)
The motion of the trapped particle is therefore the superposition of a secular motion at frequencies ω x,r and ω y,r , and of a micromotion of frequency Ω. Note that as the radial confinement frequencies are higher for lower mass species, assuming equal charges (for a = 0, ω x,r ∝ q x ∝ m -1 ), when multiple species are trapped together this causes the lower mass species to be localized closer to the axis of the trap. The micromotion has a much higher frequency and a much lower amplitude than the macromotion. Averaging over the micromotion, the macromotion is governed by an effective potential

Φ pseudo,r = m 2q ω 2 x,r x 2 + ω 2 y,r y 2 (1.21)
Note that in Eq. 1.20 a y = -a x and q y = -q x therefore:

ω 2 y,r = ω 2 x,r - a x 2 Ω 2 . (1.22)
In simulations, sometimes the effective potential in Eq. 1.21 is used instead of the exact time-dependent potential in Eq. 1.15 and this is referred to as the pseudo potential approximation. The effect of U 0 is to break the symmetry between the x and y directions. In our simulations, we set a non-zero value of U 0 in order to prevent a perfect radial symmetry which would lead to a conservation of the angular momentum along the z axis, leading to a never ending rotation of the cloud of trapped ions. In experiments this is unnecessary since the symmetry is already broken by geometrical imperfections. Note that in cold ion experiments (including ours), DC potentials are often applied individually to each electrode in order to compensate for geometrical imperfections and/or patch potentials at electrode surfaces, in order to minimize ion micromotion [START_REF] Berkeland | Minimization of ion micromotion in a paul trap[END_REF]. This is for example important in high-resolution spectroscopy applications (such as ion clocks or our H + 2 spectroscopy experiment) in order to minimize the second-order Doppler shift.

Axial confinement, along the z axis, is provided by another set of electrodes, represented in Fig. 1.5, which create a field which is confining along z and deconfining along the radial directions x and y. This deconfinement must therefore be weaker than the confinement induced by the RF field in order not to compromise the radial confinement. The electrostatic potential Φ axial near the trap centre due to the endcap electrodes is of the form

Φ axial = K 2 z 2 - x 2 2 - y 2 2 (1.23)
where K depends on the trap geometry and applied voltages. The electric field component along the trap axis is therefore -Kz and an ion of charge q and mass m feels a force -qKz which leads to a harmonic motion of angular frequency ω z q m = qK m . K = 7969.18 N m -1 C -1 in our trap (see J. Heinrich's thesis [START_REF] Heinrich | A Be + Ion Trap for H + 2 Spectroscopy[END_REF]).

Please note that ω z is proportional to q m so that the product m q ω 2 z q m is independent of q and m. Meanwhile, the stability parameters a x , a y , q x , q y are proportional to q m . In experiments where ions of different masses/charges are trapped together (like those involving sympathetic cooling), care must be taken so that all the ions in the trap are in a stable zone of the stability diagram.

In our simulations of ion trapping, rather than K it is more convenient to specify the motional frequency ω z for a given ionic species of mass m and charge q. We thus express the axial electrostatic potential in terms of ω z :

Φ axial = mω 2 z q m 2q z 2 - x 2 2 - y 2 2 . (1.24)
The total trapping potential is the sum of the electrostatic potential created by the end caps given in Eq. 1.24 and the potential created by the central Species M Q a x =-a y q x =-q y ωx 2π ωy 2π ωz 2π

Be + 9 1 5.

= 0.1 V, V 0 = 200 V, Ω = 2π • 13 MHz, r 0 = 3.5 mm and ω z = 2π • 300kHz for Q = M = 1.
electrodes given in Eq. 1.15:

Φ = U 0 + V 0 cos (Ωt) 2r 2 0 x 2 -y 2 + mω 2 z q m 2q z 2 - x 2 2 - y 2 2 , (1.25)
from which we can deduce the electric field

E = U 0 + V 0 cos (Ωt) r 2 0 (yŷ -xx) - mω 2 z q m q zẑ - xx 2 - yŷ 2 . (1.26)
The total effective potential is given as the sum of Eq. 1.21 and Eq. 1.24

Φ pseudo = m 2q (ω 2 x,r - ω 2 z 2 )x 2 + (ω 2 y,r - ω 2 z 2 )y 2 + ω 2 z z 2 (1.27)
which means that the macromotion of ions follows a harmonic motion with radial frequencies

ω 2 x = ω 2 x,r - ω 2 z 2 ω 2 y = ω 2 y,r - ω 2 z 2 (1.28)

Orders of magnitude

An example of a trapping configuration is U 0 = 0.1 V, V 0 = 200 V, Ω = 13 MHz and ω z = 300 kHz for ions of mass M = 1 (in atomic mass units) and charge Q = 1 (in units of the elementary charge). In Tab. 1.2 we give the corresponding values of trapping parameters a x and q x and secular frequencies for different ion masses in such a linear Paul trap. We can also estimate a trap depth for this trap. Indeed, using the effective potential of Eq. 1.27 and the values of Tab. 1.2 we can, for 9 Be + , calculate the height of the radial potential barrier at the distance of the electrodes, r 0 = 3.5 mm away from the centre, to be 1.18 eV in the x direction and 1.07 eV in the y direction. For the z direction we calculate a potential barrier of 0.66 eV at a distance of 0.6 cm (at the edge of the central electrodes).

RF Heating

In radio-frequency traps there is an effect which has been empirically known for a long time [START_REF] Dehmelt | Radiofrequency spectroscopy of stored ions I: Storage[END_REF] called RF heating whereby energy from the RF micromotion is transferred to the macromotion. This effect is not yet fully understood theoretically, but it has been studied numerically in a number of papers [START_REF] Hudson | Sympathetic cooling of molecular ions with ultracold atoms[END_REF][START_REF] Cetina | Micromotion-induced limit to atom-ion sympathetic cooling in paul traps[END_REF][START_REF] Chen | Neutral gas sympathetic cooling of an ion in a paul trap[END_REF][START_REF] Schowalter | Blue-sky bifurcation of ion energies and the limits of neutral-gas sympathetic cooling of trapped ions[END_REF][START_REF] Moriwaki | Collision cooling of ions stored in quadrupole radio-frequency trap[END_REF][START_REF] Baba | Sympathetic cooling rate of gas-phase ions in a radio-frequency-quadrupole ion trap[END_REF][START_REF] Marciante | Dynamique d'ions en pièges radiofréquences[END_REF], and there have been recent theoretical insights [START_REF] Nam | Analytical mean-field scaling theory of radio-frequency heating in a paul trap[END_REF]. Numerical simulations have shown that RF heating by ion-ion collisions scales quadratically with temperature [START_REF] Pollock | Statistical mechanics of dense ionized matter. II. equilibrium properties and melting transition of the crystallized one-component plasma[END_REF]. It also scales with the amplitude of the micromotion which means that ions closer to the centre/axis of the trap suffer less from RF heating in the case of a hyperbolic/linear Paul trap respectively. The fact that the RF field-free region is a line in linear traps, as compared to a single point in hyperbolic traps, implies that a larger number of ions can be cooled in linear traps with an equivalent level of micromotion. This is one of the main arguments for using the linear geometry in experiments involving cooling of more than one ion.

Coulomb Crystal

When ions are trapped in a RF trap, there is an equilibrium between the trapping force which drives the ions towards the centre of the trap and their mutual Coulomb repulsion. If the ion cloud is sufficiently cold, such as may be the case from laser cooling, it takes a crystalline structure as can be seen numerically, in Fig. 4.2, and experimentally in Fig. 1.6. The photo in Fig. 1.6 is taken by a camera which can see the 313 nm fluorescence light from the Be + ions as they interact with the cooling laser. In Fig. 1.6 we can see individual ions. If the ions were in a fluid non-crystalline form we would not be able to see individual ions as they would move around during the exposure time of the camera causing blurring. When we speak of laser-cooled ion clouds we therefore prefer the term crystal to cloud. Coulomb crystals were first predicted theoretically [START_REF] Pollock | Statistical mechanics of dense ionized matter. II. equilibrium properties and melting transition of the crystallized one-component plasma[END_REF][START_REF] Dubin | Theory of structural phase transitions in a trapped Coulomb crystal[END_REF] and observed experimentally in hyperbolic traps [START_REF] Diedrich | Observation of a phase transition of stored laser-cooled ions[END_REF][START_REF] Wineland | Atomic-ion Coulomb clusters in an ion trap[END_REF]. They were later experimentally observed in linear Paul traps [START_REF] Drewsen | Large ion crystals in a linear paul trap[END_REF][START_REF] Hornekaer | Formation process of large ion Coulomb crystals in linear paul traps[END_REF] and studied in numerical simulations [START_REF] Schiffer | Temperature, ordering, and equilibrium with time-dependent confining forces[END_REF].

It was determined [START_REF] Schiffer | Temperature, ordering, and equilibrium with time-dependent confining forces[END_REF] that crystallisation occurs approximately when Γ = 170 or greater, where

Γ = q 2 4π 0 a sw k B T (1.29)
with a sw is the WignerSeitz radius (the inter-ion distance), T the temperature of the cloud/crystal and q the charge of the ions. In our simulations a sw is of the order of ∼ 10 µm, T ∼ 1 mK which means that γ ∼ 20 000, well past the crystallisation threshold of 170.

Theoretical shape for a single-component crystal

In this section we derive formulas for the dimensions of a crystal composed of one ionic species as a function of trapping parameters and ion numbers. We make the following assumptions:

• The DC voltage U 0 applied between the RF electrodes is sufficiently small for the asymmetry between the x and y directions to be negligible, i.e. ω x = ω y .

• The ions form a uniformly charged ellipsoid with a rotational symmetry along the z axis. This assumption is justified in Ref. [START_REF] Turner | Collective effects on equilibria of trapped charged plasmas[END_REF].

• The adiabatic approximation is valid.

• The crystal's macromotion is in equilibrium at zero temperature.

Under these assumptions, the trapping potential can be approximated by the time-independent pseudopotential

V ef f = 1 2 m q ω 2 x x 2 + y 2 + ω 2 z z 2 (1.30)
The electrostatic potential inside a uniformly charged oblate (a > b) spheroid is given by [START_REF] Turner | Collective effects on equilibria of trapped charged plasmas[END_REF] V

E (x, y, z) = 1 2 0 nqa 2 b 1 √ a 2 -b 2 arcsin a 2 -b 2 a 2 - 1 4 0 x 2 + y 2 nqa 2 b 1 (a 2 -b 2 ) 3 2 arcsin a 2 -b 2 a 2 - b (a 2 -b 2 ) a 2 + 1 2 0 z 2 nqa 2 b 1 (a 2 -b 2 ) 3 2 arcsin a 2 -b 2 a 2 - 2 (a 2 -b 2 ) b (1.31)
with n the density of ions, q the charge of the ions, a the radius of the ellipsoid in the (x,y) plane and b its radius along z. A slightly modified formula is available for the potential outside the spheroid [START_REF] Turner | Collective effects on equilibria of trapped charged plasmas[END_REF]. At equilibrium the trapping force from the potential in Eq. 1.30 must compensate the electromagnetic repulsion due to the potential in Eq. 1.31. This is possible only if V ef f + V E is constant in the crystal which leads us to

ω 2 z ω 2 x = -2 arcsin 1 -1 α 2 -α 1 -1 α 2 arcsin 1 -1 α 2 -1 α 1 -1 α 2 (α > 1) (1.32)
with α = a b the aspect ratio of the spheroid. Eq. 1.32 is only valid for oblate spheroids, when α > 1. For prolate spheroids with α < 1 the formula becomes

ω 2 z ω 2 x = -2 arcsinh 1 α 2 -1 -α 1 α 2 -1 arcsinh 1 α 2 -1 -1 α 1 α 2 -1 (α < 1) (1.33)
For known trapping parameters α can be obtained by inverting Eq. 1.32 or Eq. 1.33 using software like Mathematica. From Maxwell's equations and Eq. 1.30 we find

∆V ef f = ρ 0 = qn 0 = m q ω 2 x + ω 2 y + ω 2 z (1.34)
Eq. 1.34 is simple to obtain but is interesting in that it shows that the density of the crystal depends only on the stiffness of the effective trapping potential and the mass and charge of the ions. Once α has been calculated, the volume of the spheroid is given by 4 3 πa 2 b = 4 3 πα 2 b 3 but is also given by N ions n , which allows us to express

b = 3N ions 4πα 2 n 1 3 (1.35)
When a sympathetically cooled ion is at the edge e of the crystal with no kinetic energy left, its potential is given by V E (e) + V ef f (e). Since the sum of V E and V ef f must be constant inside the spheroid for the crystal to be at equilibrium, we can say that V E (e) + V ef f (e) = V E (O) with O the origin of the coordinate system and therefore using Eq. 1.31, Eq. 1.34 and Eq. 1.35 the potential of a sympathetically cooled ion oscillating along z with an amplitude of motion equal to the z axis b of the crystal is given by

q i V edge = 1 2 m ω 2 x + ω 2 y + ω 2 z a 2 b 1 √ a 2 -b 2 arcsin a 2 -b 2 a 2 = 1 2 m ω 2 x + ω 2 y + ω 2 z 3 4πα 2 n 2 3 N 2 3 ions α 2 √ α 2 -1 arcsin α 2 -1 α 2 ≈ 3 2 mω 2 z 3 4πn 2 3 N 2 3 ions (if α ≈ 1) = 3 2 mω 2 z 1 3 q 2 N ions 4π 0 2 3 (1.36)
for oblate spheroids, with q i the charge of the sympathetically cooled ion.

In our simulations, in order to mimic the GBAR experiment we send an ion with a certain energy along the z axis towards the centre of the spheroid. It then oscillates back and forth through the crystal with the amplitude of motion allowed by its initial energy. Eq. 1.36 gives an order of magnitude on what an ion crystal can easily cool as a function of how many ions compose it. When an ion's energy is much higher, it will spend most of its time outside the crystal. But if an ion's energy is comparable to the potential at the edge of the crystal then the ion will mostly oscillate inside the crystal and interact efficiently with the laser cooled ions. q i V edge therefore gives an idea of the "energy acceptance" of an ion crystal for sympathetic cooling. Fig. 1.7 shows V edge as a function of the number of ions for a typical set of trapping parameters. From this figure one might be led to think that sympathetic cooling of ions having high initial energies in the eV range, as required for the GBAR project, could be easily obtained using a sufficiently large ion crystal (with the parameters chosen for Fig. 1.7, V edge = 1 V with N ∼ 2 × 10 6 ). However, it should be recalled that this few-eV initial energy range corresponds to an energy dispersion. Since there is no way to measure individually the position and velocity of each incoming ion (which would allow adapting the injection setup accordingly) an ideal situation, where the ions reach the edge of the ion crystal with a very low residual kinetic energy, is not achievable. At best, one may expect that the incoming ions retain their initial energy dispersion when reaching the crystal's edge. Whatever the crystal's size, it is thus relevant to study the case E > q i V edge , where the ions spend most of their time outside the crystal. In this regime, the sympathetic cooling characteristic time could be long, and numerical simulations are required to assess its feasibility.

Conclusion

In this chapter, we have introduced the key concept of ion trapping which is crucial to both the experimental and numerical work presented in this thesis. We discussed how a RF Paul trap works and gave the relevant standard formulas used in the following chapters. We saw that cold trapped ions form Coulomb crystals such as shown in Fig. 1.6 and that a phenomenon called RF heating is present when ions are trapped together in a RF trap. We also derived the potential energy at the limit of the ion crystal because it corresponds to a minimum energy an ion crystal can cool and is therefore relevant to our study of sympathetic cooling.

ION SOURCE

Introduction

The goal being spectroscopy of sympathetically cooled H + 2 ions we need a way to create these ions and capture them in the a Be + crystal without compromising the vacuum in the vessel where the trap is located. One way is to create the ions "in situ" using ionization from a neutral molecular beam, this makes the capture straightforward but care must be taken not to compromise the vacuum. Another approach is to create the ions in a separate chamber and bring them into the trap with ion optics, the vacuum can be preserved with a differential pumping system, however, this approach requires a more complicated design.

Creating ions in the trap is simpler and is favoured in all experiments. Only recently have there been attempts to load the ions externally. One major reason to do this is when the ions cannot be created in the trap. Such is the case for highly charged ions [START_REF] Schmöger | Coulomb crystallization of highly charged ions[END_REF] that can only be produced by dedicated sources. Another example is H+ antimatter ions whose capture and sympathetic cooling is one of the steps of the GBAR project (as discussed in the Introduction) and will be studied by numerical simulation in Chapter 5. In our case, H + 2 ions can be easily produced in situ at room temperature by ionization of H 2 , and therefore we chose the simplest solution.

The easiest (and cheapest) way of producing H + 2 is electron-impact ionization of H 2 . This method was used in the previous version of our experimental setup, but has the important drawback that the resulting H + 2 ion population is distributed over many ro-vibrational states, with only about 1.4% in the state we want to probe (that is ν = 0, L = 2, in order to probe the 1sσ g |ν = 0, L = 2 → 1sσ g |ν = 1, L = 2 transition), severely limiting the achievable signal-to-noise ratio [START_REF] Karr | Photodissociation of trapped H + 2 ions for REMPD spectroscopy[END_REF]. This can be dramatically improved by using a state-selective creation process called REMPI, which we briefly describe in Section 2.2. For this method we need a 303 nm laser source described in Sec. 2.4 and a neutral H 2 beam, the design of which is described in Sec. 2.3 while in Sec. 2.3.1 we present vacuum theory required to calculate pressures in the ion source design.

REMPI State Selective Ionization Process

Resonance-enhanced multiphoton ionization (REMPI) allows the creation of ions in a state selective way [START_REF] Pratt | Excited-state molecular photoionization dynamics[END_REF]. The principle of REMPI, illustrated in Fig. 2.1, is to excite the molecule to a state the wave function of which resembles the wave function of the state we want the ion to be in with an additional loosely bound electron in a Rydberg state, and then ionizing it with an addi- tional photon (which may be either at the same wavelength as for the excitation step, or at a different one). The probability to reach a given ionic state is given in first approximation by the Franck Condon principle, i.e., is proportional to the squared modulus of the molecule/ion vibrational wave function overlap. Therefore there is a high probability for the vibrational quantum number of the ionic core to be preserved in this process because the Born-Oppenheimer curves are almost parallel (e.g. Fig. 1 of Ref. [START_REF] Anderson | Resonance-enhanced multiphoton ionization of molecular hydrogen via the E,F 1 Σ + g state: Photoelectron energy and angular distributions[END_REF]). This means that an ion population could in principle be created in any vibrational state simply by tuning the excitation laser to populate the corresponding vibrational state of the Rydberg electronic level. Since REMPI is a multiphoton process (with the first excitation step, towards an excited molecular state, often being multiphotonic itself), pulsed and strongly focused excitation lasers are usually required.

The REMPI ionization of small molecules has been studied both experimentally and theoretically in the 1970s and 1980s. In the specific case of H + 2 it was studied by Pratt et al. [START_REF] Ohalloran | Photoionization dynamics of H 2 C 1 Π u : Vibrational and rotational branching ratios[END_REF] and others [START_REF] Perreault | Angular and internal state distributions of H + 2 generated by (2 + 1) resonance enhanced multiphoton ionization of H 2 using time-of-flight mass spectrometry[END_REF][START_REF] Anderson | Resonance-enhanced multiphoton ionization of molecular hydrogen via the E,F 1 Σ + g state: Photoelectron energy and angular distributions[END_REF][START_REF] Urbain | Intense-laser-field ionization of molecular hydrogen in the tunneling regime and its effect on the vibrational excitation of H + 2[END_REF]. O'Halloran et al. showed that four-photon (3+1) REMPI creation of H + 2 in 1sσ g |ν = 0, L = 2 through the C 1 Π u intermediate state with ≈303 nm light followed well the Franck-Condon principle, i.e. ionizing transitions in which the vibrational state of the resonant intermediate state is preserved in the ion are most probable. More specifically, in that work a photoelectron spectra of the Q(2) line:

H 2 X 1 Σ + g (ν = 0, L = 2) → C 1 Π u (ν = 0, L = 2)
is reported (see Fig. 7 of Ref. [START_REF] Ohalloran | Photoionization dynamics of H 2 C 1 Π u : Vibrational and rotational branching ratios[END_REF]), showing that about 90% or the ions are formed in the ground vibrational state (ν = 0) (see also Ref. [START_REF] Pollard | Rotational relaxation in supersonic beams of hydrogen by high resolution photoelectron spectroscopy[END_REF]). The population of the L = 2 state was later estimated in Ref. [START_REF] Tran | Towards high precision spectroscopy of hydrogen molecular Ions H + 2 and HD + : theory and experiment[END_REF] (under some simplifying assumptions) to be larger than 97%. These results are very important for our experiment as they would represent an increase of the useful fraction of ions by a factor of about 60 with respect to the electron-impact ionization process.

Design of the Molecular Beam Apparatus

In this section we describe the design of the ion source which aims to create H + 2 ions inside the trap. This ion source must satisfy two main criteria:

• The ion creation rate must be high enough

• The vacuum in the trap chamber must not be compromised (< 10 -10 mbar)

The standard way of creating molecular ions in a trap is to send a beam of neutrals through the trap and ionise it close to the trap centre, which is referred to as in situ creation. But sending a beam through the setup will obviously degrade the vacuum. A common way of dealing with this problem is to have differential pumping stages separated by skimmers.

A skimmer, illustrated in Fig. 2.2, has a funnel-like shape but is placed with the tip facing the beam. The skimmer lets through the part of the beam that is most perfectly going through to the centre of the trap and rejects the rest of the beam. Beam molecules can be scattered by molecules reflected from the exterior surface of the skimmer. Making the exterior angle of the skimmer small minimizes this effect. Beam molecules can also be scattered by collisions with the gas that builds up inside the skimmer through collisions with the interior surface of the skimmer. To minimize this effect requires making the interior angle of the skimmer large. The curved shape of the skimmer is a compromise between these two conflicting requirements.

Another key technique to reduce the impact on pressure is to use a fast pulsed valve which sends pulses of neutral gas synchronised with the pulses of the REMPI laser instead of continuously sending gas. In Fig. 2.3 we show a photo of our ACPV2 pulsed valve from MassSpecpecD BV.

The design (see Fig. 2.4) is then to have a reservoir chamber filled with neutral H 2 connected to a buffer chamber via a pulsed valve. This buffer chamber is connected to the trap chamber via a skimmer and has its own pumping apparatus. This buffer chamber gets most of the molecular beam and is expected to have a higher pressure than the trap chamber.

Another design is to have two buffer chambers and therefore two skimmers in order to reduce further the pressure effect of the ion source on the trap chamber. Although the double skimmer design is often used (e.g see Ref. [START_REF] Tong | Collisional and radiative effects in the state-selective preparation of translationally cold molecular ions in ion traps[END_REF]), single skimmer designs are also possible [START_REF] Kokish | Note: High density pulsed molecular beam for cold ion chemistry[END_REF] and we chose the latter because the calculation we performed in Sec. 2.3.2 led us to believe that the pressure should be good enough for our needs.

Vacuum Theory Formulas Particle Flux Through a Hole

From a chamber at pressure P and temperature T the number of particles per second Ṅ of mass m flowing through a hole of area A is given by equation 2.1.

Ṅ = AP √ 2πmk B T (2.1)
This formula can be found in chapter II of N. Ramsey's "Molecular Beams" [START_REF] Ramsey | Molecular Beams[END_REF]. This formula is only valid for effusive beams because its derivation neglects particle collisions in the hole. A beam is effusive if the mean free path λ of the H 2 is greater than the diameter of the hole w, as expressed in Eq. 2.2.

w < λ (2.
2)

The expression of the mean free path is given in equation 2. 3 ([63] chapter II) with n the density and σ the collision cross section.

λ = 1 nσ √ 2 = k B T P σ √ 2 (2.3)
In this chapter we are interested exclusively in H 2 . The collision cross section between H 2 and itself we take to be 4.52 × 10 -20 m 2 , from [START_REF] Ramsey | Molecular Beams[END_REF] chapter XIII, table XIII.I. Equation 2.1 gives the total flux through a hole. This flux is not uniform in all directions, in fact it follows a cosine law in θ, where θ is the angle with respect to the normal of the surface containing the hole. The flux per unit of solid angle d Ṅ dΩ is given by equation 2.4.

d Ṅ dΩ = AP cos (θ) π 3 2 √ 2mk B T (2.4)
We have derived this equation and it can also be found in [START_REF] Ramsey | Molecular Beams[END_REF] chapter II. The consistency of Eq. 2.4 with Eq. 2.1 can also be checked by integrating over θ from 0 to π 2 . This equation is of course valid under the condition 2.2.

Mean velocity in a beam

The mean longitudinal velocity in a beam can be derived from the Boltzmann distribution. It can be found in chapter II of Ref. [START_REF] Ramsey | Molecular Beams[END_REF] and is given by

v = 3 4 2πk B T m (2.5)
For H 2 at room temperature the mean beam velocity is 2099 m s -1 . This value will be useful to us in upcoming calculations.

Equilibrium pressure

If we have a vacuum chamber with an incoming flow of particles Q and a pumping speed D specified by the manufacturer in m 3 s -1 what is the equilibrium pressure of the chamber? If the pump is removing D m 3 s -1 then it is removing particles at a rate of P D k B T s -1 , with P the pressure in the chamber and T its temperature. At equilibrium Q matches this rate and for a given temperature T the equilibrium pressure P eq is given by equation 2.6

P eq = Qk B T D (2.6)

Ion Source Equilibrium Pressures

In the following we describe the calculation of the equilibrium pressure in the trap chamber and the ion source buffer chamber.

• For this estimate, we have assumed that the reservoir connected to a H 2 bottle, is maintained at a pressure of 1 mbar. The reservoir leaks gas into the buffer chamber through the pulsed valve shown in Fig. 2.3, which has a nozzle diameter w = 150 µm. Using Eq. 2.3. we find that the mean free path in the reservoir is λ = 648 µm, so that the effusive beam condition of Eq. 2.2 is satisfied. As will be discussed later on, in our experiments we were led to use a much higher pressure, of the order of 1 bar. This has the important consequence that the condition of Eq. 2.2 is no longer satisfied, so that the molecular beam is in the supersonic regime rather than the effusive one. However, the order-of-magnitude estimates presented here remain a useful guide.

-In the effusive beam regime, for a pressure P = 1 mbar in the reservoir, using equation 2.1 we calculate a flow of 1.9 × 10 17 s -1 through a 150 µm diameter hole.

-The pulsed valve will send 50 µs pulses 20 times per second to match the firing rate of 20 Hz of the REMPI laser. The pulsed valve will therefore be open for 1 ms per second.

-This leads to a mean flow of 1.9 × 10 14 s -1 .

-This chamber being pumped by a Leybold Turbovac 450i turbo pump specified at 200 L s -1 we calculate an equilibrium pressure of 4 × 10 -8 mbar in the buffer chamber using equation 2.6.

• This buffer chamber then leaks gas into the main vessel through the 150 µm diameter skimmer shown in Fig. 2.2.

-Using equation 2.1 we calculate a flow of 7.5 × 10 9 s -1 through the skimmer into the main vessel.

-The main vessel is pumped by a Gamma vacuum 75 S ion pump and getter ribbon. Given the specified reachable vacuums for these pumps and the fact that we are combining them, we estimate the pumping rate of H 2 to at least 200 L s -1 .

-Given a 200 L s -1 pumping rate for H 2 in the vessel, using equation 2.6 we calculate an equilibrium pressure of 1.5 × 10 -12 mbar, well below the observed limit pressure of the chamber 1 × 10 -10 mbar due to outgassing, so the leak from the buffer chamber should not limit the main vessel pressure.

-In the same fashion we can compute the total flow of particles from the H 2 beam into the vessel through the skimmer. The 150 µm diameter skimmer at 15 cm from the pulsed valve has a solid angle of 0.14 × 10 -5 sr. Taking into account the duty cycle of the pulsed valve, this leads to a flow of 4.75 × 10 7 s -1 which is negligible compared to the 1.9 × 10 10 s -1 flow from the residual pressure in the buffer chamber. We thus expect that through these two stages of differential pumping we will have a low enough pressure. In fact the estimated value is comfortably lower than the pressure goal, which gives us confidence that the experiment has tolerance for experimental imperfections and possibly over-enthusiastic manufacturer specifications of H 2 pumping rates, H 2 being a hard gas to pump.

Ion Production Rate

The H 2 reservoir at 1 mbar, through the pulsed valve, generates an effusive beam in the buffer chamber which follows the cosine law in equation 2.4. As shown in Fig. 2.4, part of this beam goes through the skimmer and into the trap where the REMPI ionization laser is focused. We focus the laser at the centre of the trap with a waist of about w 0 = 10 µm. The corresponding Rayleigh length is 1 mm, which is greater than the shadow of the skimmer which is 458 µm by Thales' theorem, the pulsed valve being 18 cm from the skimmer, itself 37 cm from the trap. The ion production area for the beam thus represents a solid angle of 1.5 × 10 -8 sr (illustrated in Fig. 2.5). Over that small solid angle we can approximate θ ≈ 0 and equation 2.4 gives us a flux of 6 × 10 16 s -1 sr -1 in the direction θ = 0. Multiplying this flux by the 1.5 × 10 -8 sr solid angle we obtain a flow of 9.2 × 10 8 s -1 through the ion production area. We can understand that F = 9.2 × 10 8 s -1 flow as a linear density of particles multiplied by an average particle speed v. Therefore F v is the linear density in that flow and the average number of particles in the approximate volume of the beam focus is given by N

= F w 0 v (2.7)
From section 2.3.1 we know that the average velocity in the beam is 2099 m s -1 and we find an average particle number N of 4.4.

We therefore expect that if we saturate the REMPI process over the whole area w 2 0 we will create N ≈ 4.4 H + 2 ions per laser pulse, of which there are 20 per second, leading to an expected ion creation rate of 87 s -1 .

However, this figure needs to be corrected for the fact that our REMPI method only creates ions from neutral H 2 in the |ν = 0, L = 2 state. At room temperature all the molecules are in ν = 0 due to the high vibration frequency. But the fraction of molecules in L = 2 needs to be evaluated. Assuming a Boltzmann distribution at T = 300 K and a rotational energy constant of B = 1.2 × 10 -21 J [START_REF]NIST H 2 page[END_REF], we can calculate the fraction as

p L=2 = (2 • 2 + 1) e - B2(2+1) k B T Z(T) (2.8)
with k B the Boltzmann constant and Z the partition function given by

Z(T) = +∞ J=0 (2J + 1) 2 -(-1) J e - BJ(J+1) k B T (2.9)
where the 2 -(-1) J factor takes into account the hyperfine degeneracy of 1 for the even J levels and 3 for the odd J levels. We find p L=2 = 0.115, leading to a corrected ion creation rate of 10 s -1 , which leads us to believe we would be able to see an ion creation signal.

Conclusion

We have shown that the ion source design should in principle be able to create a sufficient ion flow to trap some ions within a few seconds of REMPI, without compromising the vacuum and we therefore believe that the ion source design is appropriate.

REMPI Laser

A key component of our state-selective H + 2 ion source is of course the laser that will excite the 3+1 photon transition at 303 nm. We use a frequency doubled Nd:YAG Quantel Brilliant laser, producing 5 ns pulses of 532 nm light, to pump a Sirah Cobra dye laser, producing 606 nm light which is frequency-doubled to 303 nm. The Nd:YAG laser has a repetition rate of 20 Hz and an energy of up to 160 mJ per pulse. The pulse energy is adjustable by changing the delay between the flash lamps and the Q-switch. The flash lamps are used to create the population inversion in the Nd:YAG rods and the Q-switch triggers the emission of a pulse. The dye laser produces about 16 mJ of 606 nm which is doubled to about 4 mJ of 303 nm. In Fig. 2.6 we show a curve of the dependence between the Q-switch flash lamp delay and the pulse energy at 303 nm.

The dye laser uses a mixture of Rhodamine-610 and Rhodamine-640 which we have empirically optimised to maximise the output energy at 606 nm. We used a total volume of about 1 litre of ethanol with a 0.0733 g L -1 concentration of Rhodamine-610 and a 0.0267 g L -1 concentration of Rhodamine-640. With respect to the mixture used previously in the group [START_REF] Tran | Towards high precision spectroscopy of hydrogen molecular Ions H + 2 and HD + : theory and experiment[END_REF] we slightly increased the proportion of Rhodamine 640 which we found empirically to shift the peak of the energy curve towards 606 nm. In Fig. 2.7 we show the pulse energy versus wavelength as we added more Rhodamine-640.

The spatial mode of the 303 nm is not very nice and rather elliptical. In order to make it more spherical we use a telescope made up of a divergent (f 1 = -62.6 mm) and a convergent lens (f 2 = 261 mm) placed 41.5 cm from the output and 17.5 cm apart. This setup and the resulting characteristics of the beam are described in more detail in Ref. [START_REF] Tran | Towards high precision spectroscopy of hydrogen molecular Ions H + 2 and HD + : theory and experiment[END_REF].

Experimental Realisation and Testing

We built the ion source described in Sec. 2.3 and tested it with the hyperbolic trap setup shown in Fig. 1.2. The overall setup is sketched in Fig. 2.4 while in Fig. 2.8 we show the mechanical design of the cross holding the pulsed valveskimmer combination while in Fig. 2.9 we show a photo of the experimental setup. The whole ion source described above is connected to the main chamber via a flexible bellow visible on Fig. 2.9 in order to be able to align the ion source with the centre of the trap. To obtain this alignement we use 3D printed pieces to spot the centre of flanges and, using a laser, we align the ion source with the center a flange of the ion trap chamber.

Our goal was then to produce ions to test the ion source. Created ions are measured destructively by extraction from the trap into a Microchannel plate (MCP) detector. We tried creating ions with a neutral H 2 pressure of 1 mbar behind the pulsed valve but could not detect ions. We were only able to detect a weak signal with the MCP at its maximal detection efficiency (in practice adjusted by a voltage) and a H 2 pressure of around 1 bar. We were able to identify this signal (shown in Fig. 2.10) as a clear signature of the REMPI process since we were able to observe different rotational lines at the same positions as in Fig. 1 of Ref. [START_REF] Pratt | Photoionization of excited molecular states. H 2 C 1 Π u[END_REF]. It can therefore be assumed that a high with a small shift due to imperfect calibration of our laser. There was 1 bar in the H 2 reservoir, the 313 nm was not at its best, at around 3.3 mJ, the gas pulses were 200 µs long and the creation time was 2 s.

proportion of the ions is produced in a single ro-vibrational state. The line of interest for our spectroscopy experiment is the Q(2) line at λ = 303.2 nm. We believe that the main reason why the expected ion creation rate is not achieved with a backing pressure of 1 mbar is most probably because our assumption of saturating the REMPI transition over an area w 2 0 is incorrect. An example curve of the H 2 REMPI ion signal vs pulse energy measured in an ion beam experiment (without trapping of the produced ion), provided by X. Urbain from UCL, Louvain-la-Neuve [66], shows that the saturation occurs for a pulse energy of the order of 20 mJ. If we assume that our focusing is identical, given our pulse energy of 3-4 mJ, the ion creation signal is weaker than the saturation value by a factor of about 30. Although this data was obtained on a different ro-vibrational line of the 3+1 REMPI process, with a different focusing of the laser, it gives a good indication that the REMPI yield could be 1 or 2 orders of magnitude below the saturation limit.

In addition, the focusing of our laser beam is probably not as tight as in the Louvain setup, since we observed that the spatial mode is far from being a TEM 00 . This may also be a critical factor: indeed the three-photon transition probability between the X and C states of H 2 scales like 1 w 6 0 (assuming that the C state excitation is not saturated and that the ionization process by the last photon is saturated), while the effective ion production volume scales like w 2 0 so that overall the ion production rate scales like 1

w 4 0 .
It would of course be desirable to increase the laser power P since the REMPI transition probability scales like P 3 . Since taking the data in Fig. 2.10, we replaced all the optics in the dye laser and recovered some laser power but we had many struggles with our pulsed valve which we later discovered was not bakeable, and were so far unable to repeat the experiment.

With 1 mbar we had calculated that the ion trap chamber pressure would be unaffected but using 1 bar caused the pressure in the ion trap chamber to spike up to ∼1 × 10 -9 mbar.

At 1 bar instead of 1 mbar the mean free path, calculated according to Eq. 2.3, is 0.65 µm instead of 650 µm. As a consequence, the Knudsen number [START_REF] Halwidl | Development of an Effusive Molecular Beam Apparatus[END_REF], which is the ratio of the mean free path to the opening of the valve (150 µm), becomes much smaller than 1, indicating that the beam is supersonic [START_REF] Miller | Free jet sources[END_REF] instead of effusive. Collisions of molecules in the beam can no longer be neglected and it is expected that the beam undergoes rotational cooling, reducing the proportion of H 2 in the L = 2 level, therefore reducing the proportion of neutral molecules that are susceptible to undergo the REMPI transition. However experimental data [START_REF] Pollard | Rotational relaxation in supersonic beams of hydrogen by high resolution photoelectron spectroscopy[END_REF] with a slightly smaller 70 µm valve diameter suggests that the decrease of L = 2 population is not dramatic: from ≈ 12% to ≈ 9%. Another point which deserves further study is the influence of the supersonic character of the beam on its angular distribution and therefore on its intensity close to the beam centre.

Conclusion and Perspectives

We have obtained experimental evidence for state-selective H + 2 ion production by the 3+1 REMPI process and subsequent trapping in our hyperbolic trap. However this required using a much higher reservoir pressure than originally planned. The large discrepancy between the observed ion production yield and our initial estimates may essentially be due to the characteristics of our 303 nm laser beam (pulse energy and focusing), which are probably far from saturating the REMPI transition.

One way to improve the situation greatly could be to sum 606 nm and 303 nm light in a non-linear crystal to obtain 202 nm light and perform 2+1 REMPI through the H 2 E,F 1 Σ + g intermediate state [START_REF] Perreault | Angular and internal state distributions of H + 2 generated by (2+1) resonance enhanced multiphoton ionization of H 2 using time-of-flight mass spectrometry[END_REF][START_REF] Anderson | Resonance-enhanced multiphoton ionization of molecular hydrogen via the E,F 1 Σ + g state: Photoelectron energy and angular distributions[END_REF]. Being a twophoton transition this transition should have a much greater cross section than the 3+1 REMPI we attempted. As a consequence, the REMPI transition could be saturated with a notably larger beam waist and therefore a much larger ion production volume.

If possible, it would be very helpful to upgrade the non-bakeable pulsed valve. Not being able to bake the buffer chamber contaminates the trap chamber with water and compromises vacuum. We were lucky to have been able to bake the experiment and perform the experiment shown in Fig. 2.10 as all further bakings led to a leaky valve and it took a while for us to realise the valve was not bakeable, as we had assumed any vacuum component would be bakeable to at least 100 • C.

Since taking the data of Fig. 2.10 we have connected the ion source to the trap chamber via a flexible bellow shown in Fig. 2.9. This has allowed us to align the ion source with the centre of the trap. This might improve the signal significantly but we have so far been unable to repeat the experiment.

H + 2 PHOTODISSOCIATION

For our H + 2 high-precision REMPD spectroscopy experiment, we plan to detect the 1sσ g |ν = 0, L = 2 → 1sσ g |ν = 1, L = 2 transition by dissociating the 1sσ g |ν = 1, L = 2 state with a 213 nm laser. A prerequisite for this is to test that the photodissociation of H + 2 works as expected. We will study this process with H + 2 ions created by electron-impact ionization of H 2 , because in this case the vibrational population distribution is well known. The motivation would then be to study the photodissociation of H + 2 ions produced by the REMPI ion source described in Chapter 2. This would allow us to characterize the vibrational distribution of the produced ions, and thus to check that a high percentage of ions are in the vibrational ground state v = 0 which is mostly unaffected by the 213 nm dissociation laser, as can be seen in Fig. 3.1.

213 nm Laser Source

The 213 nm laser is a commercial Xiton Impress 213 nm laser. It is a frequencyquintupled pulsed Nd:YVO 4 laser (λ = 1064 nm, similarly to a Nd:YAG laser) pumped by a temperature-controlled diode laser. Internally an AOM is used to Q-switch the 6.6 ns pulses at a frequency of a few kHz to a few tens of kHz. A maximum output mean power of about 170 mW is obtained for a repetition frequency of 12.5 kHz. The fifth harmonic is obtained by generating a second harmonic, summing the first harmonic to generate third harmonic and finally summing the third and second harmonic. All the required non-linear crystals are temperature controlled to an optimal value.

Model

Model for the photodissociation of an H + 2 ion cloud created by electron impact

When H + 2 is created from neutral H 2 using electron impact, it is created in a non-state-selective way in many different ro-vibrational states. The photodissociation of 1sσ g depends heavily on the vibrational state. In Tab. 3.1 we give the photodissociation cross-sections, calculated in our team [START_REF] Kilic | Le problème coulombien à trois corps : états liés, résonances et interaction avec la lumière vers la spectroscopie vibrationnelle dans H + 2[END_REF], for the different vibrational states as well as the vibrational distribution of H + 2 created by electron impact of neutral H 2 using the average values of [START_REF] Busch | Photodissociation of H + 2 and D + 2 : Experiment[END_REF] and [START_REF] Weijun | Vibrational population of H + 2 after electroionization of thermal H 2[END_REF].

Our spectroscopy experiment relies on the fact that the dissociation cross section is 72 times lower for ν = 0 than for ν = 1 as shown in Tab. 3.1 and Fig. 3.1, allowing us to detect the ν = 0 to ν = 1 transition by selectively dissociating ν = 1. However, looking at the cross sections of higher vibrational states (Tab. 3.1) one can see that ν = 4, ν = 7, ν = 10 and ν = 12 have lower cross sections than ν = 1. These states could therefore contribute background and noise to our REMPD signal. However, the population of these higher vibrational states is low in the case of electron-impact ionisation as can be seen in Tab. 3.1 and should be completely negligible with state-selective REMPI creation. Indeed we expect 90% of ions in ν = 0 and almost all other ions in ν = 1 [START_REF] Tran | Towards high precision spectroscopy of hydrogen molecular Ions H + 2 and HD + : theory and experiment[END_REF]. Note that the 313 nm Be + cooling laser has a much higher probability of dissociating these higher vibrational states. For example the ν = 7 state has a 5 × 10 -22 m -2 [START_REF] Dunn | Jila report[END_REF] dissociation cross section at 313 nm. The cooling laser has a power of about 1 mW compared to the dissociation laser's total power of 100 mW, so, assuming similar beam waists for the 313 nm and 213 nm lasers in the final experiment, we can estimate that the dissociation of ν = 7 by the 313 nm is an effect comparable to the dissociation by the 213 nm laser.

We will now calculate the decay of each vibrational state, having a photodissociation cross section σ ν . To do this we make the following assumptions:

• The hyperbolic trap contains a non-laser-cooled spherical cloud of H + 2 with a temperature T 1000K, and a Gaussian radial density profile of radius ∆r. This distribution was verified experimentally [START_REF] Knight | Laser scanning measurement of the density distribution of confined 6 Li + ions[END_REF][START_REF] Church | Radiative cooling of an electrodynamically contained proton gas[END_REF].

• The 213 nm laser beam is Gaussian and focused with a waist w 0 = 280 µm

(measured with a beam profiler) on the center of the trap.

• The characteristic thermalisation time of the ion cloud is much smaller than the photodissociation characteristic time, so that the radial density profile of the ion cloud can be considered constant during irradiation by the 213 nm laser. We will show that the time required to dissociate a significant fraction of the ions in the cloud is of the order of one second.

• The Rayleigh length of the Gaussian beam is large compared to ∆r which , and their photodissociation cross-sections at 213 nm. The populations are taken as the averages of the values given in [START_REF] Busch | Photodissociation of H + 2 and D + 2 : Experiment[END_REF] and [START_REF] Weijun | Vibrational population of H + 2 after electroionization of thermal H 2[END_REF]. In the following we neglect the population of ν > 12 vibrational states. The cross sections are from [START_REF] Kilic | Le problème coulombien à trois corps : états liés, résonances et interaction avec la lumière vers la spectroscopie vibrationnelle dans H + 2[END_REF].

allows us to neglect the variation of the beam waist with position in the cloud.

By the definition of cross section, and because the dissociation is a linear onephoton process, we can write the fraction of remaining ions in the vibrational state ν after one laser pulse as

f ν = e -T 0 σν φ(t )dt = e -σν φ photons T (3.1)
with φ photons the mean flux of photons during time T the repetition period of the laser. This can be rewritten as

f ν = e -σν Iλ hc T (3.2)
with I the mean irradiance of the beam. The irradiance in a Gaussian beam is of the form

I = I 0 e -2 ρ 2 w 2 0 (3.3)
with ρ the distance to the beam propagation axis and w 0 the waist of the laser. The total power of the beam is then

P total = +∞ 0 I(ρ) 2πρdρ = πI 0 w 2 0 2 (3.4)
and therefore

I 0 = 2P total πw 2 = 2E πw 2 0 T (3.5)
with E the energy of a pulse. By combining Eq. 3.5, Eq. 3.3 and Eq. 3.2 the remaining fraction for ions located at a distance ρ from the beam axis ions can be expressed as By combining Eq. 3.6 and Eq. 3.7 we can write the remaining fraction of ions in the cloud after one laser pulse as

f ν = e -2λσν E
f ν = 1 π 3 2 ∆r 3 +∞ 0 +∞ -∞ e -ρ 2 +z 2 ∆r 2 e -2λσν E πhcw 2 0 e -2 ( ρ w 0 ) 2 2πρdρdz
The integral over z is Gaussian and gives a factor √ π∆r so

f ν = 2 ∆r 2 +∞ 0 e -ρ 2 ∆r 2 e -2λσν E πhcw 2 0 e -2 ( ρ w 0 ) 2 ρdρ (3.8)
Finally, the total remaining fraction of ions after N pulses pulses, taking into account all the vibrational states, is given by

f diss = 12 ν=0 p ν f N pulses ν = 12 ν=0 p ν   2 ∆r 2 +∞ 0 e -ρ 2 ∆r 2 e -2λσν E πhcw 2 0 e -2 ( ρ w 0 ) 2 ρdρ   N pulses (3.9)

Trap losses

In the course of our experiments, we observed that the decay of the ion population under the effect of the photodissociation laser was in some cases significantly faster than predicted by Eq. 3.9, signalling the existence of an additional loss mechanism. In order to understand our data (Fig. 3.8), we revisited the model of Ref. [START_REF] Karr | Photodissociation of trapped H + 2 ions for REMPD spectroscopy[END_REF] describing ion losses in our hyperbolic trap ; firstly, in the absence of the photodissociation laser, and secondly in presence of laser irradiation at 213 nm.

In this section we aim to model the time evolution of the ion number in the hyperbolic Paul trap after they are created by electron impact. In most trapped ion experiments, trap losses are dominated by collisions or chemical reactions with background gas. In the case of H + 2 they mainly originate from the reaction H + 2 + H 2 → H + 3 + H [START_REF] Allmendinger | Experimental study of the ion-molecule reaction H + 2 + H 2 → H + 3 + H at low collision energies[END_REF]. H 2 is the most abundant species in a UHV environment and is constantly outgassed by the stainless-steel walls of the vacuum chamber or by the trap electrodes. In many cases, reactions with background gas lead to a simple exponential decay law, but here the situation is made more complicated by the fact that our electron gun causes a pressure "spike" (probably due to outgassing from the stainless steel parts induced by electron impact) which relaxes over time. This leads to a non-trivial time evolution of the trapped ion number. We make the following assumptions for our modelisation:

• The electron gun is on for a given time to create N 0 ions at time t 0 = 0.

• During this ion creation time the pressure increases to a value P 0 = P (t 0 ).

• The ions are trapped during a time t, extracted and counted.

• When the electron gun is off, the pressure decays exponentially to a limit pressure P ∞ with a time constant β. This β parameter is expected to be at least of the order of the vacuum chamber volume (in litre) divided by the pumping speed for H 2 (< 63 L s -1 specified by the manufacturer). Therefore β > 0.3 s, estimating the chamber volume at 20 L.

• The ion loss rate α in the trap depends linearly on the pressure as α = α 0 + γP . This implies that the ion decay rate does not depend on its internal state. The rate of the reaction with H 2 actually depends on the vibrational level ν of the H + 2 ion, but only at the level of a few percent per vibration quantum [START_REF] Eaker | A surface hopping quasiclassical trajectory study of the H + 2 +H 2 and (H 2 +D + 2 ) systems[END_REF][START_REF] Baer | A three-dimensional quantum mechanical study of the H 2 + H + 2 system: Calculation of reactive and charge transfer cross sections[END_REF][START_REF] Pollard | Stateselected reactive scattering. i. H + 2 +H 2 → H + 3 +H[END_REF], and this effect can be neglected here. The α 0 term aims to represent the effects of losses such as those due to instabilities from trap anharmonicities caused by, for example, the holes that are drilled through the electrodes.

The evolution of the ion number is determined by the following differential equation

dN (t) dt = -(α 0 + γP (t))N (t) (3.10)
while the pressure's time dependence can be written as

P (t) = P ∞ + (P 0 -P ∞ ) e -(t-t 0 ) β . (3.11)
with P 0 the pressure at time t 0 . Inserting Eq. 3.11 into Eq. 3.10 we obtain

dN N = -α 0 + γP ∞ + γ (P 0 -P ∞ ) e -(t-t 0 ) β dt (3.12)
which can be integrated to obtain

N (t -t 0 ) = N 0 e (-α 0 -γP∞)(t-t 0 )+γβ(P 0 -P∞) e - (t-t 0 ) β -1 . (3.13)
The remaining ion fraction F after time t -t 0 can be expressed as

F (t -t 0 , P 0 , P ∞ ) = N (t -t 0 ) N 0 = e (-α 0 -γP∞)(t-t 0 )+βγ(P 0 -P∞) e - (t-t 0 ) β -1
.

(3.14)

Trap Losses due to the 213 nm Laser

As explained above, there is experimental evidence of ion losses induced by the 213 nm laser, but due to some other process than photodissociation. These losses are probably due to the pressure increase we observe in the trap chamber during laser irradiation. This leads us to suppose that stray UV light hits some stainless steel parts, induces outgassing of H 2 and thus increases the ion loss rate due to the reaction of H + 2 with H 2 . To model this effect, we make the following assumptions:

• The laser repetition period τ = 80 µs is much smaller than the pressure evolution time constant, so that we can approximate the effect of the laser on the pressure as continuous.

• During laser irradiation, the pressure increases linearly with at a rate δP τ , where δP represents the pressure increase induced by one laser pulse.

• After being created the ions are trapped from t 0 to t 1 , exposed to the photodissociation laser from t 1 to t and then extracted for counting.

• The laser photodissociates ions, leaving a fraction f diss (t-t 1 ) in the trap, as calculated in Eq. 3.9 with N pulses = t-t 1 τ . Between t 0 and t 1 , the ion number evolves following Eq. 3.13 and according to Eq. 3.11 the pressure at time t

1 is P (t 1 ) = P ∞ + (P 0 -P ∞ ) e -(t 1 -t 0 ) β .
Between t 1 and the extraction time t the ion number evolves under the effects of photodissociation losses and laser-induced trap losses. Those two effects are multiplicative. Since we assume that the trap loss rate does not depend on ν, we can factorise this term for all ν states.

Pressure evolution

For t > t 1 , the pressure in the chamber evolves under the effect of pumping and outgassing following the differential equation

dP (t) dt = - (P (t) -P ∞ ) β + δP τ (3.15)
which leads to

P (t) = P L + (P 1 -P L ) e -(t-t 1 ) β (3.16)
with P 1 the pressure at time t 1 and P L = P ∞ + βδP τ the pressure limit reached when the dissociation laser is shined continuously. The pressure evolution with laser, given by Eq. 3.16, is similar to the result without laser of Eq. 3.11 with P ∞ replaced by P L and P 0 replaced by P 1 . Consequently, for t > t 1 the evolution of the ion number under the effect of trap losses (not including photodissociation losses) is given by Eq. 3.14 with the same substitutions, i.e.

N (t) = F (t -t 1 , P 1 , P L ) N (t 1 ).

Ion number evolution

In practice, we first create ions, and then trap them for a given time t trap while shining the laser during t diss (between t = t trap -t diss and t = t trap ). With these notations, we have t 1 -t 0 = t trap -t diss and t -t 1 = t diss . In order to measure ion losses due to photodissociation, we perform two such sequences, one with the UV laser on, and immediately afterwards, one with the laser off. We measure the remaining ion numbers, N U V ON and N U V OF F , and calculate their ratio R = N U V ON N U V OF F . In the absence of laser-induced trap losses, trap losses would be identical in both sequences and cancel out, so that R is simply equal to the fraction of unphotodissociated ion f diss , calculated in Eq. 3.9.

With the effect of laser-induced trap losses included, we have

N U V ON = N 0 f diss (t) × F (t trap -t diss , P 0 , P ∞ ) × F (t diss , P 1 , P L ) N U V OF F = N 0 × F (t trap , P 0 , P ∞ ) (3.17) so that R = f diss (t) F (t trap -t diss , P 0 , P ∞ ) F (t diss , P 1 , P L ) F (t trap , P 0 , P ∞ ) (3.18)
Using the multiplicative property

F (t trap , P 0 , P ∞ ) = F (t trap -t diss , P 0 , P ∞ ) F (t diss , P 1 , P ∞ ) (3.19)
this simplifies to

R (t diss ) = f diss (t) F (t diss , P 1 , P L ) F (t diss , P 1 , P ∞ ) . (3.20)
Finally, using Eq. 3.14 and the expressions of P 1 and P L , one obtains after some algebra:

R (t diss ) = f diss (t) e -γ(P L -P∞) t diss -β 1-e - t diss β (3.21)
In this equation the first factor represents dissociation losses, and the second one corresponds to trap losses induced by the UV laser.

Experimental Testing of Pressure and Trap Losses

In this section we experimentally test the different models of trap losses (both without and with the UV laser) presented in Sec. 3.2.2 and Sec. 3.2.3.

Pressure evolution without laser

By creating ions and monitoring the pressure after creation, one can test Eq. 3.11 and determine the limit pressure P 0 right after ion creation, the pressure P ∞ and the relaxation time constant β. As can be seen in Fig. 3.3, the situation is slightly more complicated since an exponential decay does not yield a very good fit.

A double exponential decay law works much better as can be seen in Fig. 3.3 with two very different time constants 2.43 s and 23 s. A possible interpretation is the following. The ion pump has a pumping speed of 63 L s -1 (for N 2 at 1 × 10 -6 mbar, less at 1 × 10 -9 mbar and even less for H 2 ) and the vessel volume is a few tens of litres. The resulting pumping time constant is lower bounded by 0.33 s because the gauge is far from the pump and therefore sees a reduced pumping speed. This could correspond to the 2.43 s time constant. During e-gun emission, trap parts are heated up and outgas. When e-gun emission is turned off, heating stops, temperature, outgassing and pressure decrease with a long (thermal) time constant which is consistent with the 23 s we observe.

One can also fit the beginning of the pressure decay with a single exponential, as shown in Fig. 3.4 by truncating the data at 10 s. Fig. 3.4 shows that short term decay gives β = 2.5 ± 0.1 s which is compatible with the 2.43 s ± 0.03 s of the double exponential fit. Because the photodissociation measurements were performed with a total trapping time of 3 s (1 s of creation and 2 s of trapping), we can indeed model the pressure variation as in Eq. 3.3, with a single exponential decay.

Pressure Evolution with Laser

By shining the laser continuously and monitoring the pressure, we can test Eq. 3.16. Experimental results are plotted in Fig. 3.5. Again, Eq. 3.16 is not fully satisfactory for a global fit (see left part) but yields a good fit of the short-term data (truncated at t = 10 s, see right part). In the latter case a time constant of 5.2 s is found. This is significantly longer than the pumping time constant β ≈ 2.4 s we found in Sec. 3.3.1, indicating that UV-induced outgassing is not instantaneous, maybe resulting from heating of metal pieces or ablation.

Repeating the experiment for different 213 nm laser powers we can determine the limit pressure in the trap chamber as a function of laser power. In Fig. 3.6 we show the results along with a linear fit, confirming the linear dependence between P L -P ∞ and laser power deduced from Eq. 3.15. Indeed from Eq. 3.15 we deduced P L -P ∞ = β δP τ and it isn't unreasonable to assume that δP τ might be proportional to laser power: P L -P ∞ = ηP total .

Trap Losses without Laser

In order to study trap losses, we create ions until time t 0 , trap them until time t when they are extracted and counted. According to Eq. 3.13 the ion number is expected to decay with a law of the form depending on 4 parameters A = N 0 e -βγ(P 0 -P∞) , B = α 0 +γP 0 , C = βγ (P 0 -P ∞ ) and β. In Fig. 3.7, we show the H + 2 ion signal versus trapping time after creating ions for 2 s with 1.05 A in the e-gun filament. Every point of data is averaged over 4 measurements. The data is fitted with Eq. 3.22 (fitting parameters given in caption) and the value of β we find is consistent with that obtained from the pressure decay in Sec. 3.3.1. The fit follows the experimental data well which confirms Eq. 3.13.

N (t -t 0 ) = Ae -B(t-t 0 )+Ce -(t-t 0 ) β (3.22)

Electron-Gun H + 2 Lifetime with Laser

Having checked our trap loss model in Sec. 3.3, we are ready to interpret our photodissociation experimental results. The fraction of non-photodissociated ions is measured in the following way:

• In a first sequence, the UV laser is off. Ions are created, trapped during a time t trap = 2 s, extracted and counted.

• The second sequence is identical, except for the fact that the UV laser is on during a time t diss at the end of the trapping period.

The ratio of measured ion signals F = N U V ON N U V OF F is plotted in Fig. 3.8 as a function of the interaction time t diss . We have shown experimentally in Fig. 3.6 that P L scales linearly with the laser power P total therefore according to Eq. 3.21, defining κ = γη, we can fit our data with a law of the form with which we fit the data in Fig. 3.8. To do this fit we used a bootstrap procedure of fitting with a variable β, κ and ∆r with constant measured P total , followed by a fit with the obtained values of β, κ and ∆r with variable P total , followed by a fit with the obtained P total and a variable β, κ and ∆r... for 10 iterations. This procedure converged but the β we found was very small (negligible compared to the interaction time) so we repeated the procedure with β = 0 s to obtain the fits of Fig. 3.8. The P total found by this procedure are similar to the ones we measured as shown in Tab. 3.2.

F (t diss ) = f diss (t) e -κP total t diss -β 1-e - t diss β (3.23)

Conclusion

We have shown in Fig. 3.8 that we were able to understand the dissociation rate of H + 2 created non-state-selectively by electron impact. The same type of work will need to be done on ions created by our REMPI ion source in order to double-check that the ions are indeed produced in a state-selective manner in ν = 0. We can also conclude that the principle of our REMPD measurement is sound. Indeed, we need to be able to dissociate ν = 1 without affecting ν = 0 to detect the ν = 0 → ν = 1 transition. With ∆r = 1.13 mm, w 0 = 280 µm, E = 7.5 µJ (we measured 67.5 mW at 9 kHz), Eq. 3.8 predicts that after 2 s of dissociation we would have dissociated 2% of ν = 0 and 74% of ν = 1, giving a ∼74% detection efficiency for H + 2 REMPD spectroscopy.

GPU NUMERICAL SIMULATIONS OF SYMPATHETIC COOLING

Introduction

The GBAR experiment aims to study gravity on H+ ions produced at CERN. These ions are produced with a high energy spread and must therefore be cooled down in order to study the comparatively weak effect of gravity. Unfortunately H+ cannot be laser cooled and the GBAR collaboration plans to cool down the H+ ions via sympathetic cooling by a laser-cooled crystal of Be + ions.

While the H+ ion is sympathetically cooled by the Be + it will be exposed to the 313 nm cooling laser which can cause photodetachement of the H+ . The photodetachement cross section of H -in the UV region has been calculated [START_REF] Broad | One-and two-electron photoejection from H -: A multichannel J-matrix calculation[END_REF] and measured [START_REF] Smith | Photodetachment cross section of the negative hydrogen ion[END_REF] allowing us to extrapolate a photodetachement cross section at 313 nm σ ≈ 2 × 10 -17 cm 2 . The cooling laser beam intensity is typically I ∼ 1 mW mm -2 , corresponding to a photon flux Φ = Iλ hc = 1.57 × 10 17 cm -2 s -1 . This leads to a photodissociation rate σΦ = 3.14 s -1 which means that the sympathetic cooling of the H+ ion must take place in less than 1 s. When ions are created outside the trap and guided in with ion optics they are loaded in with a certain amount of kinetic energy. One can understand that if this energy is too big, the ion will go through the crystal of ions that was meant to sympathetically cool it so fast that it will barely interact and exchange energy with it. In fact a cancer treatment method, ion beam therapy [START_REF] Schlaff | Bringing the heavy: carbon ion therapy in the radiobiological and clinical context[END_REF], is based on this fact that ions will not exchange much energy with matter when their speed is high, but as the beam slows down, it interacts more and more and abruptly gives all its energy to a very localized part of the matter.

The exchange of energy between ions and macroscopic matter is well understood in plasma physics [START_REF] Miyamoto | Fundamentals of plasma physics and controlled fusion[END_REF] but the sympathetic capture and cooling of externally produced ions in a mesoscopic crystal of a few thousand laser cooled ions is not. In fact it has only recently been achieved experimentally [START_REF] Schmöger | Coulomb crystallization of highly charged ions[END_REF]. Part of my work was to numerically simulate this system to better understand it and determine optimal parameters for the sympathetic cooling of H+ in the GBAR experiment. In this chapter we describe the physical system we solve, our implementation, the variable time step scheme we developed and the tests we performed to ascertain the validity of our simulations, which we performed on graphics processing units (GPU) for performance reasons. In chapter 5 we present results from our simulations.

Notation Meaning q

Charge in Coulombs m Mass in kilograms Q Charge in atomic units M Mass in atomic units ω z q m

The axial confinement frequency felt by an ion ωz

The axial confinement frequency felt by an ion with Q M = 1

Tab. 4.1: Table of Notations.

Numerical Model

In this section we detail the exact equations we use in our numerical model of the physical system. We solve Newton's equations for each ion.

Trapping Force

In the parameters of a simulation we can choose between the exact timedependent potential of an ideal linear trap, the pseudopotential approximation, or the trapping potential of a "real" trap based on Simion potential files.

Exact Trapping Force

In our parameter files (example given in appendix A) we define Ω the angular frequency of the trapping potential, V 0 the amplitude of the radiofrequency field, U 0 the dissymmetrising tension, r 0 the distance between the centre of the trap and the RF rods and ωz the angular frequency of the trapping potential along the z axis for a charge to mass ratio of 1 in atomic units. Those parameters are linked or identical to physical trapping parameters discussed in chapter 1. The force on an ion at time t in position (x, y, z) with charge q and mass m is given by Eq. 4.1 -4.3:

F x = m Q M ωz 2 2 -q U 0 + V 0 cos (Ωt) r 2 0 x (4.1)
F y = m Q M ωz 2 2 + q U 0 + V 0 cos (Ωt) r 2 0 y (4.
2)

F z = -m Q M ωz 2 z (4.3)
with Q and M the ion's charge and mass respectively, in atomic units.

Pseudopotential approximation

Under the pseudopotential approximation the trapping force no longer depends on time and the force on an ion becomes

F x = - 1 2 m Ω 2 2 q 2 x 2 + a x - Q M ωz 2 x (4.4) F y = - 1 2 m Ω 2 2 q 2 x 2 -a x - Q M ωz 2 y (4.5) F z = -m Q M ωz 2 z (4.6)
with q x and a x the stability parameters of the ion, given in Eq. 1.17.

Simion Potential Files

Our code can use Simion potential files generated from a model of the trap electrodes. This is useful to simulate trapping while including anharmonicites due to the geometry of the electrodes. We can also use this to precisely simulate the injection of the ions which is planned to be done by temporarily grounding one of the end cap electrodes to lower the potential barrier. Simion generates one potential file for each electrode, corresponding to a 3D grid of potential values assuming that the electrode under consideration is set to a potential of 10 000 V, while the others are grounded. We can then use this to deduce a 3D electrostatic potential grid for any electrode voltage configuration. The value of the field at any point is then calculated by interpolating the values of nearby grid points. We use an order 2 interpolation of the potential, requiring the use of the 27 nearest grid points, which we differentiate to obtain a first order interpolation of the electric field.

This functionality has been implemented in the code and tested, but unfortunately due to time constraints it has not been used for the simulations presented here.

Coulomb Force

We calculate the Coulomb force on each ion as the sum of the Coulomb forces upon it due to each other ion. The force on ion 1 with charge Q 1 due to ion 2 with charge Q 2 in units of the elementary charge is given by

F 12 = - q 2 e Q 1 Q 2 4π 0 d 2 r12 (4.7)
with d the distance between them, q e the elementary charge and r12 the unit vector going from ion 1 to 2.

Interaction with the Cooling Laser

In this section we describe the interaction between the Be + ions and the cooling laser at 313 nm. We use a Monte Carlo approach which reproduces the Doppler limit while being close to the real physical process, by opposition to fictitious cooling forces. We model the Be + ions as two level systems with a ground state, labelled as 1, and an excited state, labelled as 2, and our goal is to find expressions for the probability per timestep δt of absorption B 12 δt, stimulated emission B 21 δt and spontaneous emission A 21 δt. The laser has a Rabi frequency Ω R and a detuning δω. In most cases, it is parallel to the trap axis and passes Denoting Ω R the Rabi frequency of the laser and δω the detuning of the laser with respect to the 313 nm cooling transition, k the wave vector of the photons and v the velocity of the ion, we can express the absorption probability [START_REF] Cagnac | Physique atomique: introduction à la physique quantique et structure de l'édifice atomique[END_REF]

B 12 δt = Ω 2 R Γ + 4(δω-k•v) 2 Γ δt (4.9)
which can be rewritten as

B 12 δt = 1 2 Γs 1 + 4 δω-k•v Γ 2 δt (4.10)
by introducing the dimensionless saturation parameter s = 2

Ω 2 R Γ 2 .
The formula of [START_REF] Cagnac | Physique atomique: introduction à la physique quantique et structure de l'édifice atomique[END_REF] is valid for a plane wavefront, in our case we need to take into account the Gaussian nature of our beam. In the notations of Fig. 4.1 we can write

B 12 δt = 1 2 Γs 1 + 4 δω-kv Γ 2 w 0 w (CH) 2 e -2 MH 2 w(CH) 2 δt (4.11)
where s is now the saturation parameter at the center of the beam, at the beam waist (point C in Fig. 4.1). We have to express MH and CH in the coordinates of the trap. To do this we first have to change our coordinate system from having its origin in O to having its origin in C. We write (x , y , z ) the coordinates of M in this new coordinate system. Calling φ and θ the azimuthal and polar angles of the orientation of the laser we find

CH 2 = ((x cos φ + y sin φ) sin θ + z cos θ) 2 .
(4.12)

From the Pythagorean theorem we can write

MH 2 = CM 2 -CH 2 = x 2 + y 2 + z 2 -((x cos φ + y sin φ) sin θ + z cos θ) 2 .
(4.13) To express w (CH), Gaussian optics [START_REF] Kogelnik | Laser beams and resonators[END_REF] tells us that

w 0 w (CH) 2 = 1 1 + CH z 0 2 , (4.14)
with the Rayleigh length z 0 = πw 2 0 λ where λ is the wavelength of the cooling laser.

The expression of the absorption probability becomes

B 12 δt = 1 2 Γs 1 + 4 δω-kv Γ 2 1 1 + CH z 0 2 e -2 MH 2 w 2 0 1 1+ ( CH z 0 ) 2 δt (4.15)
We can then express the stimulated emission probability [START_REF] Cagnac | Physique atomique: introduction à la physique quantique et structure de l'édifice atomique[END_REF] as

B 21 δt = g 1 g 2 B 12 δt (4.16)
with g 1 the degeneracy of the ground state and g 2 the degeneracy of the excited state with g 1 = 1 and g 2 = 1 in our two level model of the ion. Using Eq. 4.8, Eq. 4.15 and Eq. 4.16 at every timestep, if an ion is in the ground (resp. excited) state we can compute the probability of photon absorption (resp. spontaneous and simulated emission). We draw random numbers between 0 and 1 to decide if a transition happens and, if it does, we apply the appropriate velocity kicks due to the momentum of the absorbed or emitted photon hk. The recoil velocity due to absorbing or emitting a photon is given by

v recoil = h λm lc = 1.3 m s -1 (4.17)
with m lc the mass of the laser-cooled beryllium. The recoil velocity in Eq. 4.17 is applied to the ion in the direction of the laser if a photon is absorbed, in the opposite direction if a photon is emitted due to stimulated emission and in a random direction if a photon is emitted due to spontaneous emission. We use the inverse sampling [START_REF] Sillitoe | Numerical Simulations Of Ion Cloud Dynamics[END_REF] method to draw random unit vectors uniformly on a sphere.

Integration Algorithm

We chose to simulate the physical model described in section 4.2 using the Velocity Verlet algorithm [START_REF] Verlet | Computer "experiments" on classical fluids. I. thermodynamical properties of lennard-jones molecules[END_REF] defined by Velocity Verlet has the advantage of calling the force only once, which is important to us due to the O(N 2 ) scaling of the Coulomb interaction. It is also symplectic (phase space surface is conserved) and time reversible. The Velocity Verlet algorithm is also valid when the timestep varies which is important to our variable timestep method presented in Sec. 4.6. However symplecticity and time reversal are lost when using a variable timestep with Velocity Verlet and different methods, requiring the resolution of implicit equations, would have to be used in order to recover it [START_REF] Huang | The adaptive verlet method[END_REF][START_REF] Hairer | Variable time step integration with symplectic methods[END_REF].

x (t + dt) = x (t) + v (t) dt + 1 2 a (t) dt 2 a (t + dt) calculated from x (t + dt) v (t + dt) = v (t) + 1 2 (a (t) + a (t + dt)) dt

Formation of Crystals

In order to simulate the dynamics of sympathetic cooling we first have to form the laser-cooled ion crystal that will do the cooling. In this section I describe the procedure we have developed to do this. We first initialise ions randomly in a cylinder as shown in the left part of Fig. 4.2, then we simulate with the trapping force, the Coulomb interaction and a fictitious friction force in order to reach equilibrium more quickly. After a short amount of time the ions arrange themselves in a crystalline structure as shown in the right part of Fig. 4.2. At this point, the ions are unrealistically cold, being well below the Doppler limit, because of the excessive friction force and the absence of a fluctuating force corresponding to spontaneous emission. In addition the ions are slightly displaced with respect to a more realistic treatment of cooling because the radiation pressure force is missing. So the last step of the crystal preparation is a step of interaction with the cooling laser in order to bring the crystal to equilibrium. In Fig. 4.3 we show the transition between fictitious friction force and laser cooling. After the simulation shown in Fig. 4.3 the position,velocity and acceleration of every particle is saved to a file so we can use this crystal as a starting point for any number of sympathetic cooling simulations.

We might for example from the same crystal try to cool some H+ ions , q x = 0.0223 for the Be + and ωz = 250 kHz. Characteristic lengths are given along with linear density n. By courtesy of J. Heinrich [START_REF] Heinrich | A Be + Ion Trap for H + 2 Spectroscopy[END_REF].

starting at different initial energies, while at a given energy we will run the simulation many times changing only slightly the initial position of the H+ ion in order to determine an average cooling time. Indeed, experience has taught us the cooling time fluctuates a lot between simulations, although they start from almost identical initial conditions. These large fluctuations in cooling times are probably due to the fact the cooling is very dependent on the occurence of rare collisions with a low impact parameter (close to head-on collisions). It is necessary to repeat each simulation a large number of times in order to get a reliable estimate of the cooling time.

In Fig. 4.4 we show a two-component crystal of Be + and H + 2 which is an example of what might be used for our H + 2 spectroscopy experiment.

Example Curves and Definitions

In this section, we show some example curves produced by our simulations. In most of our simulations we place one sympathetically cooled ion at rest close to the axis of the linear trap, with a specified angular tolerance seen from the centre of the trap. The ion's initial position is determined by its total initial potential energy (trapping potential + Coulomb interaction potential with the Be + crystal) which is specified in the parameter file along with a tolerance (typically 1%). This is implemented in a very simple way by randomly generating initial positions until the desired potential energy and angle is reached within the tolerances. The ion then goes back and forth through a laser-cooled crystal and will hopefully eventually be captured and cooled by the crystal.

By "capture" we refer to the moment where the sympathetically cooled ion loses enough energy such that the amplitude of its motion doesn't allow it to exit the crystal any more.

We can then look at many different quantities such as fluorescence photons by the Be + , but what interests us the most is usually the trajectory of the sympathetically cooled ion and its kinetic energy.

When we look at kinetic energy, we average over the RF period in order to eliminate micromotion, and we sometimes further average over a time that is much longer than the macromotion period in order to smooth the curve over macromotion variations. Indeed, the sympathetically cooled ion's kinetic energy goes down by orders of magnitudes when it reaches the maximal amplitude of its harmonic oscillator type motion and this type of variation does not interest us very much because the kinetic energy was merely converted to potential energy. We give an example of these types of curves in Fig. 4.5.

We can see in Fig. 4.5 that the H+ oscillates in and out of the crystal (reaching up to 700 m s -1 along the z axis, compared to ∼ 4 m s -1 of micromotion for a Be + radially 0.1 mm away from the trap axis from Eq. 1.19) without much radial amplitude in its oscillation at first, but as it gets pushed around by the Be + some of its axial (z) energy is converted to radial (x,y) energy. Over time its oscillation amplitude decreases and this is confirmed by the decreasing energy. Eventually, after about 7 ms the oscillation amplitude is low enough so that the H+ ion no longer gets out of the crystal and it takes a position at one end of the crystal, which has an axial extension of about 1 mm. The simulation is stopped because the program sees that the ion's energy has gone below the target energy specified in the parameter file (0.5 K in this case), the final stages of the cooling are not simulated and the program moves on to the next simulation with a different seed of the random number generator resulting in a slightly different initial position, as allowed by the angular and energy tolerances.

The sympathetically cooled ion tends to go to the opposite end of the crystal with respect to the laser beam direction because of the radiation pressure felt by the Be + and not by the sympathetically cooled ion.

In many cases, in order to study the practical feasibility of sympathetic cooling in the GBAR experiment, we are mostly interested in the dynamics of the sympathetic cooling process, i.e. by one number: the time it takes for the sympathetically cooled ion to go from its initial energy, oscillating back and forth through the laser-cooled crystal, to being captured inside the crystal. We refer to this as the capture time. In order to have reproducible numbers we have, from looking at many curves, decided to define our capture time as the time it takes for the sympathetically cooled ion to go below 2.5 K in averaged macromotion energy. Note that the 2.5 K limit for the capture time is different to the 0.5 K limit on stopping the simulation. After this, the ion will of course be further cooled closer to the Doppler limit, but these final stages of cooling, as well as the final temperature reached by the ion, are much less critical in the H+ is sent into a crystal of 1919 Be + . On top we show the trajectory of the H+ and on the bottom we show the averaged (right) and unaveraged (left) macromotion kinetic energy of Be + in purple and the H+ in green. The energy given is per ion which explains the small bumps in the purple curve (1919 ions) while the green one has bigger ones (1 ion). a x = 4.7 × 10 -6 , q x = 0.052 for the H+ , ωz = 65 kHz and other parameters as in Fig. 4.3. LC stands for the laser-cooled species, the beryllium in our case. SC stands for the sympathetically cooled species, H+ in this case.

context of the GBAR project. Mainly because once the sympathetically cooled ion is cooled to ∼ 2.5 K, it no longer leaves the cloud and this configuration is known to work experimentally and numerically. But also because for initially hot ions this phase is much shorter, and therefore negligible compared to the initial capture; and finally, because an exact knowledge of the final temperature is not crucial for the GBAR experiment. Indeed the sympathetically cooled ions will be transferred to another linear Paul trap for further sympathetic cooling, by a single Be + ion, to the motional ground state before performing the free-fall experiment.

Timestep Criteria

In any numerical simulation aiming to solve differential equations, the choice of the integration timestep is critical in order to find a trade off between simulation time (long timesteps) and precision (short timesteps). For instance, in the case of an harmonic oscillator, the Verlet algorithm is known to be stable for Ωδt < 1 and unstable otherwise. In our case, several criteria set limits on the timestep we can use. The one coming from the Coulomb interaction leads to a variable timestep scheme which we published [START_REF] Sillitoe | H+ Sympathetic Cooling Simulations with a Variable Time Step[END_REF].

Radio Frequency Trapping

The timestep should be short enough to properly describe the radiofrequency trapping potential of the linear Paul trap. One can naturally formulate the criterion that

δt 2π Ω ⇒ δt < β2π Ω , β 1 (4.19)
Because of this timestep requirement to represent the RF potential correctly cold ion physics simulations are often performed under the pseudo-potential approximation, allowing the use of longer timesteps at the cost of not being able to describe effects caused by the RF field like micro-motion and RF heating. This can be correct to study the equilibrium properties of ion crystals but cannot simulate dynamical behaviours. In Sec. 4.6.5 we discuss our choice for β.

Coulomb Collision

To understand sympathetic cooling our numerical simulations need to correctly describe the Coulomb collision of ions. One can understand from Fig. 4.6 that too long of a timestep will completely miss the collision of two ions as they go through each other. We therefore need a timestep criteria that can correctly describe the collision of two ions so that, when integrated into our simulations, it allows us to describe the N-body Coulomb interactions of our ion cloud.

Criterion Based on Minimum approach distance

To correctly describe a Coulomb collision one idea is to compare the distance δd, travelled by the ions during δt, to their minimum approach distance d min and impose δd d min (4.20)

Let us consider the collision of the two ions in the centre of mass reference frame. Let us call v their relative speed and a their relative acceleration. If we call δd the variation of distance between the two ions during one timestep, for our Velocity Verlet algorithm we have

δd ≤ vδt + 1 2 aδt 2 (4.21)
We can therefore write criterion 4.20 as

vδt + 1 2 aδt 2 = αd min (4.22)
with α << 1. Solving for δt we get

δt = -v + √ v 2 + 2aαd min a (4.23)
To evaluate d min , we assume a head-on collision and we note that using energy conservation

1 2 µv 2 + q 1 q 2 4π 0 d = q 1 q 2 4π 0 d min ⇒ d min = d 1 + 2π 0 µv 2 d q 1 q 2 (4.24)
with q 1 and q 2 the charges of the particles, d the distance between them and µ their reduced mass. In Fig. 4.7 we show that when simulating the collision of a projectile ion on a stationary target ion of equal mass there is a critical timestep below which the collision is well described which depends on the energy of the projectile ion. For too long a timestep almost no energy is exchanged as the ions pass through each other. But below a critical timestep (a few are marked by arrows) the collision is well described and energy is properly exchanged. Note the intermediate regime where the energy of the ion after the collision is sometimes higher than it was before due to the ions coming closer than their minimum approach distance because of an inappropriately long timestep. We can choose an appropriate constant timestep at the beginning of our simulations by using the following arguments. In our simulations, we inject a single hot ion with an initial speed v ini into a laser-cooled crystal whose ions have velocities that are smaller by several orders of magnitude. In such conditions one may safely assume that in the course of the simulation, no ion will acquire a velocity significantly larger than v ini . Then an estimated lower bound of the minimum approach distance d min can be obtained by setting v = v ini and d → ∞ in Eq. 4.24. With the thus-obtained d min , Eq. 4.23 with the same replacement v = v ini provides a "safe" value for a constant time step 1 .

It is easy to show that the timestep defined in this way is proportional to v -3 ini or E -3/2 ini . Since our main goal is to simulate the sympathetic cooling of ions having large initial energies E ini , this criterion leads to extremely short timesteps such as δt 1 × 10 -13 s for a 1 eV H+ in Be + , and prohibitively long computation times. Indeed a timestep of 1 × 10 -13 s would require 10 10 timesteps for 1 ms of simulation, which is the order of magnitude of our simulation times. This led us to investigate a more efficient variable timestep criterion as described in the next paragraph.

Criterion Based on Current Distance

A less stringent criterion consists in imposing that the variation of distance δd between two ions in a time step should be much smaller than the current distance d:

δd d ⇒ δd = αd, α 1 (4.25)
By the same calculations as above we get

δt = -v + √ v 2 + 2aαd a (4.26)
Eq. 4.26 gives a variable timestep criterion by using as a timestep the minimum δt calculated for all ion pairs. This variable timestep criterion is computationally expensive because it adds another O (N 2 ) computation to the already expensive Coulomb interaction calculation. It slows down the code by a factor of ∼ 3 and it leads to similar timesteps as the criterion described in Sec. 4.6.2 when the ions come close to their minimal approach distance but gives much longer timesteps at further distances, decreasing simulation times. In Fig. 4.8 we show that this variable timestep criterion leads to correct descriptions of Coulomb interactions. To do this we simulate a head-on (impact parameter = 0) 1D collision between a projectile ion and a stationary target ion from a distance d 0 = 1 mm with an initial energy E in . We simulate the interaction according to the variable timestep scheme, stop the simulation once the distance between the ions has reached 1 mm and record the energy E out of the projectile ion after the collision. We repeat this process for 20000 energies from 0.01 eV to 1 eV in a geometric progression and record the highest ratio Eout E in . The analytical result of this equal mass collision (in the limit d 0 → ∞) is that Eout E in should be 0, so by taking the highest ratio among 20000 simulations we are effectively looking at the worst case. We then repeat this procedure for different values of α in order to produce Fig. 4.8. Fig. 4.8 shows that as α goes towards 0, Eout E in goes towards the analytical result of 0. This shows that the results of the simulations tend towards the correct ones when α tends to zero. However this observation alone does not allow us to choose α because we lack a criterion to determine the maximal admissible error E out /E in . To do that we need to simulate the physical situation we are interested in (sympathetic cooling of a hot ion) and determine the required values of α and β to obtain converged result on the hot ion's energy loss. These tests are described in Sec. 4.6.5.

Interaction with the cooling laser

The interaction with the cooling laser imposes that the timestep be much smaller than the inverse of the natural life time of the transition 5.15 × 10 -8 s such that the probability of having more than one event per timestep is negligibly small. The interaction with the cooling laser also imposes that the timestep not be too small compared to the optical period because for such time scales the description of the interaction with the cooling laser should be done in the framework of the Optical Bloch Equations [START_REF] Rouse | Superstatistical velocity distributions of cold trapped ions in molecular-dynamics simulations[END_REF]. These criteria lead to Eq. 4.27.

δt << 5.15 × 10 -8 s δt >> 1 × 10 -15 s (4.27)

Energy Conservation Test

Energy conservation is a very important test of the validity of our integration algorithm. Care should be taken in the formulation of the algorithm; for example the leapfrog algorithm, which is an alternative formulation involving velocities computed at half timesteps, is not compatible with a variable timestep scheme. The velocity Verlet algorithm described above, where all quantities are defined at the same time t has no such problem. In order to check this, we decided to perform an energy conservation test on one ion oscillating back and forth in a crystal of "ghost" ions with no charge, allowing us to verify that the ion's energy is preserved by our integration with a variable timestep. Although energy conservation is not guaranteed due to the presence of a time-dependent (RF) external field, it is shown [START_REF] Landau | Mechanics[END_REF] that the ion's motion is stable for appropriate trapping parameters. The ion's amplitude of motion has an upper and a lower bound. As the ion oscillates its amplitude changes due to being in different phases of the micromotion when it reaches the top of its oscillation, but given enough time the same amplitude will be reached again.

In Fig. 4.9 we show an energy conservation test on an harmonic oscillator. Indeed, the motion along the z axis is caused by a time-independent harmonic potential and isn't coupled to the motion along x and y because there is only one ion. We see a variation of the ion's oscillation amplitude on the eighth digit over 10 ms time scales. We interpret this as merely numerical noise as when we had bugs or used Verlet integration the variations were many orders of magnitude larger.

In Fig. 4.10 we show an energy conservation test on the radial motion of an ion. As we have mentioned we expect the oscillation amplitude to vary due to micromotion but with sufficient time any oscillation amplitude should be reproduced and we find that this is the case whereas it clearly was not in bugged previous versions of the program.

Choice of α and β

Having verified the validity of our variable timestep integration algorithm, we still need to choose the values of α and β such that the RF field and Coulomb Fig. 4.9: An H+ ion is placed very close to the trap axis, about 2 mm away from the centre of the trap. The ion is trapped by the exact RF force but goes through a crystal of "ghost" Be + with no charge so that the timestep varies as it would in normal conditions. On the left we show the ion's z trajectory while on the right we zoom on the curve to show the amplitude drift with time. Other parameters as in Tab. 4.2. , q x = 0.05 for the Be + , ωz = 300 kHz and other parameters as in Fig. 4.3.

collisions are described at a sufficient level of precision. In order to find out what this required precision level actually is, we decided on a pragmatic approach. Since the main information we will extract from simulations is the capture time (as defined in Sec. 4.5), we chose α and β in order to have converged results on this quantity.

Note that this criterion depends on the initial energy of the hot ion. Indeed, as this energy increases, a larger number of collisions is required to cool the ions, which implies that each collision should be represented with higher accuracy, leading to a stricter criterion. In order to obtain safe values of α and β, we chose a high initial energy of 625 meV i.e. the highest energy we have used in our simulations. In Tab. 4.2 we show a test on the sympathetic cooling of a 625 meV mass 8 ion in a crystal of 1919 Be + . As α decreases we can see that the cooling time stays the same within the error bar of the standard deviation. This indicates that for the purposes of simulating capture times of sympathetic cooling, α = 1 5 is strict enough. A different way of testing that the sympathetic cooling process is correctly described is to study the energy transfer from the hot ion to the Be + ions. To performed this study, we turned off the laser cooling which would rapidly wipe out excess energy as it is transferred to the Be + . In Tab. 4.3 we test the convergence of the final energy as a function of α and β the parameters of the timestep criteria described in Eq. 4.26 and Eq. 4.19 respectively. α = 1 5 and β = 1 60 seem good enough for the result to be approximately converged because the mean final energy in Tab. 4.3 is the same (within error bars) for α = 1 5 , β = 1 60 and α = 1 30 , β = 1 300 . Note that β = 1 30 instead of 1 60 has much greater consequences on the convergence of the final energy than using α = 1.1 instead of 5.

In view of the results of Tab. 4.2 and Tab. 4.3 as well as some other tests not presented here (e.g: spot checking of results with stricter parameters) we decided to simulate using α = 1 5 and β = 1 60 . We also spot-checked results by simulating with smaller values of α and β when we saw unexpected results but didn't find noticeable discrepancies.

We stress that as a result of the criterion we chose, the validity of the chosen values of α and β have only been demonstrated in the strict framework of Tab. 4.3: A 625 meV mass 8 ion is launched into a crystal of 1919 Be + with the cooling laser turned off and we are interested in the final energy of the Be + after the mass 8 ion has deposited its energy into the crystal. We expect the result to be converged for small enough parameters α and β of the variable timestep criteria described in Eq. 4.26 and Eq. 4.19 respectively. a x = 5 × 10 -5 and q x = 0.05 for the Be + , ωz = 65 kHz and other parameters as in Fig. 4.3.

our studies on the capture time. They might not be valid for other purposes. For example, in Fig. 4.11 we show that for precise simulations of RF heating, β should be less than about 1 800 . However, these RF heating tests were performed with a Be + crystal that is much hotter than the Doppler limit (slightly more than 1 K). In our simulations, such as the ones presented in chapter 5, we used crystals close to the Doppler limit, so that the RF heating rate, which has been shown to follow a quadratic scaling law with temperature [START_REF] Schiffer | Temperature, ordering, and equilibrium with time-dependent confining forces[END_REF], becomes negligibly small. Furthermore what we care most about is the exchange of energy between the sympathetically cooled ion and the crystal and not a very precise description of RF heating for the Be + . Out tests such as those of Tab. 4.3, showed that β = 1 60 was enough for our sympathetic cooling simulations.

Timestep Orders of Magnitude

The choice of β = 1 60 imposes that δt <1.28 × 10 -9 s to properly describe the RF field if the trap has a frequency of 13 MHz.

In Tab. 4.4 we give orders of magnitude, taken from actual simulations, of the minimum and average timesteps during a simulation of sympathetic cooling of an H+ in a crystal of Be + as a function of the energy of the H+ ion. The d min criterion would, in theory, simulate with the minimum timestep at all times whereas the d criterion results in much longer timesteps on average. In practice the d min criterion chooses a timestep significantly lower than even the minimum timestep reached by the d criterion because the minimum approach distance in Eq. 4.24 is computed assuming a head on collision whereas in practice this is very rare. As can be seen in Tab. 4.4 the average timestep is very close to the maximum timestep allowed to represent the RF trapping field because most of the time the sympathetically cooled ion is outside of the crystal or passing in between ions. The timestep only goes down when it comes very close to one of the ions of the crystal.

Energy [meV] Minimum Timestep

Implementation

The code was written from scratch in CUDA to closely reproduce already existing Fortran code. All calculations are done in double precision on the GPU, data is only brought back to the CPU to write in log files. The Coulomb interaction is calculated exactly in O N 2 time following Nvidia's Nbody implementation [START_REF] Nyland | Fast N-Body Simulation with CUDA[END_REF]. With 15360 ions, or more, the code can exploit ≈100% of the Titan Black GPU's power giving us access to ≈1TFlop.

Low ion number performance optimisation

The O N 2 Coulomb interaction calculation is computationally expensive and requires about 15360 ions or more to expose enough parallelism to saturate the Titan Black GPU. However, the code includes an optimisation to allow lower number of particles without losing too much performance compared to the optimal number of 15360. Instead of assigning one thread per ion we assign several threads per ion, which share the summation of the forces due to the other ions, in order to reach a total of 15360 threads. In Tab. 4.5 we can see that performance of the Coulomb interaction calculation decreases much slower with decreasing numbers of ions with this performance optimisation. We chose to do most of our simulations with 1920 ions because the performance is ≈ 73% of the optimum but the 8 fold reduction of the ion number causes a 64 fold reduction in the O N 2 computational cost.

Multiple GPU implementation

The code can use all the GPUs on a computer in order to increase performance.

To do this an extra level of parallelism must be exposed, parallelism between the GPUs. To split the work between the GPUs we assign each GPU an equal fraction of the ions. Each GPU computes the Coulomb force on the ions it is responsible for due to all the other ions. Using the calculated Coulomb force, each GPU applies the Verlet Velocity update to its ions and then the updated velocities, positions and accelerations must be shared between the GPUs after the timestep, each GPU giving what it updated to the others. We implemented the data sharing both by passing through the CPU memory and by direct peer-to-peer transfers between the GPUs. At the beginning of the execution, the program empirically tests the speed of the two methods and picks the fastest. On our hardware data sharing through CPU memory is always the faster method.

For low ion numbers such as 1920 we didn't use multiple GPUs because there isn't really enough work available for two Titan Blacks, making two GPUs actually slower than one. We instead preferred to run different simulations on different GPUs. Although the multi-GPu feature reached close to perfect scaling, doubling performance with two GPUs, we didn't end up using this feature of the code very much because we preferred running faster simulations with lower ion numbers.

Available Hardware

In order to run simulations we had two machines available, with two GPUs each. One machine with two Titan Black GPUs and one machine with two Titan GPUs. The Titan and Titan Black GPUs are identical except for the Titan Black having 15 streaming multiprocessors versus the Titan's 14. Because of this we ran simulations with 1792(14 • 128) ions on the Titan GPUs and simulations with 1920(15 • 128) ions on the Titan Black GPUs.

The Titan GPU is from the Nvidia's Kepler generation of GPUs which is two generations older than the current Pascal GPUs. While our Titan GPUs can provide about 1 TFlop of double precision computation, the current Pascal GPUs can provide about 5 TFlop and the Volta generation coming later in 2017 will be able to provide about 7.5 TFlop.

The Titan GPUs are an exception in that they are from Nvidia's gamer line of GPUs, making them a five times cheaper than the usual Tesla(professional) line at around 1000 euros when they were brand new and at the time of this writing, 400 euros second hand. The Titan GPU was providing full double precision performance whereas Nvidia's gaming line of GPUs usually has chips with significantly lower double precision performance. Indeed, GPUs have at most half the performance in double precision as they do in single precision, with gaming GPUs typically having a twenty-fourth of the single precision performance in double precision. If one is satisfied with single precision computation then Nvidia's gaming line of GPUs can provide a lot of it. But for typical scientific computing in double precision, one needs to be careful that the GPU has good double precision performance. And, at least in the case of the Titan, one had to check an "Enable double precision" box in Nvidia-Settings in order to get the advertised performance.

SYMPATHETIC COOLING SIMULATION RESULTS

In this chapter we present findings from running the code for approximately the last 6 months of my PhD, previous results being tainted by bugs. Not only are the simulations long, easily taking days for a single simulation at higher energies or ion numbers, but the parameter space is large as well. Indeed, for example, we can vary the trapping parameters Ω, U 0 and V 0 and the laser parameters of angle and detuning, we can add mass-intermediate ions to the Be + along with all the combinations of these parameters. We were surprised by the large influence those parameters could have and in this chapter we present limited results from our partial exploration of this large parameter space. We present general findings about sympathetic cooling along with findings more directly relevant to GBAR's sympathetic cooling needs.

Some preliminary results on the sympathetic cooling of H + 2 and H+ ions had been previously obtained in the team using a CPU simulation code with a constant timestep scheme [START_REF] Hilico | Preparing single ultra-cold antihydrogen atoms for free-fall in gbar[END_REF]. However, at the time, due the limited computing power of CPUs, the simulations were limited to relatively low initial energies of the sympathetically cooled ion. In addition it was not possible to explore the parameter space significantly, nor to repeat the simulations a sufficient number of times to get reliable values of the capture time and error bars. The code presented in the previous chapter allows for a much deeper study of the problem, both due to the computing power of GPUs and to the variable timestep scheme we have implemented.

To begin with, we present a theoretical model of sympathetic cooling using plasma physics knowledge which we compare to our results in order to, for instance, attempt to give scaling laws between capture time, energy and ion numbers.

Theoretical Model of Sympathetic Cooling

In this section we derive an approximate theoretical model of sympathetic cooling using knowledge from plasma physics. In plasma physics exchanges of energy with macroscopic amounts of matter are well understood using meanfield approaches but the problem we are interested in is more mesoscopic and is further complicated by the fact that the sympathetically cooled ion is not in the crystal at all times, oscillating in and out if it has enough energy.

In this section we assume cooling by a single component crystal.

Energy Loss

According to Eq. (2.23) of Miyamoto's "Fundamentals of Plasma Physics and Controlled Fusion" [START_REF] Miyamoto | Fundamentals of plasma physics and controlled fusion[END_REF] the energy relaxation of an ion by collisions in a non-neutral plasma can be written

dE dt = - E τ (5.1) with 1 τ = q 2 sc q 2 lc n lc ln (Λ) 2π 2 0 m sc m lc v 3 (5.2)
with q sc the charge of the sympathetically cooled ion, q lc the charge of the laser cooled ions, n lc the density of laser cooled ions, m sc the mass of the sympathetically cooled ion, m lc the mass of the laser cooled ions, v the velocity of the sympathetically cooled ion, assuming that v is much larger than the thermal velocities of the plasma particles, and ln (Λ) the "Coulomb logarithm" of plasma physics which arises from the integral over possible values of the impact parameter. Note the strong 1 v 3 dependence of the relaxation time. In the initial stages of the cooling, the ion oscillates back and forth along the trap axis with a period T z . The speed of the ion reaches its maximum in the crystal where one can make the approximation that it is constant because the electric potential within must be constant, otherwise the laser-cooled crystal would not be in its equilibrium configuration. As the ion oscillates in the trap with a period T z , it approximately spends a fraction of dt = 2L vTz dt of its time inside the crystal. To account for this, it is natural to modify Eqs. 5.1 and 5.2 by rescaling time by a fraction vTz 2L :

dE dt = - q 2 sc q 2 lc n lc ln (Λ)EL π 2 0 m sc m lc v 4 T z = - m sc q 2 sc q 2 lc n lc ln (Λ)L 4π 2 0 m lc T z E (5.3)
Here E = 1 2 m sc v 2 is the kinetic energy of the ion while it is inside the ion crystal. Equation 5.3 shows that energy decay rate has an even more pronounced 1 v 4 dependence on v due to the ion's oscillations in and out of the crystal. This explains the tendency that we see, in for example Fig. 5.4, of the sympathetically cooled ion to lose its energy faster and faster until at some point it gets captured by the crystal very quickly compared to the total duration of the simulation.

Cooling time

Using Eq. 5.3 which describes the energy loss of the sympathetically cooled ion, we will derive some scaling laws on the cooling time of the sympathetically cooled ion.

Integrating Eq. 5.3 on both sides and neglecting the dependence of ln (Λ) on v we get

E 2 -E 2 0 2 = - m sc q 2 sc q 2 lc n lc L ln (Λ) 4π 2 0 m lc T z (t -t 0 ) (5.4)
Considering that the ion is captured when E = 0, the capture time t capt is given by

t capt = 2T z π 2 0 m lc m sc Lq 2 sc q 2 lc n lc ln (Λ) E 2 0 (5.5)
with E 0 = 1 2 m sc v 2 0 the initial kinetic energy of the sympathetically cooled ion inside the crystal and v 0 its initial velocity. We therefore expect

t capt ∝ E 2 0 (5.6)
and since the length of the laser cooled crystal L ∝ N 1 3 we expect

t capt ∝ N -1 3 (5.7)
We should note that the model derived here which resulted in the scaling laws of Eq. 5.6 and Eq. 5.7 does not take into account RF heating and is onedimensional, i.e. it assumes that the sympathetically cooled ion stays on the axis of the trap. In reality we know that the ion converts energy of motion between the 3 directions and will therefore acquire a radial oscillatory motion; this will modify the fraction of time spent inside the crystal and therefore the cooling time.

Coulomb Logarithm Value

In plasma physics textbooks, the derivation of collision cross sections involves integration over possible values of the impact parameter, leading to introduce the Coulomb logarithm, ln (Λ). Most areas of plasma physics deal with socalled "weakly coupled plasmas", where the density n is such that nλ 3 D 1, in other words there are many particles in a sphere of radius λ D , where

λ D = 0 k B T lc q 2 e n lc (5.8)
is the Debye length, the screening length in plasmas. In weakly coupled plasmas, it is natural to set the upper cutoff in the integral over the impact parameter as equal to the Debye length, leading to [START_REF] Miyamoto | Fundamentals of plasma physics and controlled fusion[END_REF] ln (Λ)

≈ λ D b 2 1 r dr ≈ ln 2λ D b (5.9)
where b is a minimal impact parameter (related to the minimal approach distance). From Eq. 5.8 we find that in a laser-cooled ion crystal, T ≈ 1 mK and n ≈ 1 × 10 12 m -3 . λ D is typically of the order of one micron, much smaller than the inter-ion distance which is of the order of tens of microns. We are thus in the opposite case of "strongly coupled plasmas". It is hard to estimate the value of the Coulomb logarithm in this regime. The Coulomb logarithm for cold Yb + ions in a radiofrequency trap was studied experimentally in [START_REF] Chen | Measurement of the Coulomb logarithm in a radio-frequency paul trap[END_REF], but this study focused on collisions between the laser-cooled ions while we are interested in collisions of a hot ion with the laser-cooled ones.

To conclude, because of the uncertainty on the Coulomb logarithm it is hardly possible to make quantitative predictions on the capture time. However, this model did provide us with the scaling laws of Eq. 5.6 and Eq. 5.7, under the assumption that the velocity dependence of ln(Λ) can be neglected.

Optimal Shape of the Crystal

The crystal responsible for the sympathetic cooling can take different shapes depending on trapping conditions. Because we initialise the sympathetically cooled ion with all its energy along the z axis one might think that a crystal which is particularly elongated along the z axis would be best to slow down the sympathetically cooled ion. However, as can be seen in Fig. 4.5, while the ion initially has energy mostly along the z axis, this energy is converted to radial (x,y) energy much quicker than total energy decay because of Coulomb interactions in the Be + crystal which push the ion in all directions.

For the case of a single component Be + crystal we investigated what the optimal shape should be for sympathetic cooling of a mass 8 ion. We chose mass 8 because this is an advantageous mass ratio with respect to the Be + compared to the 1 to 9 ratio of H+ , allowing for faster simulations. In Fig. 5.1 we can see that the optimal ellipticity ωz ωx , with ω x the secular trapping frequency, is close to 1, with a crystal slightly stretched out along the z axis, but close to a sphere. We suspect that the optimal shape at higher energy is the sphere because the capture time increases with energy, so that the initial conditions matter less and less, with the sympathetically cooled ion's energy ending up spread out evenly along all 3 axes on average.

Because of this finding we decided to do most of our simulations with spherical crystals but we cannot be sure that this is the optimal shape in all possible situations as it may depend on other parameters. Tab. 5.1: Data of Fig. 5.1. a x = 5.2 × 10 -5 and q x = 0.058 for the Be + .

Ineffectiveness of single component Be + crystal

For the purposes of the GBAR collaboration and more specifically for the capture trap we need to find a configuration where we can capture and cool ideally up to 1 eV in less than 1 s. For this reason we tried using Be + crystals to sympathetically cool H+ but our simulations have revealed that this configuration is inefficient in that the cooling works only up to a few meV of initial energy; at higher energy it becomes too slow or even non-existent.

In Fig 5 .2 we show the capture time as a function of the energy of the H+ . The capture time is already of the order of 10 ms for a few meV and in fact at around 3 meV the cooling does not always work as sometimes the H+ heats up instead of cooling, in Fig. 5.3 we show a curve from such a case and in Sec. 5.4 we discuss why we believe this happens. In effect, in the conditions of Fig. 5.2, cooling times for energies higher than ∼ 3 meV are infinite.

Impact of RF Heating

In Fig. 5.3 we show an example of a case where the ion we want to sympathetically cool actually heats up. This phenomenon is intriguing so we tried to understand what was causing this. We believe the explanation is that there is a competition between RF heating and sympathetic cooling. Since RF heating is linked to collisions one may assume that when the sympathetically cooled ion collides with an ion of the crystal, it gives part of its energy to that ion but also receives energy from the RF field, causing a competition between these two processes. To test our hypothesis we simulated the cooling of H+ under the pseudopotential approximation where RF heating isn't present.

We found that under the pseudopotential approximation we were able to cool H+ with much higher initial energies. Indeed, in Fig. 5.4 we show the capture of a 21 meV H+ by a crystal of 1919 Be + under the pseudopotential approximation. We've found that at this energy the average capture time is about 3.2 ms in the pseudopotential approximation while cooling doesn't occur with the exact RF potential. Furthermore, with the exact RF potential, it takes ≈10 ms to capture a few meV as shown in Fig. 5.2 which is more time than it takes to cool a 21 meV H+ ion in the pseudopotential approximation. We believe this is a strong indication that there is a competition between RF heating and sympathetic cooling, and that this competition is easily lost by sympathetic cooling in the case of H+ being cooled by a crystal of Be + , probably because of the unfavourable mass ratio of 1 to 9. We did not notice as strong of an RF heating effect with more favourable mass ratios. Note that one might try to reduce the adverse effect of RF heating by using a lower trapping parameter q x . The dependence of sympathetic cooling on q x is discussed in Sec. 5.8.

Improvement Using a Two Component Crystal

An idea [START_REF] Zhang | Molecular-dynamics simulations of cold single-species and multispecies ion ensembles in a linear paul trap[END_REF] to improve sympathetic cooling efficiency is to use auxiliary ions with a mass in between that of the sympathetically cooled ion and the Be + . Indeed, we have found that using a two component crystal of Be + and HD + was a lot more effective for the sympathetic capture and cooling of H+ . The HD + is sympathetically cooled by the Be + and serves as an intermediary of better mass ratio, 3 to 1, for the sympathetic cooling between the Be + and the H+ . Numerically we confirm that this idea works and in Fig. 5.5 we show capture times of H+ in a crystal of 1441 Be + and 350 HD + . Keeping in mind the 1 s limit on the lifetime of the H+ in the capture trap Fig. 5.5 shows that in this configuration the highest acceptable energy for the H+ is about 125 meV. In Fig. 5.6 we show an example curve of a 125 meV H+ being cooled as in the conditions of Fig. 5.5. In Fig. 5.6 we can see that the auxiliary HD + ions remain hotter than the Be + ions, which is expected because they are being cooled indirectly by the Be + ions while being heated up by the H+ ion and RF heating. In the left part of Fig. 5.7 we show the crystal used in the simulations of Fig. 5.5. The reasons why we used a crystal of such peculiar shape are explained in the next section. H+ sent into a crystal of 1441 Be + and 350 HD + with stability parameters q x = 0.2 and a x = 0.0014 for the H+ , ωz = 250 kHz and other parameters as in Fig. 5.2. The error bars are the standard deviation divided by the square root of the number of simulations. The solid line is a fit of the form c 1 x 2 + c 2 with c 1 = 0.015 ± 0.004 ms meV -2 and c 2 = -2 ± 1.5 ms. Eq. 5.5 predicts c 1 = 0.046 ms meV -2 with ln Λ = 1 and taking the mass, density and charge of the auxiliary ions. The dashed line is an identical fit that does not weight error bars with c 1 = 0.0274 ± 0.0003 ms meV -2 and c 2 = -11 ± 2 ms. We comment on the fit in Sec. 5.9.

Trapping Parameter a x Disfavouring Orbits

To do the simulations in Fig. 5.5 we had to break the rotational symmetry of the trap quite strongly using a high DC voltage U 0 , leading to a high trapping parameter a x . Indeed we found that a lot of simulations would get stuck in a configuration where the H+ would orbit around the crystal indefinitely when the a x parameter was lower. A strong asymmetry, with secular frequencies for the Be + of f x = 117 kHz, f y = 19.5 kHz and f z = 83.3 kHz, seems to prevent those indefinite orbits from occurring in this situation. Due to their smaller mass the auxiliary HD + ions are much more symmetrically trapped with secular frequencies of f x = 323 kHz, f y = 253 kHz and f z = 144 kHz while for the H+ f x = 938 kHz, f y = 872 kHz and f z = 250 kHz. For reasons unknown to us, we did not observe sympathetically cooled ions going into stable orbits in other simulations such as those with single component crystals. In Fig. 5.8 we show example curves of the H+ reaching a stable orbit with identical trapping conditions except that the dissymetrising voltage U 0 = 0.275 V is slightly lower than the value of 0.3 V used to obtain the results of Fig. 5.5-5.7.

Be + /Auxiliary ion balance

Having in mind the goal of finding the best sympathetic cooling configuration for GBAR, we tried to improve upon the results of Fig. 5.5 by increasing the proportion of HD + ions to Be + ions. Hoping that the HD + 's better mass ratio of 3 to 1 would cool the H+ faster. Macromotion energy of the three species involved. a x = 0.0013 for the H+ and other parameters as in Fig. 5.2. Note that the HD + becomes cooler than the average Be + ion because the HD + , having a lower mass than the Be + , crystallises close to the axis of the trap and therefore experience a smaller RF field. The temperature of the HD + is actually given by the temperature of the coolest Be + ions located closer to the trap axis. Keeping the total number of ions in the two component crystal to 1791, we tried a crystal of 991 Be + and 800 HD + . We found that at an initial energy of 62.5 meV the cooling is indeed more efficient. If we take the error bar of the average capture time to be the standard deviation divided by the square root of the number of simulations, we find that the capture time for the 991 Be + /800 HD + crystal is 43.1 ms ± 5.7 ms against 96 ms ± 16 ms for the 1441 Be + /350 HD + crystal.

However, we found that having this many auxiliary ions made the configuration unstable for the capture of 125 meV H+ . The crystal would suffer runaway heating and loss of the ions, preventing the cooling of the H+ . In Fig. 5.9 we show an example curve of this happening with conditions otherwise identical to those used to generate the data of Fig. 5.5. In Fig. 5.9 we see that the auxiliary ions in blue heat up as the H+ in green loses its energy. Unfortunately the Be + isn't able to cool down the auxiliary HD + fast enough and the HD + gets so hot that it completely melts the Be + crystal past a point of "no return" as the laser detuning is not appropriate to cool such hot Be + and the RF heating increases due to its temperature dependence. After around 500 ms the purple curve goes to 0 which corresponds to the loss of all the Be + in the trap. In contrast, the crystal of 1441 Be + and 350 HD + (Fig. 5.6) behaved much better as, although it was slow, the 10 simulations we performed all led to the capture of the H+ . We also found similar behaviour when using H + 2 as an auxiliary ion instead of HD + . We believe this is because H + 2 has a worse mechanical coupling with Be + and therefore is cooled less efficiently which makes it more likely to manifest runaway heating. In Fig. 5.10 we show an example of sympathetic cooling not working due to the crystal heating up more and more. In Fig. 5.11 we show the capture time and its standard deviation as a function of the energy of the In light of our findings we abandoned the idea of using H + 2 as an auxiliary ion and stuck to the proportion of 350 HD + to 1441 Be + . But we noted that a higher proportion can significantly speed up cooling if it is stable which leads us to believe this could be an important parameter, experimentally and numerically, to find a better configuration for the cooling of high energy H+ as required by GBAR.

Capture Time q x Dependence

In view of the RF heating process evidenced in Sec. 5.4 it is obviously important to investigate the influence of the trapping parameter q x . We have found a strong dependence of the capture times on this parameter q x . We mostly studied the case of sympathetic cooling by a single-component crystal of Be + , also varying the mass of the sympathetically cooled ions: m = 1 u ( H+ : see Tab. 5.2), m = 2 u (H + 2 : see Tab. 5.3) and m = 8 u (see Tab. 5.4). In order to assess the role of RF heating in the observed behaviour, we ran simulations under the pseudopotential approximation for m = 2 u (Tab. 5.5) and m = 8 u (Tab. 5.6). Finally, we also investigated the case of a two-component crystal of Be + and HD + (Tab. 5.7). All the data of Tab. 5.2-5.7 is plotted in Fig. 5.12-5.14.

Before commenting on these curves, it is interesting to derive the predicted dependence in the framework of the model presented at the beginning of this chapter. Assuming a spherical crystal (ω x = ω y = ω z = ω), which is close to what we use (except for the two-component Be + /HD + crystal, see Fig. 5.7), according to Eq. 5.5 the capture time scales like

t capt ∝ T z n lc L (5.10)
Using the fact that T z ∝ 1 ω ∝ 1 qx (neglecting a x in Eq. 1.20), n lc ∝ ω 2 ∝ q 2

x (Eq. 1.34), and L ∝ n

-1 3 lc ∝ q -2 3
x , one obtains

t capt ∝ q -7 3
x .

(5.11)

Quantitatively, we have

t capt = 32 √ 2π 2 4 3 0 m 1 3 lc 3 2 3 Ω 7 3 q 7 3
x m sc N

1 3
Ions q 2 sc q 2 3 lc ln (Λ)

E 2 0 .
(5.12)

For the sympathetic cooling of a m = 8 u ion (Fig. 5.12) one observes a decrease of the capture time with increasing q x . A q -7 3

x power-law fits the data well (Fig. 5.15). The curve obtained under the pseudopotential approximation (all other parameters being identical) is in good agreement, indicating that RF heating is a negligible effect in this case.

However, in case of a m = 1 u ion (Fig 5 .13), the variations of the capture time vs q x are qualitatively very different. For a H+ sympathetically cooled by a single-component Be + crystal, reveals an optimum around q x = 0.1. It is worth noting that the study of the capture time dependence on energy (Fig. 5.2), for which the value q x = 0.052 was used, was thus not done in optimal conditions (however, the optimal q x value may depend on the initial energy of the sympathetically cooled ions). In the case of cooling by a twocomponent Be + /HD + crystal, an optimum is found at q x = 0.2 (although less pronounced than in the previous case), which is the value chosen for the study of the capture time dependence on energy (Fig. 5.6). This behaviour could be interpreted as resulting from a competition between the t capt ∝ q -7 3 x law and RF heating, which becomes increasingly important with increasing q x .

For a m = 2 u ion (Fig 5 .14), an optimum is found at q x = 0.25 which could be interpreted along the same lines. In order to assess the effect of RF heating, we performed simulations under the pseudopotential approximation, which reveal a very different and surprising behaviour. Above a threshold value (q x ∼ 0.16) the capture time suddenly increases by a large factor, and becomes much longer than with the time-dependent trapping field. Although the exact reasons of this behaviour are still unclear, we have noticed that the parameter U 0 can have a crucial influence on the capture time. Indeed the q x = 0.155 point in Tab. 5.5 was done with U 0 = 0.005 V and U 0 = 0.05 V. We found that with U 0 = 0.005 V the capture times were segregated into two groups of a few ms and ∼ 30 ms capture times. This corresponds to the effect already discussed in Sec. 6.6: for small values of U 0 , the system is close to rotationally symmetric, and the sympathetically cooled ion often ends up orbiting around the Be + crystal. The stepwise dependence of the capture time could be due to Tab. 5.2: Capture time and its standard deviation of a 3.125 meV H+ in a crystal of 1919 Be + for different q x parameters. The q x parameter given is that of the H+ . As the radial confinement parameter q x was changed, the axial confinement parameter ωz was adapted to maintain a spherical shape of the crystal. Other parameters as in Fig. 5 and its standard deviation in a cloud of 1919 Be + for different q x parameters. The q x parameter given is that of the H + 2 . As q x is varied ωz is adapted as well to maintain a spherically shaped cloud. a x = 1.2 × 10 -4 with other parameters as in Fig. 5.2.

the fact that as q x increases while U 0 is kept fixed, the x-y asymmetry (defined e.g. as 2 |ωx-ωy| ωx+ωy ) decreases; the orbiting of the sympathetically cooled ion might appear below some threshold value of the asymmetry. Further simulations are required in order to check this hypothesis.

Having found a strong dependence of capture times on the trapping parameter q x we believe that any sympathetic cooling experiment should optimise this parameter empirically and possibly using simulations.

Capture Time vs Energy Scaling Law

A rough scaling law between capture time and the initial energy of the sympathetically cooled ion would be very useful, allowing extrapolations to higher energies. Eq. 5.5 predicts a quadratic scaling law but doesn't take into account RF heating for example. We chose to fit our data with a law of the form c 1 E 2 + c 2 , where c 2 is expected and is found to be negative, in order to take into account the finite cooling time from our arbitrary 2.5 K limit to the Doppler limit. In Fig. 5.2 (Be + / H+ ), Fig. Tab. 5.4: Capture time of a 21.25 meV mass 8 ion and its standard deviation in a crystal of 1919 Be + for different q x parameters. The q x parameter given is that of the mass 8 ion (a x = 3 × 10 -5 ). As q x is varied ωz is adapted as well to maintain a spherically shaped crystal. Other parameters as in Fig. 5 Tab. 5.5: Capture time of a 21.25 meV H + 2 and its standard deviation in a crystal of 1791 Be + for different q x parameters, under the pseudopotential approximation. The q x parameter given is that of the H + 2 . As q x is varied ωz is adapted as well to maintain a spherically shaped cloud. a x = 1.2 × 10 -4 with other parameters as in Fig. 5 Tab. 5.6: Capture time of a 21.25 meV mass 8 ion and its standard deviation in a crystal of 1791 Be + for different q x parameters, under the pseudopotential approximation. The q x parameter given is that of the mass 8 ion. As q x is varied ωz is adapted as well to maintain a spherically shaped crystal. The <60 µs capture times are due to a 60 µs resolution of the data which was designed for ms time scales. Other parameters as in Fig. 5 Tab. 5.7: Capture time of a 31.25 meV H+ and its standard deviation in a crystal of 1441 Be + and 350 HD + for different q x parameters. The q x parameter given is that of the H+ . As q x is varied ωz and a x are adapted to try and maintain a similarly shaped crystal. Other parameters as in Fig. 5.2. Fig. 5.12: Capture time versus q x for a mass 8 ion with an initial energy of 21.25 meV. When the radial stability parameter q x is changed, the axial trapping is adapted to keep the shape of the crystal constant. The "BeM8 pseudo" curve corresponds to the data of Tab. 5.6, "BeM8" to Tab. 5.4.

Error bars are taken to be the standard deviation divided by the square root of the number of simulations.

Fig. 5.13: Capture time versus q x for an H+ ion. When the radial stability parameter q x is changed, the axial trapping is adapted to keep the shape of the crystal constant. The "BeH" curve corresponds to the data of Tab. 5.2 (3.125 meV) while "BeHHD" curve corresponds to the data of Tab. 5.7 (31.25 meV). Error bars are taken to be the standard deviation divided by the square root of the number of simulations. Fig. 5.14: Capture time versus q x for an H + 2 ion initial energy of 21.25 meV. When the radial stability parameter q x is changed, the axial trapping is adapted to keep the shape of the crystal constant. The "BeH2 pseudo" curve corresponds to the data of Tab. 5.5 while the "BeH2" curve corresponds to the data of Tab. 5.3. Error bars are taken to be the standard deviation divided by the square root of the number of simulations. x law. For the pseudopotential data we find c 1 = 0.005 ± 0.002 ms and c 2 = -1.8 ± 0.1, and for the exact potential data we fixed c 2 = - 7 3 and find c 1 = 6.9 × 10 -4 ± 8 × 10 -5 ms. We could not fit while weighting the error bars because of the undefined error bars on some data points (see caption of Tab. 5.6). Eq. 5.12 predicts 3.17 × 10 -3 ms with ln Λ = 1. of the energy of the sympathetically cooled ion. Mostly looking at the fits without weighting the error bars because the first points have artificially low error bars due to simulations being faster and therefore more numerous. We find the best agreement in Fig. 5.11 and Fig. 5.5 with more mediocre fits in the 2 other figures. The fits of Fig. 5.2 and Fig. 5.16 pass through the first few points of data very well but are under the points at higher energy indicating a higher than quadratic dependence. The fact that our model underestimates the increase of the capture time with energy could be explained by the fact that it doesn't take into account RF heating, which is known to increase with energy and by conversion of axial energy into radial energy, potentially causing the hot ion to orbit around the crystal. Fig. 5.3 for example is in contradiction with the quadratic scaling law because it is an example of seemingly infinite cooling time as the ion heats up instead of cooling down.

The quadratic scaling law between capture time and energy of Eq. 5.5 is correct in some situations, where the assumptions used in its derivation are roughly verified and could be used as a rule of thumb but we have shown its shortcomings. In the case that interests us the most, as it is the most promising scheme for the GBAR experiment, Fig. 5.5, we have a good agreement and we will make use of it in Sec. 5.12.

Lack of Improvement using a Hot Cloud

In the plasma physics literature [START_REF] Zwicknagel | Stopping of heavy ions in plasmas at strong coupling[END_REF] one can find the idea that the stopping power of a plasma depends on its temperature. The optimum depends on the exact conditions but order of magnitude is that the speed of the sympathetically cooled ion should be of the same order as the thermic velocity v th = k B T m with m the mass of the ions making up the laser-cooled crystal and T its temperature. In Fig. 3, on page 58 of [START_REF] Zwicknagel | Stopping of heavy ions in plasmas at strong coupling[END_REF] one can see that the stopping power shows an optimum at a few v th and decreases monotonously on both sides. In our case, the ∼mK Be + 's thermal velocity is much lower than that of the hotter and lighter H+ so this line of reasoning led us to investigate the effect of increasing the temperature of the Be + crystal, for example by increasing the detuning and/or intensity of the cooling laser.

We tried this to the point of having Be + barely one order of magnitude cooler than the sympathetically cooled H+ but were unable to find notable decreases in capture times. This is due to the fact that previous studies of the stopping power of a plasma treated a very different situation, where the plasma particles (typically electrons) are much lighter than the particle to be cooled. In our case, requiring the velocity of the sympathetically cooled ion to be of the same order of magnitude as the thermal velocity of the laser-cooled ions essentially comes to requiring it to be already colder than the laser-cooled ions.

Detection through Fluorescence Signal

The sympathetically cooled ion does not interact with the cooling laser and therefore it cannot be seen by a camera. Ensembles of sympathetically cooled ions usually appear as a dark area within the laser-cooled crystal ; single ions within a few-ion string are also easily identified. However, the identification of a single "dark" ion in a crystal of several thousand laser-cooled Be + ions is a challenging task. In this section we discuss a possible detection scheme revolving around the idea of detecting a variation of fluorescence from the crystal due to the heating influence of the sympathetically cooled ion. Such a diagnostic would be extremely useful during the running of the GBAR experiment.

If the sympathetically cooled ion were to significantly heat up the crystal one may be able to detect a noticeable increase in average fluorescence count, or even a decrease if the ions get so hot that the laser detuning doesn't allow for much interaction. However, in our configuration where a light ion deposits low amounts of energy every time it passes through a crystal composed of thousands of ions, we were not able to notice such an effect in the "raw" fluorescence signal. In Figure 5.17 we give an example of a fluorescence curve, which appears to be just noise around a mean value.

However if we look at the Fourier transform of the data in Fig. 5.17 we find a peak at around 30 kHz. The two very thin peaks at 16 kHz and 34 kHz are unaffected by the absence of the sympathetically cooled ion but the peak at 30 kHz has a much lower amplitude if the sympathetically cooled ion isn't oscillating through the crystal. We also found that this peak was present in cooling of H+ in the crystal of 1469 Be + and 450 H + 2 from Fig. 5.11. This leads us to conclude that this peak is indeed an indication of the presence of the H+ which could in principle be accessible experimentally.

In order to be able to take the Fourier transform of "raw" fluorescence signals it was important for us to have fluorescence data equally spaced in time. In order to achieve this the code would, when necessary, take a smaller timestep in order to be able to save fluorescence data which was as equally spaced in time as double precision would allow (16-18 digits).

We have identified the 30 kHz frequency to the macromotion frequency of the Be + along the y axis by taking the Fourier transform of a Be + 's trajectory. The y macromotion frequency for one ion in the trap would be 19.5 kHz but that is disturbed by the space charge of the other ions. Approximating the HD + crystal in between the Be + (shown in Fig. 5.7) to be a spheroid and neglecting the space charge of the other Be + ions we calculate a peak at 33 kHz which is close to the observed 30 kHz. The motion along the y axis modulates the fluorescence although the laser is along the z axis because our code takes into account the Gaussian profile of the beam. Motion along y therefore modulates the cooling laser intensity the Be + sees. Our interpretation is that the H+ heats up the crystal, exciting this mode, which then influences the height of this peak in the Fourier transform of the fluorescence data.

Interestingly when the stability parameter a x is lowered, the peak moves to higher frequencies and its magnitude seemed to decrease, indicating that the strength of this peak may be specific to the peculiar trapping parameters used, which produced the crystals depicted in Fig. 5.7. Perhaps there are configurations that lead to a bigger peak.

We found that looking at the Fourier transform of different samples of fluorescence data leads to significantly different peak heights. We believe this is due to the random variations of the energy exchange between the sympathetically cooled ion and the ion crystal. With this intuition in mind, we wrote code to determine the height of the peak as a function of time. This is done by a kind of sliding Fourier transform: for every point of fluorescence data, the Fourier transform is performed on 1 ms of data around that point and the height of the peak is recorded as the highest value seen between 29 kHz and 33 kHz. The first and last 0.5 ms of the data are ignored as the 1 ms window around them doesn't exist. It should be noted that with 1 ms of data the two thin peaks at 16 kHz and 34 kHz are not resolved. In Fig. 5.18 we show such plots of the height of the 30 kHz peak as a function of time.

We see in Fig. 5.18 that without the H+ depositing energy in the crystal, the height of the peak stabilises to a baseline value. But in the presence of the H+ the height of the peak varies dramatically up to 15 times the baseline. It might therefore be feasible, in the GBAR experiment, to perform Fourier transforms of the fluorescence data in real time in order to detect the presence of the H+ . However it will have to be determined whether or not this diagnostic works experimentally as our data assumes no noise and 100% collection efficiency.

Capture Time vs Ion Number Scaling Law

Eq. 5.7 predicts that the capture time scales like N - 1 3 Ions . Meanwhile Eq. 1.36 tells us that the "energy tolerance" of an ion crystal, i.e. the maximum energy Both plots are given in units of the baseline height of the peak when the crystal is in equilibrium. For a short time the peak is higher than the baseline on the left-hand plot because of the non-stationary state the crystal is in before reaching its equilibrium. Parameters are identical to Fig. 5.17 with the data of the right-hand plot being from that figure .   of an incoming ion for capture to be feasible, should scale like N Ions . A scaling law between ion number and capture time would be very interesting for GBAR in order to extrapolate what capture times might be with a larger crystal than what we use numerically. Indeed, we have to keep ion numbers relatively low in our simulations because of the O(N 2 Ions ) scaling law on simulation times. We explored this scaling law in the case that was most relevant to GBAR. We used the highest energy with the best cooling configuration we found and redid the simulations with more ions while trying to keep the configuration close to the original. Results are given in Tab. 5.8 and plotted in Fig. 5. [START_REF] Major | Charged Particle Traps: Physics and Techniques of Charged Particle Field Confinement[END_REF]. We find that the scaling between capture time t capt and N ions is approximately t capt ∝ N - 1.39 Ions which is very different from Eq. 5.7. If we extrapolate to a crystal of 36500 ions it would take 6.5 ms to capture a 125 meV H+ . Further assuming that the capture time scales quadratically with energy according to Eq. 5.6, and as verified in Fig. 5.5, we can extrapolate that a 1 eV H+ would be captured in about 400 ms i.e. the same time as required with an initial H+ energy of 125 meV and N ions = 1791. This is in agreement with the N Tab. 5.8: Capture time of a 125 meV H+ in crystals composed of different ion numbers with approximately 20% HD + and 80% Be + . Other parameters as in Fig. 5.2. As the ion numbers were increased we had to decrease a x to avoid the elongated shape shown in Fig. 5.7 from touching the electrodes. One simulation, shown in Fig. 5.20, with 2883 Be + and 700 HD + did not lead to capture and was left out of the statistics. One simulation with 5829 Be + and 1850 HD + also didn't lead to capture as the H+ went into orbit. We believe we didn't find optimal trapping parameters for the two configurations with higher ion counts as iterations are prohibitively computationally expensive. Ions with the fitting parameters found to be c 1 = 1.3 × 10 7 ±9 × 10 6 ms and c 2 = -1.38±0.09. ωz = 250 kHz, q x = 0.2 for the H+ and other parameters described in Tab. 5.8.

CONCLUSION

REMPI Ion Source

In this thesis we have successfully demonstrated the trapping of H + 2 ions produced in a state-selective way by multiphoton ionization (REMPI). However, the trapped ion signal was orders of magnitude lower than expected, forcing us to use a supersonic beam with a backing pressure of 1 bar instead of the planned effusive beam with a backing pressure of 1 mbar. We believe this is due to not saturating the REMPI transition through lack of optical power and/or inadequate focusing due to the spatial mode of the 303 nm laser. We have made improvements to the ion source but have been unable to test it since then. In particular, a flexible bellow allowed us to align the H 2 molecular beam with the centre of the trap and the pulsed valve is now closer to the skimmer which should increase ion production.

Further work could involve switching from the 3+1 REMPI process at 303 nm to a 2+1 photon process at 202 nm through the E, F state of H 2 [START_REF] Perreault | Angular and internal state distributions of H + 2 generated by (2 + 1) resonance enhanced multiphoton ionization of H 2 using time-of-flight mass spectrometry[END_REF]. Indeed, the larger transition strength of the X → E, F two-photon transition (as compared to the X → C three-photon transition) would allow a similar transition probability with a larger focal spot and thus a much larger production volume. This could be achieved rather easily in our laser setup by summing the 606 nm light with the with the already-existing second-harmonic 303 nm light currently used for REMPI in a single pass through a non-linear crystal. It may also be desirable to replace the pulsed valve which we found out was not bakeable. Finally, this ion source will need to be mounted next to the main chamber with the linear Paul trap in order to carry out the sympathetic cooling of H + 2 in a crystal of laser-cooled Be + ions, and finally the high-resolution spectroscopy experiment.

We have shown that we were able to understand the dissociation of H + 2 created by electron-impact ionization. Further work will involve doing the same tests on ions created by our REMPI ion source in order to double-check that the ions are indeed created in ν = 0.

Sympathetic Cooling Simulations

We have presented our work of writing, testing and exploiting CUDA code to simulate the dynamics of sympathetic cooling on GPUs. We have shown that the capture trap of the GBAR experiment can hardly use a single component Be + crystal because the maximum energy such a crystal can capture and cool isn't very high. Indeed, we have shown that there is a competition between cooling and RF heating which is easily lost for such an unfavourable mass ratio of 1 to 9. A somewhat unexpected and counter-intuitive result is that in such situations the (supposedly) sympathetically cooled ion ends up being heated. Instead, we have shown that a highly rotationally-asymmetric crystal composed of about 80% Be + and 20% HD + is a promising solution. The HD + serves as an intermediate mass to improve the mechanical coupling between the H+ and the laser-cooled Be + . We found that the high rotational asymmetry of this configuration prevented the H+ from going into orbit around the crystal. We showed that this configuration could also allow experimental detection of the H+ via the Fourier transform of the Be + fluorescence signal. This could be a very useful experimental tool, especially in the early stages of the experiment. We have shown that the stability parameter q x strongly influences cooling times and briefly discussed the optimal shape of the cooling crystal, finding an optimum around the sphere. We discussed scaling laws between capture time, initial energy of the sympathetically cooled ion and ion numbers in the cooling crystal. The t capt ∝ E 2 law between capture time and initial energy derived from a simple plasma physics model is verified in some cases but, as the model doesn't take into account, among other things, RF heating, it is too optimistic and the dependence sometimes seems to be stronger than quadratic at the highest energies we investigated. However, the dependence of the capture time on the ion number predicted by the same model (t capt ∝ N -1 3 ions ) is found to be strongly underestimated as in the case of H cooled by a Be + /HD + crystal we found something closer to t capt ∝ N -4 3 ions . Extrapolating from the t capt ∝ E 2 and t capt ∝ N -4 3 ions laws we estimate that GBAR's goal of cooling H with an energy dispersion of ≈ 1 eV might well be possible with a crystal of ∼ 36500 ions, but only carrying out the experiment will bring a definitie conclusion.

There are many things that could be studied with the code we developed. For example, our latest simulations have shown the importance of dissymetrising the trap by applying a DC voltage U 0 to the RF electrodes. The optimization of H+ cooling by a single-component crystal should probably be redone in light of this new finding to gain a deeper understanding of the influence of the trapping parameters on sympathetic cooling. A very important step will consist in using our code to study and optimize the H+ ion injection procedure, exploiting the code's feature of importing Simion potential file. One way to inject is to temporarily ground one of the endcap electrodes and it will be important to study how the ion crystal in the trap will be affected by this operation.

Finally, the ion detection method via the Fourier transform of the fluorescence signal we propose could probably be improved, by looking for ways to increase the height of the peak in order to make it easier to detect with experimental noise and low collection efficiency. These insights may not only lead to define a satisfactory experimental arrangement for the capture of antimatter ions in GBAR, but be of general use for any experiment aiming to trap and cool ions produced in external sources.

INTRODUCTION TO GPU PROGRAMMING

For my PhD I had to learn GPU programming in order to exploit the computational power of GPUs for my simulations. This is not part of standard education for a physicist and in this chapter I hope to transmit some of this knowledge. Here are some, among many, references one could give to learn C++ [START_REF] Lippman | C++ Primer[END_REF][START_REF] Stroustrup | The C++ Programming Language[END_REF] and CUDA [START_REF] Nvidia | CUDA Udacity course[END_REF][START_REF] Cook | CUDA Programming: A Developer's Guide to Parallel Computing with GPUs[END_REF].

B.1 Introduction

GPU stands for Graphics Processing Unit because their original purpose is to do the computationally expensive graphics rendering operations for video games. The computational cost of video games is usually dominated by graphics rendering operations. In a typical computer the GPU is one order of magnitude faster than the CPU in terms of floating point operations per second (flops). This computational power is the reason more and more code is being run on the GPU. A notable example is the rising field of machine learning which has grown to a point where GPU manufacturers are making architecture changes specifically targeted to that field such as half-precision 16 bit floating point operations. We can also exploit the power of the GPU for the purposes of physics simulations.

B.1.1 What makes the GPU faster than the CPU?

To answer this question it is important to understand that computational power is a question of power efficiency. A chip is allowed approximately 200 W of power that can be dissipated by the cooling system and increasing power consumption is prohibitively expensive for data centres anyway.

Suppose you have 1 CPU core and you want to double the flops output. One way is to double the clock frequency of the chip, passing through it twice as much current and therefore, by Joule heating, quadrupling the power consumption. An other way is to simply use 2 CPU cores, doubling the power consumption. Making chips with lower clock frequencies and higher number of cores is therefore more energetically efficient and at equal power consumption allows for more computational power.

In a standard desktop computer nowadays a CPU has 2 to 4 cores. More expensive CPUs for clusters have higher core counts reaching numbers such as 12. For comparison the latest chip from Nvidia has 3584 cores. Higher core count at lower clock frequencies is the essence of how GPUs achieve higher power efficiency. However this changes the programming paradigm from single core sequential code to massively parallel code, which depending on the task might be hard to achieve, some problems being inherently sequential.

The point is also made that the GPU is optimised for throughput while the CPU is optimised for latency. The CPU is much quicker to execute one simple instruction than the GPU as it takes a few µs to execute anything on a GPU but the GPU can, per second, execute a lot more instructions than a CPU.

B.1.2 Faster than GPUs

There are two types of chips that can be even more energy efficient than GPUs.

ASIC

The Application Specific Integrated Circuit (ASIC) is a chip designed for one and only one task. It is therefore a lot more efficient at this task. Notable examples of the use of ASICs are crypto-currency mining and the GRAPE [START_REF] Makino | Astrophysics with grape[END_REF] series of astrophysics NBody calculating ASICs. ASICs can be an order of magnitude or two faster than a GPU.

FPGA

The Field Programmable Gate Array (FPGA) is somewhere in between the general purpose CPU and the application specific ASIC. FPGAs are chips with reprogrammable connections. After it has been programmed for a certain task the FPGA behaves like an ASIC. FPGAs have a high performance per Watt but they tend to be lower powered, in the end making them slower than a GPU and requiring several in parallel to make up for the lower power which is what was done with GRAPE-9 [START_REF] Fukushige | Hierarchical tree algorithm for collisional N-body simulations on grape[END_REF].

B.1.3 GPU programming languages

The architecture of the GPU being so different from that of a CPU, GPU programming has to be done with specific languages.

CUDA

CUDA is a GPU programming language which is proprietary to Nvidia and only works for Nvidia GPUs and therefore does not work for AMD GPUs. It is based on C++ with language extensions allowing control of the GPU.

OpenCL

OpenCL is not proprietary, works for all GPUs, CPUs and Xeon Phis and is capable of extracting about as much performance from a Nvidia card as CUDA. It is based on C with language extensions allowing to control GPUs. It is a lot more verbose and although OpenCl code can be transferred from a multi-core CPU to a GPU, one should not expect the performance to do so, the architectures are different and therefore the code needs to be different in order to fully exploit the different hardware.

OpenACC/OpenMP

OpenACC is a set of compiler directives with the goal of being able to transfer CPU code to the GPU by adding "#pragma" directives in the parts of the code that should be run on the GPU. This works but typically doesn't allow the same level of performance as a lower level CUDA/OpenCL implementation. OpenMP now also supports such directives since version 4.0.

B.2 Introduction to C++

In order to introduce CUDA programming I must first introduce C++ as CUDA is based on C++. Here I assume basic programming knowledge as a physicist usually knows Matlab, Python and/or Fortran.

B.2.1 Variable declaration

Unlike Fortran, a variable can never be declared implicitly in C++, the Fortran equivalent of Implicit None. A variable is declared by the following syntax:

double a ;
where "double" is a double precision floating point number but could have been any other type and "a" is the name I chose for the variable. Variable declaration and initialisation can be done on the same line using the following syntax:

double a =1;
where 1 is the value I arbitrarily chose for a's initial value. An equivalent syntax is available since C++11. double a { 1 } ; This syntax using braces also works for arrays and objects and is therefore recommended since C++11 in order to unify syntax among other things although = is still the most widespread syntax for obvious reasons.

Variable scope

In C++ a variable has a scope which is limited to the nearest braces {}. This is often the braces delimiting the function you are coding in. Outside these braces the variable doesn't exist and using it will result in a compilation error. If a variable is declared outside any function it is referred to as a global variable which can be accessed anywhere (except the GPU in the case of CUDA).

Constexpr

Since C++11, which is used in the most recent versions of CUDA, one can declare a variable with the prefix constexpr. Which means that the variable is a constant known at compile time and is therefore usable on both the GPU and CPU because it is equivalent to copy pasting the value in the code. c o n s t e x p r double PI = 3 . 1 4 ;

B.2.2 Functions

While in FORTRAN there are routines and sub-routines, in C++ there are only functions. A function has a return type and parameters and is declared using the following syntax: int MyFunction ( parameter1 , parameter2 ) { int r e s u l t ; //Work g o e s h e r e return r e s u l t ; } Here I gave an example of a function which returns an integer. A function may not actually return anything and this is specified to the compiler with the void keyword.

void MyFunction ( parameter1 , parameter2 ) { cout << " H e l l o World" << e n d l ; }

B.2.3 Pointers and References

A pointer is an elementary object in C/C++ and assembly language. One can think of it as merely an integer number N which "points" to byte number N in memory. Pointers are necessary for the program to remember where it has stored data in RAM. The reason a "32 bit" machine can only have 4Gb of RAM is because a 32 bit pointer can only have 2 32 ≈ 4 × 10 9 different values. A pointer has a type associated to it. For example, int* is a pointer to an integer. References also have to do with memory indexes. The reference operator is denoted as &. If I have a variable a then &a is a reference to a which can be thought of as its memory address. If I have an integer i=2 and I assign a pointer to be equal to its reference, that pointer is now pointing to i and I can change the value of i through that pointer. To do that I can use the dereference operator *. If I have a pointer p then *p is the value of whatever p was pointing to. int i =2; int * p ; p=&i ; * p=3; // i now has a v a l u e o f 3

B.2.4 Arrays

In C/C++ an array A is represented by a pointer to its first element. Knowing the type of the pointer the compiler knows which memory address to go to if you request the Nth element of the array A[N-1], the -1 being due to 0-based indexing. In order to allocate memory for an array you have to call malloc() or new. And to free this memory you have to call free() or delete.

int * A=new int [ 1 0 ] ; // A l l o c a t e an a r r a y o f 10 i n t e g e r s . . . // Work w i t h a r r a y

Cores, Warps

Cores are in groups of 32 known as Warps. The cores in a warp share an instruction fetching unit, which has the disadvantage that all cores in a warp must be executing the same instruction per clock cycle. If half of the cores in a warp want to do something different from the other half in a given clock cycle, then half of the cores will be able to execute their instruction while the other half will wait for the next clock cycle, in effect dividing performance by two. This is known as the Single Instruction Multiple Threads (SIMT) architecture. The advantage of this architecture is that it allows the chip to contain more calculation units at the cost of the instruction fetching units. When GPU code is slowed down by cores in a warp wanting to do different instructions we speak of warp divergence. The most common cause of this is having "if" statements in the code as different threads will go down different branches.

Thread

A thread is an independent part of a program that can be run in parallel with other threads. A multi-threaded program is a program where the workload has been spread among different threads. These different threads can then be run in parallel on the different available physical cores in order to speed up the program. If a GPU has, say, 1000 cores, one needs at least 1000 threads to be running simultaneously in order to provide work for all the cores. However, a GPU is much more efficient if you can provide each core with several threads because GPU cores are designed to swap from thread to thread as a given thread may be waiting for a memory transfer to finish. This is why we need 15360 threads to get the best performance out of our Titan Black GPU which only has 2880 cores. Every thread in a block (the notion of block will be defined later on) has a unique index which determines which part of the computation it is responsible for.

Kernel

A kernel is a C++ function that can be executed on the GPU, it may call other functions. In order to exploit the GPU a CUDA code must declare at least one kernel and launch it.

Streaming Multiprocessor

A GPU is split into streaming multiprocessors each having an equal share of cores and cache memory. The notion of streaming multiprocessor is linked to the notion of block because a block of threads can only be run on a streaming multiprocessor.

Block

A kernel is always run on a certain number of blocks of threads. Threads in a block can access a "Shared Memory" cache and can synchronize among each other to make sure all the threads in a block have reached a chosen point in the computation before proceeding. Every block has a unique index which determines which part of the computation it is responsible for.

A block of threads can only be run on the resources of a streaming multiprocessor. If a kernel is run on less blocks than there are streaming multiprocessors on a GPU then part of the GPU will not have work to do. If a kernel is run on more blocks than there are streaming multiprocessors, the excess blocks may be run simultaneously or not depending on the limitations of the GPU in question. For example if running two blocks would require more shared memory than a streaming multiprocessor has, then the two blocks will be run sequentially and not in parallel. In any case it is usually a good idea to have a number of blocks that is a multiple of the number of streaming multiprocessors the GPU has because if your GPU has 15 streaming multiprocessors and you launch a kernel with 16 blocks the execution time may well be twice as long because one streaming multiprocessor gets double the work, the work isn't shared optimally for parallel execution.

In our case a Titan Black has 15 streaming multiprocessors. This is why in the ion sympathetic cooling simulations described in chapter 4 the optimal ion number is 15360 corresponding to 15 blocks of 1024 threads. Each streaming multiprocessor has 192 cores, making up the 2880 mentioned total.

Grid

For mathematical convenience when one launches a kernel on a certain number of blocks of threads one does not just specify the number but rather the layout in up to three dimensions. A kernel is launched on a grid of blocks of threads. The grid could have for example a width of 5 in the x dimension, 6 on the y dimension and 1 on the z dimension for a total of 30 blocks.

Global Memory

This is the most abundant memory on the GPU, typically counted in gigabytes but is also the slowest type of memory. The GDDR memory which composes the RAM in a GPU has a higher bandwidth but also a higher latency than the DDR RAM of a CPU. The latency of the global memory is of about 300 clock cycles.

Constant Memory

The GPU has 64 KB of constant memory. It is located in the global memory and therefore has the same latency but when something is fetched from constant memory the result is shared between threads in a warp. This memory is therefore useful to store constant numerical values that are regularly used during a calculation.

L1 Cache / Local Memory

The GPU has a relatively small amount of "local memory" or L1 cache. Every thread gets a share of local memory that will typically be used to store a few auxiliary values in a computation. If a thread runs out of local memory, local memory is "spilled" into global memory and this can have a performance cost. The order of magnitude of bytes of local memory a thread can have is 100.

Local memory is about 100 times faster than global memory.

Shared Memory

Shared memory is as fast as local memory but is shared between a block of threads. This can allow blocks of threads to efficiently work together on a problem.

Texture Memory

Memory accesses on a GPU are fastest when neighbouring threads access neighbouring bytes in memory. Texture memory allows for this efficiency to take place but in two dimensions which is very useful for computer graphics because pixels are on a plane. We do not use texture memory in our code.

Synchronicity and Asynchronocity

When the host gives a task to the device it may wait for the device to finish the task before proceeding with the rest of its code or not. If the CPU waits for the device to finish the task then the call is referred to as synchronous and if the CPU does not wait then the call is referred to as asynchronous.

Pinned memory

In CUDA we make an important distinction between pinned and non-pinned memory. Pinned memory is page-locked, meaning that it is necessarily in the RAM and therefore cannot be in the hard drive swap. Pinned memory accelerates transfers to and from the GPU greatly compared to its non-pinned counterpart.

B.3.2 Programming Syntax

Memory Allocation

Before being able to do anything we need to allocate some memory to work with. You can allocate page-locked memory on the host using cudaMallocHost and paged memory using new or malloc. You can allocate global memory on the device using cudaMalloc. Page-locked memory allows for much faster data transfers between the GPU and the CPU and allows the use of asynchronous memory transfers between host and device.

Stream

The purpose of streams is for the programmer to decide what will be run sequentially and in parallel. One launches kernels and data transfers in a stream. Everything that is launched in a stream takes place sequentially while the different streams may run in parallel. One first has to create a stream object and then assign it to a GPU. cudaStream t stream0 ; c u d a S e t D e v i c e ( 0 ) ; // Choose t h e f i r s t GPU a v a i l a b l e on t h e machine cudaStreamCreate(& stream0 ) ; One can ask the host process to wait for the completion of tasks started in a stream by calling cudaStreamSynchronize ( stream0 ) ;

Memory Transfers

In any CUDA program you have to transfer memory back and forth between the host and the device in order to give the GPU the data it should work on and to get the results back. This can be done using cudaMemcpy and cud-aMemcpyAsync. cudaMemcpy is a synchronous call and cudaMemcpyAsync is asynchronous. Furthermore, cudaMemcpyAsync only works if the host memory is page-locked.

You can also use cudaMemcpy to transfer memory from the CPU to the CPU and from the GPU to the GPU.

Kernel

A Kernel is declared just like a standard C++ function but with the global prefix. A kernel can't return a value so its type is always void. For example: This syntax allows allocating shared memory with a size known at compile time. Shared memory can also be allocated dynamically with a size which is only known at run time as a parameter of the kernel launch but we did not use dynamic shared memory allocation in our code.

Launching a Kernel

Launching a kernel is done with a triple chevron syntax as follows: MyKernel<<<gridDimension , blockDimension , s i z e , stream >>>( parameter1 , parameter2 , . . . ) ;

The size and stream parameters are optional. If no stream parameter is passed then the kernel is launched in the default stream which can be synchronised by calling cudaStreamSynchronize() without giving it a stream parameter. The size parameter allows dynamic allocation of shared memory.

The gridDimension and blockDimension parameters can be dim3 CUDA objects (object composed of 3 integers x,y,z) or scalar integers. Specifying a scalar integer is equivalent to specifying a dim3 object with 1 as the x and y dimensions.

Local memory

Local memory is declared implicitly as you would declare any variable in the scope of a kernel or a device function. In the following example d e v i c e void MyFunction ( ) { int MyInt ; . . . } MyInt will be stored in local memory unless there is no more of it and local memory is spilled into global memory.

Device and Host functions

In order for a function to be used on the GPU it must be declared with the prefix device . A function to be used on the CPU is declared with the prefix host but this is usually redundant because if you declare a function without either of these prefixes the compiler will assume the function is destined for the CPU. A function can be usable on both CPU and GPU by prefixing it with both host and device .

B.3.3 CUDA examples

In this section we will see a minimal working example of CUDA code. We will make a program that initialises the values of an array to 1,2,3... and recovers this array on the CPU. To do this on the CPU we would loop over the elements of the array. In CUDA code we can also use loops but a common paradigm is to, instead of looping, spread out work among threads. Therefore, instead of looping over 15360 array elements we can spawn 15360 threads each responsible for initialising the value of an array element. In this way, it is a common occurrence when porting CPU code to the GPU to have loops disappear. This minimalist version of the Coulomb computation can be optimised greatly by using shared memory to cache data and therefore reduce access to global memory which is slow. It is hardcoded for 15360 ions and targets a GPU with 15 streaming multiprocessors. If this code was deployed on a GPU with 14 multiprocessors the run time would be twice as long as 1 streaming multiprocessor would have to process two blocks of threads.

Fig. 0. 2 :

 2 Fig. 0.2: Illustration of our spectroscopy experiment. H + 2 in state 1sσ g |ν = 0, L = 2 is excited to 1sσ g |ν = 1, L = 2 by a 9.17 µm laser and then dissociated by a 213 nm laser. The purpose of the dissociation is to determine whether or not the excitation to the 1sσ g |ν = 1, L = 2 level occurred. By sweeping the frequency of the excitation laser we expect maximal loss of ions at resonance.

Fig. 0. 3 :

 3 Fig. 0.3: Schematic of the GBAR experiment's different steps to produce H+ , cool it and measure its free fall.

Fig. 1 . 1 :

 11 Fig. 1.1: Linear Paul trap built by our team [31]. The electrodes are made of molybdenum held together by vacuum compatible ceramics for electric insulation. Diametrically opposed electrodes are separated by 2r 0 = 7 mm. The small rod electrodes are 12 mm long.

Fig. 1 . 2 :

 12 Fig. 1.2: Photo of the hyperbolic Paul trap from our team. Stainless steel hyperbolic trap made of a ring electrode with an hyperbolic shape and two endcaps (the top one is missing). The white ceramics ensure electric insulation between the different electrodes. Note the holes in the ring electrode used for molecular beam intake and UV irradiation.

Fig. 1 . 3 :

 13 Fig. 1.3: Representation of the electrodes of a hyperbolic Paul trap. The ring electrode is at potential V = U 0 + V 0 cos (Ωt) while the endcap electrodes are grounded. The curvature of the electrodes is hyperbolic.

Fig. 1 . 4 :

 14 Fig. 1.4: Left: Stability diagram of a hyperbolic Paul trap. The red and blue areas respectively represent the stability regions for the motion along the axial (z) and radial (x, y) directions. 3D trapping is therefore achieved in the intersection of these two areas. Right: Stability diagram of a linear Paul trap.Here the red and blue areas correspond to the x and y directions respectively. By courtesy of J. Heinrich[START_REF] Heinrich | A Be + Ion Trap for H + 2 Spectroscopy[END_REF].

Fig. 1 . 5 :

 15 Fig. 1.5: Schematic representation of a linear Paul trap. A set of 4 central electrodes confine the ions radially while another set of 4 end cap electrodes provide the axial confinement with a voltage U end ≥ 0. The electrode radius r e ≈1.15r 0 because this geometry was found to minimize anharmonicities of the potential[START_REF] Reuben | Ion trajectories in exactly determined quadrupole fields[END_REF]. By courtesy of J. Heinrich[START_REF] Heinrich | A Be + Ion Trap for H + 2 Spectroscopy[END_REF].

Fig. 1 .

 1 4 on the right we show the first stability region of the linear Paul trap which is what most experiments use. If |a x |, |a y |, |q x |, |q y |1 the adiabatic approximation is valid and[START_REF] Ghosh | Ion Traps[END_REF]

Fig. 1 . 6 :

 16 Fig. 1.6: False colour photo of a Be + Coulomb crystal taken in our team's linear Paul trap.

Fig. 1 . 7 :

 17 Fig. 1.7: Value of the potential, for an ion of charge 1 e, at the edge of a spherical crystal composed of ions of charge 1 e as a function of the number of ions, with ωz = 2π • 250 kHz using Eq. 1.36.

Fig. 2 . 1 :

 21 Fig. 2.1: Illustration of the principle of REMPI. A molecule in its ground electronic state is excited to a Rydberg state resembling an ionic core in a desired ro-vibrational state plus a loosely bound electron. This electron is then expelled with an extra photon and the ionic core in the desired ro-vibrational state is left.

Fig. 2 . 2 :

 22 Fig. 2.2: Left: Copper skimmer with an opening of 150 µm mounted onto a CF-63 to CF-40 adaptor flange. Right: Schematic representation of the skimmer.

Fig. 2 . 3 :

 23 Fig. 2.3: ACPV2 pulsed valve (made by MassSpecpecD BV) mounted onto a CF-63 flange with gas and BNC feedthroughs. The diameter of the valve is 150 µm, pulse duration is between 5 µs and 200 µs and the valve can be operated at up to 5 kHz.

Fig. 2 . 4 :

 24 Fig. 2.4: Schematic of the H + 2 ion source design. Pulses of H 2 gas pass through a skimmer to be photo-ionized by the REMPI laser and state-selectively produce H + 2 in the trap.

Fig. 2 . 5 :

 25 Fig. 2.5:The molecular beam passes through an approximate area considered to be the focus of the REMPI laser beam where 3+1 photo-ionization will have a significant probability of occurring.

Fig. 2 . 6 :

 26 Fig. 2.6: Mean energy per pulse of our REMPI laser at 303 nm versus the Q-switch flash lamp delay in the Nd:YAG pumping the dye laser. The longer the delay between the flash lamps and the Q-switch the more time the population inversion has to relax due to spontaneous emission. 190 µs is the minimum allowed by the laser.

Fig. 2 . 7 :

 27 Fig. 2.7: Energy curve of the dye laser, before frequency doubling, as we added millilitres of Rhodamine-640 solution to the previously used dye mixture. The solution of the "+11" curve (0.088 g L -1 of Rhodamine-610 and 0.032 g L -1 of Rhodamine-640) is what we used in our experiments until we broke a dye cell. Afterwards we used a more dilute solution (concentrations given in text).

Fig. 2 . 8 :

 28 Fig. 2.8: Ion source design. A pulsed valve is mounted onto a flange with gas and electrical feedthroughs. The pulsed valve is aimed towards a skimmer which only lets through the part of the beam aligned with the centre of the trap (unrepresented).

Fig. 2 . 9 :

 29 Fig. 2.9: Photo of the experimental setup of the ion source. To the left of the flexible bellow is the hyperbolic trap (not shown).

Fig. 2 . 10 :

 210 Fig. 2.10: Trapped ion signal (in arbitrary units) as the wavelength of the REMPI laser is swept. The peaks can be identified with those in Fig. 1 of Ref. [65]with a small shift due to imperfect calibration of our laser. There was 1 bar in the H 2 reservoir, the 313 nm was not at its best, at around 3.3 mJ, the gas pulses were 200 µs long and the creation time was 2 s.

Fig. 3 . 1 :

 31 Fig. 3.1: Dissociation cross section for the different vibrational states of H + 2 versus the wavelength of the dissociation laser. Calculations done in Ref. [70].

Fig. 3 . 2 :

 32 Fig. 3.2: Illustration of the r, ρ and z coordinates. The solid circle represents the spherical ion cloud the dissociation laser beam going through it along the z axis.

3 2 ∆r 3 e -( r ∆r ) 2 ( 3 . 7 )Figure 3 .

 332373 Figure 3.2 shows the definition of the spherical (r), cylindrical (ρ) and axial (z) coordinates and illustrates the relation r 2 = ρ 2 + z 2 .By combining Eq. 3.6 and Eq. 3.7 we can write the remaining fraction of ions in the cloud after one laser pulse as

Fig. 3 . 3 :

 33 Fig. 3.3: Evolution of the pressure in the chamber after H + 2 electron-impact creation (5 s with a current I=1.04 A applied to the electron gun), fitted by a simple (green) and double (blue) exponential decay law. The green fit is of the form c 1 e -t β + c 2 with β = 3.19 ± 0.03 s, c 1 = 2.78 × 10 -9 ± 3 × 10 -11 torr and c 2 = 2.76 × 10 -10 ± 1 × 10 -12 torr. The blue fit is of the form c 1 e -t t 1 + c 2 + c 3 e -t t 2 with t 1 = 23 ± 1.5 s, t 2 = 2.43 ± 0.03 s, c 1 = 2.13 × 10 -10 ± 9 × 10 -12 torr, c 2 = 2.43 × 10 -10 ± 2 × 10 -12 torr and c 3 = 3.11 × 10 -9 ± 3 × 10 -11 torr.

Fig. 3 . 4 :

 34 Fig. 3.4: Zoom on the beginning of Fig. 3.3 up to 10 s (right). The solid line is a fit of the form c 1 + e -t β + c 2 with β = 2.5 ± 0.1 s, c 1 = 3.2 × 10 -9 ± 1 × 10 -10 torr and c 2 = 3.8 × 10 -10 ± 2 × 10 -11 torr.

3. H + 2 48 Fig. 3 . 5 :

 24835 Fig. 3.5: Pressure increase when 13.5 mW of the 213 nm laser is applied without interruption. Left: full record fitted with an equation of the form c 1 e -t c 2 + c 3 with c 1 = -2.43 × 10 -10 ± 2.5 × 10 -12 torr, c 2 = 7.3 ± 0.1 s and c 3 = 4.319 × 10 -10 ± 5 × 10 -13 torr. Right: short term behaviour fitted with an equation of the form c 1 e -t c 2 + c 3 with c 1 = -2.46 × 10 -10 ± 6 × 10 -12 torr, c 2 = 5.2 ± 0.35 s and c 3 = 4.12 × 10 -10 ± 7 × 10 -12 torr.

Fig. 3 . 6 :

 36 Fig. 3.6: Limit pressure in the trap chamber after continuous exposition to different powers of the 213 nm laser dissociation laser. The solid line is a fit of the form ηP total + c 1 with η = 3.6 × 10 -12 ± 1 × 10 -13 torr mW -1 and c 1 = 1 × 10 -10 ± 1 × 10 -11 torr.

Fig. 3 . 7 :

 37 Fig. 3.7: After 2 s of creation by electron gun, H + 2 ions are extracted for detection by the MCP. Ion signal versus extraction time is plotted and fitted (solid line) with Eq. 3.22 with A = 32 ± 22, B = 0.13 ± 0.06 Hz, C = 0.6 ± 0.7 and β = 2.5 ± 2.9 s. The 213 nm dissociation laser is off.

Fig. 3 . 8 :

 38 Fig. 3.8: Remaining fraction of ions created by electron impact as a function of time for different attenuations of the 213 nm laser. Each data point is the average of 16 measurements. The error bar is the standard error on the average, i.e. the standard deviation of the measurements divided by 4. The solid lines are fits according to Eq. 3.23 with κ = 0.003 98 mW -1 , ∆r = 1.13 mm, the powers are given in Tab. 3.2 and we assumed β = 0 s.

Fig. 4 . 1 :

 41 Fig. 4.1: An ion at point M is projected to point H on the beam axis. O is the origin of our system of coordinates and C is the position of the waist.

(4. 18 )Fig. 4 . 2 :

 1842 Fig. 4.2: Left: Ion cloud initialised randomly in a cylinder. Right: Ions in their approximate equilibrium position after the cylinder is simulated for 0.1 ms with a fictitious friction force and no laser interaction.

Fig. 4 . 3 :

 43 Fig. 4.3: Macromotion kinetic energy per ion as a function of time with the fictitious friction force replaced by the proper laser interaction at around 0.1 ms. We see that under the fictitious friction force the energy of the 1919 Be + in purple goes below the Doppler limit. At around 0.1 ms the fictitious friction force is replaced by the proper laser interaction and the cloud immediately gains two to three orders of magnitude in energy. With time the temperature stabilises close to the Doppler limit at a few mK. The laser is aligned along the z axis, detuned by -Γ with Γ = 19.4 MHz the natural width of the cooling transition and laser irradiance is at 1.5 times saturation with a waist of 1 mm. Ω = 13 MHz, r 0 = 3.5 mm, a x = 5.2 × 10 -7 and q x = 0.011 for the Be + and ωz = 125 kHz.

Fig. 4 . 4 :

 44 Fig. 4.4: Coulomb crystal composed of 1469 Be + and 450 H + 2 with a x = 5.25 × 10 -7, q x = 0.0223 for the Be + and ωz = 250 kHz. Characteristic lengths are given along with linear density n. By courtesy of J. Heinrich[START_REF] Heinrich | A Be + Ion Trap for H + 2 Spectroscopy[END_REF].

Fig. 4 . 5 :

 45 Fig. 4.5: A 2.5 meVH+ is sent into a crystal of 1919 Be + . On top we show the trajectory of the H+ and on the bottom we show the averaged (right) and unaveraged (left) macromotion kinetic energy of Be + in purple and the H+ in green. The energy given is per ion which explains the small bumps in the purple curve (1919 ions) while the green one has bigger ones (1 ion). a x = 4.7 × 10 -6 , q x = 0.052 for the H+ , ωz = 65 kHz and other parameters as in Fig.4.3. LC stands for the laser-cooled species, the beryllium in our case. SC stands for the sympathetically cooled species, H+ in this case.

Fig. 4 . 6 :

 46 Fig. 4.6: In a) the timestep is short enough to describe the collision of two ions. In b) the timestep is too long and the two ions go through each other.

Fig. 4 . 7 :

 47 Fig. 4.7: Energy of a projectile ion after a collision with a stationary target ion of equal mass as a function of simulation timestep for different energies of the projectile. Arrows mark the critical timestep below which the collision is approximately well described.

Fig. 4 . 10 :

 410 Fig. 4.10: The situation is identical to Fig. 4.9 except that the ion has more radial motion energy and we look at the conservation of that energy by looking at the amplitude of motion along x. Other parameters as in Tab. 4.2.

Fig. 4 . 11 :

 411 Fig. 4.11: Test of RF heating on a cloud of 1920 Be + for different values of the β parameter in the timestep requirement of Eq. 4.19 (the corresponding values of 1 β are given in the figure). The laser cooling is turned off and we see the average energy of the ions go up slowly with time. The crystal starts at about 1 K in order to have much more noticeable RF heating than we would at the Doppler limit because of RF heating's quadratic temperature dependence. Other parameters as in Tab. 4.2.

Fig. 5 . 1 :

 51 Fig. 5.1: Capture time t of a 62.5 meV mass 8 ion as a function of the ellipticity ωx ωz of a spheroid Be + crystal. An ellipticity of 1 corresponds to a sphere, superior to 1 corresponds to a prolate crystal and inferior to 1 corresponds to an oblate (pancake-shaped) crystal. Ω = 13 MHz, r 0 = 3.5 mm. Data and other parameters in Tab. 5.1.

Fig. 5 . 2 :

 52 Fig. 5.2: Capture time of an H+ by a crystal of laser-cooled Be + as a function of the energy of the H+ . The error bars are taken to be the standard deviation divided by the square root of the number of simulations. At 3.125 meV cooling is not always successful and only data from successful tries is used. The solid line is a fit of the form c 1 • E 2 + c 2 with c 1 = 0.86 ± 0.05 ms meV -2 and c 2 = -0.29 ± 0.02 ms. Eq. 5.5 predicts c 1 = 0.73 ms meV -2 assuming ln Λ = 1. The dashed line is an identical fit that does not weight error bars with c 1 = 1.4 ± 0.2 ms meV -2 and c 2 = -0.8 ± 1.3 ms. We comment on the fit in Sec. 5.9. Ω = 13 MHz, r 0 = 3.5 mm, a x = 4.7 × 10 -6 , q x = 0.052 for the H+ , ωz = 65 kHz. The laser is aligned along the z axis, detuned by -Γ with Γ = 19.4 MHz the natural width of the cooling transition, and laser irradiance is at 1.5 times saturation with a waist of 1 mm.

Fig. 5 . 3 :

 53 Fig. 5.3: Energy of an H+ starting at 3.125 meV sympathetically cooled by a crystal of 1919 Be + . After 120 ms the H+ has actually heated up to about 8.6 meV while the Be + has also heated up by a factor of about 2. Other parameters as in Fig. 4.5.

Fig. 5 . 4 :

 54 Fig. 5.4: Energy of a 21 meVH+ in green as it gets captured by a crystal of 1919 Be + in purple under the pseudopotential approximation. a x = 0.00047, q x = 0.31 for the H+ , ωz and other parameters as in Fig.5.2.

Fig. 5 . 5 :

 55 Fig. 5.5: Average capture time as a function of initial energy for an H+ sent into a crystal of 1441 Be + and 350 HD + with stability parameters q x = 0.2 and a x = 0.0014 for the H+ , ωz = 250 kHz and other parameters as in Fig. 5.2. The error bars are the standard deviation divided by the square root of the number of simulations. The solid line is a fit of the form c 1 x 2 + c 2 with c 1 = 0.015 ± 0.004 ms meV -2 and c 2 = -2 ± 1.5 ms. Eq. 5.5 predicts c 1 = 0.046 ms meV -2 with ln Λ = 1 and taking the mass, density and charge of the auxiliary ions. The dashed line is an identical fit that does not weight error bars with c 1 = 0.0274 ± 0.0003 ms meV -2 and c 2 = -11 ± 2 ms. We comment on the fit in Sec. 5.9.

Fig. 5 . 6 :

 56 Fig. 5.6: Energy of a 125 meV H+ in green being sympathetically cooled by a crystal of 1441 Be + in purple and 350 HD + in blue with parameters as in Fig. 5.5. "INT" refers to the auxiliary ions with an intermediate mass in between that of the laser-cooled and the sympathetically-cooled species used to improve mechanical coupling.

Fig. 5 . 7 :

 57 Fig. 5.7: Left: Image of the asymmetric Be + , HD + crystal which was used in the simulations of Fig. 5.5. Right: Image of the asymmetric Be + , H + 2 crystal which was used in the simulations of Fig. 5.11. Be + is in blue, H + 2 in black and HD + in green.

Fig. 5 . 8 :

 58 Fig. 5.8: An H+ ion with an initial energy of 31.25 meV is sympathetically cooled by a crystal of 1441 Be + and 350 HD + but it gets into a stable orbit around the crystal. Top: Trajectory of the H+ (r = x 2 + y 2 ). Bottom:

Fig. 5 . 9 :

 59 Fig. 5.9: Energy of a 125 meV H+ in green being sympathetically cooled by a crystal of 991 Be + in purple and 800 HD + in blue. The purple curve going to 0 around 500 ms corresponds to loss of all the Be + in the trap. Other parameters as in Fig. 5.5.

Fig. 5 . 10 :

 510 Fig. 5.10: Energy of a 62.5 meV H+ in green being sympathetically cooled by a crystal of 1469 Be + in purple and 450 H + 2 in blue. Other parameters as in Fig. 5.5.

Fig. 5 . 11 :

 511 Fig. 5.11: Average capture time as a function of initial energy for anH+ sent into a crystal of 1469 Be + and 450 H + 2 with other parameters as in Fig.5.5. At 62.5 meV the cooling does not always work and the capture time given is for the cases when it succeeded. The error bars are the standard deviation divided by the square root of the number of simulations. The solid line is a fit of the form c 1 x 2 + c 2 with c 1 = 0.0069 ms meV -2 ± 0.0002 ms meV -2 and c 2 = -1.16 ± 0.07 ms. Eq. 5.5 predicts c 1 = 0.042 ms meV -2 with ln Λ = 1 and optimistically taking the mass, density and charge of the auxiliary ions. We comment the fit in Sec. 5.9.

Fig. 5 . 15 :

 515 Fig. 5.15: The data of Tab. 5.6 (purple) and the data of Tab. 5.4 are fitted with a c 1 q c 2x law. For the pseudopotential data we find c 1 = 0.005 ± 0.002 ms and c 2 = -1.8 ± 0.1, and for the exact potential data we fixed c 2 = -7 3 and find c 1 = 6.9 × 10 -4 ± 8 × 10 -5 ms. We could not fit while weighting the error bars because of the undefined error bars on some data points (see caption of Tab. 5.6). Eq. 5.12 predicts 3.17 × 10 -3 ms with ln Λ = 1.

Fig. 5 . 16 :

 516 Fig. 5.16: Average capture time of mass 8 ions in a crystal of 1919 Be + as a function of the initial energy with a q x = 0.25 and a x = 3 × 10 -5 for the mass 8 ions, ωz = 2.51 MHz and other parameters as in Fig. 5.2. The error bars are taken to be the standard deviation divided by the square root of the number of simulations. The solid line is a fit of the form c 1 E 2 + c 2 with c 1 = 9 × 10 -5 ± 1 × 10 -5 ms meV -2 and c 2 = -4 × 10 -3 ± 7 × 10 -3 ms. Eq. 5.5 predicts c 1 = 2.3 × 10 -4 ms meV -2 with ln Λ = 1. The dashed line is fit by the same formula without weighting error bars with c 1 = 1.25 × 10 -4 ± 7 × 10 -6 ms meV -2 and c 2 = 1 × 10 -4 ± 1.5 ms.

Fig. 5 . 17 :

 517 Fig. 5.17: Fluorescence photons(top) and its Fourier transform(bottom) from Be + in a crystal of 1441 Be + and 350 HD + as it cools down a H+ ion with an energy of 125 meV. Fluorescence points are 1 × 10 -7 s apart. Trapping parameters are identical to Fig. 5.5.

Fig. 5 . 18 :

 518 Fig. 5.18: Height of the 30 kHz peak as a function of time in the absence (left) and presence (right) of the H+ .Both plots are given in units of the baseline height of the peak when the crystal is in equilibrium. For a short time the peak is higher than the baseline on the left-hand plot because of the non-stationary state the crystal is in before reaching its equilibrium. Parameters are identical to Fig.5.17 with the data of the right-hand plot being from that figure.

Fig. 5 . 19 :

 519 Fig. 5.19: Capture time of a 125 meV H+ in crystals composed of different ion numbers with approximately 20% HD + and 80% Be + . The full line is a fit of the form c 1 N c 2Ions with the fitting parameters found to be c 1 = 1.3 × 10 7 ±9 × 10 6 ms and c 2 = -1.38±0.09. ωz = 250 kHz, q x = 0.2 for the H+ and other parameters described in Tab. 5.8.

Fig. 5 . 20 :

 520 Fig. 5.20: One of the simulations for a 125 meV H+ in a crystal of 2883 Be + and 700 HD + , in which heating occurs instead of sympathetic cooling. Parameters given in Tab. 5.8. As the heating occurred Be + ions started heating up and touching the electrodes (therefore being removed from the simulation), wich aggravated the situation.

  g l o b a l void MyKernel ( parameter1 , parameter2 , . . . ) { . . . } Shared Memory Shared memory can be declared inside a kernel or a device function with the following syntax: s h a r e d int B u f f e r [ 1 0 2 4 ] ;

  g l o b a l void I n i t A r r a y ( int * d t ) { // D e f i n e t h e k e r n e l int i n d e x=b l o c k I d x . y * blockDim . x+t h r e a d I d x . x ; d t [ i n d e x ]= i n d e x ; } int main ( ) { int * h t , * d t ; cudaMemcpy ( d i o n , h i o n , 1 5 3 6 0 * s i z e o f ( i o n ) , cudaMemcpyHostToDevice ) ; Coulomb<<<15,1024>>>( d i o n , d f ) ; c u d a D e v i c e S y n c h r o n i z e ( ) ; // Get f o r c e s back from GPU, p r i n t them , c h e c k them . . . }

  Trapping parameters and secular frequencies for different ionic species of interest for us with masses M, in atomic mass units, and charge Q, in units of the charge of the electron, in a linear Paul trap with U 0

		25 × 10 -4 0.0525 239 kHz	229 kHz	100 kHz
	HD +	3 1 1.57 × 10 -3 0.157	718 kHz	708 kHz	173 kHz
	H + 2	2 1 2.36 × 10 -3 0.236	1.08 MHz 1.07 MHz 212 kHz
	H+	1 1 4.72 × 10 -3 0.472	2.17 MHz 2.16 MHz 300 kHz
	Tab. 1.2:

  Average and standard deviation of the capture time for a 625 meV mass 8 ion in a crystal of 1919 Be + as a function of the α parameter in the variable timestep criterion described in Eq. 4.26 and the β parameter of Sec. 4.6.1. a x = 5 × 10 -5

	1 α	1 β	Mean Capture Time Standard Deviation N Curves
	5	15	13 ms	4 ms	33
	15	30	11 ms	5 ms	32
	5	60	14 ms	4 ms	32
	30	60	16 ms	5 ms	30
	60	60	13 ms	4 ms	26
	1000 300 13 ms	4 ms	14
	Tab. 4.2:

  The timestep cannot be longer than 1.28 × 10 -9 s because of the timestep constraint on describing the RF potential discussed in Sec. 4.6.1. Simulations done with α = 5 in Eq. 4.26. Other parameters as in Fig.4.5.

	Ion Number With optimisation Without optimisation
	15360	1	1
	7680	0.95	0.68
	3840	0.91	0.38
	1920	0.73	0.2
	960	0.41	0.09
	480	0.15	0.04
	240	0.04	0.03
	Tab. 4.5: Performance of Coulomb interaction calculation with and without the low
	ion number performance optimisation normalised to the performance of
	the 15360 ion case.
			[s] Average Timestep
	0.625	1.28 × 10 -9	1.28 × 10 -9 s
	1.25	1.28 × 10 -9	"
	1.875	4 × 10 -10	"
	2.5	2 × 10 -10	"
	3.15	1 × 10 -10	"
	6.25	8 × 10 -11	"
	31.5	2.4 × 10 -11	"
	62.5	1 × 10 -11	"
	312.5	5.3 × 10 -12	"
	625	3 × 10 -12	"
	1250	1.4 × 10 -12	"

Tab. 4.4: Minimum and average timestep required to simulate sympathetic cooling of an H+ in a crystal of Be + as a function of the energy of the H+ .

  .2.

	q x	ωz	t capt	STDEV N Simulations
	0.1	250 kHz 14.3 ms 16.5 ms 69
	0.15 385 kHz 7.62 ms 9.57 ms 116
	0.2	500 kHz 2.69 ms 1.55 ms 40
	0.25 625 kHz 2.18 ms 1.97 ms 46
	0.3	750 kHz 2.36 ms 2.54 ms 119
	0.4	1 MHz	3.92 ms 9.44 ms 29
	Tab. 5.3: Capture time of a 21.25 meV H + 2

  [START_REF] Heiße | High-precision measurement of the proton's atomic mass[END_REF].6 (Be + /HD + / H+ ), Fig.5.11 (Be + /H + 2 / H+ ) and Fig.5.16 (Be + /mass 8) we fit capture times as a function

	q x	ωz	t capt	STDEV N Simulations
	0.025	250 kHz	7.65 ms	3.24 ms	17
	0.0375 385 kHz	3.19 ms	2.12 ms	41
	0.05	500 kHz	1.41 ms	0.496 ms 180
	0.0625 625 kHz	0.759 ms 0.216 ms 154
	0.075	750 kHz	0.467 ms 0.105 ms 595
	0.1	1 MHz	0.268 ms 0.073 ms 52
	0.125	1.25 MHz 0.164 ms 0.043 ms 310
	0.15	1.5 MHz	0.1 ms	0.024 ms 512
	0.175	1.75 MHz 78.6 µs	22.4 µs	154
	0.2	2 MHz	58.7 µs	16.8 µs	155
	0.225	2.25 MHz 45.4 µs	12.7 µs	592
	0.25	2.51 MHz 37.5 µs	11.2 µs	578
	0.275	2.75 MHz 29.5 µs	9.03 µs	155
	0.3	3 MHz	22.5 µs	6.73 µs	152
	0.375	3.75 MHz 15.4 µs	4.74 µs	65
	0.5	5 MHz	8.47 µs	3.21 µs	51
	0.625	6.25 MHz 7.8 µs	2.83 µs	1116

  .2.

	q x	ωz	t capt	STDEV N Simulations
	0.05	135 kHz 35.6 ms 11.1 ms 34
	0.1	250 kHz 6.55 ms 1.07 ms 24
	0.125	310 kHz 4.02 ms 1.63 ms 40
	0.15	385 kHz 5.33 ms 7.5 ms	29
	0.155	390 kHz 4.34 ms 8 ms	81
	0.1625 400 kHz 15.6 ms 10.5 ms 42
	0.175	440 kHz 26 ms	11.4 ms 25
	0.2	500 kHz 29.4 ms 8.24 ms 5
	0.25	625 kHz 24.9 ms 3.42 ms 5
	0.3	750 kHz 13.7 ms 1.79 ms 9

  .2.

	q x	a x	ωz	t capt	STDEV N Simulations
	0.0125 3 × 10 -5 135 kHz	23.8 ms	3.96 ms	52
	0.025	3 × 10 -5 250 kHz	3.85 ms	0.72 ms	36
	0.0375 3 × 10 -5 385 kHz	2.29 ms	1.28 ms	58
	0.05	3 × 10 -5 500 kHz	1.94 ms	0.743 ms 70
	0.0625 3 × 10 -5 625 kHz	1.58 ms	0.286 ms 80
	0.075	3 × 10 -5 750 kHz	0.621 ms 0.130 ms 174
	0.2	3 × 10 -6 2 MHz	<60 µs		40
	0.25	6 × 10 -6 2.51 MHz <60 µs		74
	0.375	3 × 10 -6 3.75 MHz <60 µs		71
	0.5	6 × 10 -6 5 MHz	<60 µs		44

  .2.

	q x	a x	ωz	t capt	STDEV N Simulations
	0.15	8.3 × 10 -4 190 kHz 17.7 ms 16.8 ms 5
	0.175 1.1 × 10 -3 215 kHz 16.6 ms 7.88 ms 24
	0.2	1.4 × 10 -3 250 kHz 10.1 ms 6.12 ms 40
	0.225 1.7 × 10 -3 280 kHz 17.4 ms 27.6 ms 21
	0.25	1.9 × 10 -3 300 kHz 14.8 ms 26.8 ms 133
	0.3	1.4 × 10 -3 375 kHz 17.3 ms 19.7 ms 9

Photodissociation

Eq. 4.23 can also be used as a variable timestep criterion by using as a timestep the minimum of all the δt found for all ion pairs, this criterion is discussed in some more detail in our publication[START_REF] Sillitoe | Numerical Simulations Of Ion Cloud Dynamics[END_REF]

ACKNOWLEDGEMENTS

I would like to thank Laurent Hilico and Jean-Philippe Karr for supervising my PhD work and for their guidance in the writing of the manuscript, as well as my fellow PhD candidates Johannes Heinrich, who started at the same time as I did, and Thomas Louvradoux, who started two years later, for their involvement in our research. I would also like to thank the rest of the researchers and the technical and administrative personnel of the LKB. Looking back on my PhD, what is most memorable to me is that being surrounded by hard-working, knowledgeable and passionate people was inspiring.

This work was financed by the ANR BESCOOL 13-ISO4-002-01 contract.

CONTENTS

B.2.5 Passing by reference

When you want to pass parameters to a function you can pass either by copy (also known as passing by value) or by reference. For example

. . . } is a function which returns an integer and takes an integer by copy as input. The function

. . . } takes that same integer by reference as input. There are two consequences to this, one is that copy operations can be avoided, which can be expensive for large objects and the other is that a variable passed by reference can be modified by the function outside of its scope. Whereas the function which was passed the integer by copy can only modify the local value of that integer. For example int i =2; f u n c t i o n c o p y ( i ) ; // i s t i l l has a v a l u e o f 2 f u n c t i o n r e f ( i ) ; // i now has a v a l u e o f 3

Passing by reference avoids an often unnecessary copy operation. But sometimes passing by copy is the desired behaviour and passing by reference of small basic data types like doubles is unnecessary because the size of the double is the same size as the reference to the double anyway. It is recommended to pass basic data types by copy and larger objects by reference, unless copy is the desired behaviour.

B.2.6 Object Oriented Programming

Classes allow the programmer to make custom data types to simplify his code. For example instead of having 3 integer x,y and z variables one can define a vec class containing the 3 integers. The code will compile to the same executable but the objects can help simplify the code and avoid bugs. For example instead of doing an operation on x,y and z on three slightly different lines one can define custom vec operators. For example to avoid copy pasting errors when writing repetitions of three lines and to write things as standard vectorial operations as we would write them in mathematical equations.

B.2.7 Main

In C++ a program needs to have a function called main with an integer return type. Main is executed at the start of the program and does whatever you want it to, including calling your other functions. Main returns an integer to the operating system when its execution is over to signal the proper termination of the program.

B.3 Introduction to CUDA

In this section I hope to introduce the novice reader to CUDA programming with a lot of the terminology and concepts, leading up to a minimal working example of CUDA code.

B.3.1 Terminology

In this section I introduce the terminology used in CUDA. A lot of these concepts are present in OpenCL merely with a different name, as the underlying hardware architecture is the same.

Host and Device

The GPU is also sometimes referred to as an accelerator and works under instructions from the CPU in a master-slave relationship. The CPU is referred to as the host and the GPU is referred to as the device. cudaMallocHost (& h t , s i z e o f ( int ) * 15360) ; // s i z e o f i s a s t a n d a r d c++ f u n c t i o n cudaMalloc (& d t , s i z e o f ( int ) * 15360) ; I n i t A r r a y <<<15,1024>>>( d t ) ; // Launch t h e k e r n e l // Wait f o r t h e k e r n e l we l a u n c h e d i n t h e d e f a u l t stream t o f i n i s h c u d a D e v i c e S y n c h r o n i z e ( ) ; // Get t h e i n i t i a l i s e d a r r a y back from t h e GPU cudaMemcpy ( h t , d t , s i z e o f ( int ) * 15360 , cudaMemcpyDeviceToHost ) ; } In the code above we first declare two integer pointers d h and d t. We then allocate memory using cudaMallocHost and cudaMalloc which modify the value of the pointers to point towards the newly allocated memory. Next we declare our kernel Init Array, in which every thread of every block in the grid is responsible for initialising a different element of the array based on its index in its block and the index of its block in the grid. We then launch the kernel, wait for the end of its execution and transfer the array from GPU memory to CPU memory.

Note that in the naming of h t and d t we have used the standard practice in CUDA of distinguishing host and device memory with the prefixes d and h. Now that we have given a minimal working example of CUDA code we give a few extra examples. The following is a kernel which multiplies the previously created array by a constant k passed as a parameter.

Abstract

The high-resolution ro-vibrational spectroscopy of the H + 2 molecular ion by resonance-enhanced multiphoton dissociation (REMPD) requires a state-selective source of H + 2 . In this thesis we present work on a functional state selective H + 2 ion source using resonance enhanced multiphoton ionisation (REMPI) with a 303 nm pulsed laser. The second part of the thesis presents numerical simulations of sympathetic cooling in linear RF traps, whose main application is the GBAR project (Gravitational Behaviour of Antihydrogen at Rest) which involves sympathetic cooling of an antimatter H+ ion by laser-cooled Be + ions. We developed a GPU code using a variable timestep technique enabling a fast description of Coulomb interactions. We discuss the influence of RF heating and scaling laws between cooling times, initial energy and ion numbers in the cooling crystal. We show that the H+ sympathetic cooling step of GBAR could be feasible using a rotationally asymmetric two-component Be + /HD + crystal which appears more effective than a single-component Be + crystal. We find that the H+ ion's capture by this crystal could be detected experimentally by Fourier analysis of the fluorescence data.

Résumé

La spectroscopie ro-vibrationelle de haute résolution de l'ion moleculaire H + 2 par REMPD requiert une source sélective en état interne. Dans cette thèse nous présentons notre travail de conception et de réalisation d'une telle source utilisant la photo-ionisation multiphotonique résonante (REMPI) de H 2 à l'aide d'un laser pulsé à 303 nm. Dans un deuxième temps nous présentons nos simulations numériques de refroidissement sympathique dans un piège de Paul linéaire, avec pour application principale le projet GBAR qui implique le refroidissement sympathique d'un ion d'antimatière H+ par des ions Be + refroidis par laser. Nous avons dévelopé un code GPU utilisant un pas de temps variable permettant de décrire les interactions coulombiennes de façon efficace. Nous discutons de l'influence du chauffage RF et de lois d'échelles entre le temps de capture, l'énergie initiale et le nombre d'ions dans le cristal. Nous montrons que le refroidissement sympathique de H+ requis pour GBAR pourrait fonctionner avec un cristal dissymétrique de Be + /HD + qui semble plus efficace que le Be + seul. Nous montrons qu'avec un tel cristal la capture du H+ pourrait être détectée expérimentalement par analyse de Fourier des données de fluorescence.