Analysis methods and code for very high-precision mass measurements of unstable isotopes
Jonas Karthein, Dinko Atanasov, Klaus Blaum, David Lunney, Vladimir Manea, Maxime Mougeot
IINTRODUCTION — One of the fundamental ob-servables of the atomic nucleus is the nuclear binding energy. This quantity directly results from the difference of the nuclear mass and the sum of the individual nuclear constituent masses (neutrons and protons). While measurements of stable systems allow for long preparation and measurement times, the precision for short-lived nuclides is strongly limited by the production rate at the radioactive ion beam (RIB) facilities (down to a few atoms per hour), as well as the half-lives (down to a few tens of milliseconds). However, studying these exotic systems often provides exclusive insights into the nuclear struc-ture, the weak interaction or rapid nucleosynthe-sis processes compared to stable atoms. This particularity motivates technical advances in the field of high-precision mass spectrometry. Over the past four decades, Penning traps have emerged as the experimental tool of choice for high-precision measurements of atomic masses [Bl06, Hu17]. The ISOLTRAP spectrometer has pio-neered the technique for radioactive species [Kl13] and Penning trap facilities now exist at all radioactive ion beam facilities worldwide [Lu18]. Further developments in the measurement tech-nique using ISOLTRAP [Bl02, Ge07] allowed us to recently set a new world record for the highest precision on a mass of a short-lived (T ½ < 1 h) iso-tope [Ka19a]. But this experiment also demon-strated the technique’s intrinsic precision limita-tion, requiring hundreds of ions per measure-ment cycle, which is not feasible for the most short-lived isotopes of interest for nuclear physics. ABSTRACT — We present a robust analysis code developed in the Python language and incorporat-ing libraries of the ROOT data analysis framework for the state-of-the-art mass spectrometry method called phase-imaging ion-cyclotron-resonance (PI-ICR). A step-by-step description of the dataset construction and analysis algorithm is given. The code features a new phase-determination approach that offers up to 10 times smaller statistical uncertainties. This improvement in statistical uncertainty is confirmed using extensive Monte-Carlo simulations and allows for very high-preci-sion studies of exotic nuclear masses to test, among others, the standard model of particle physics. Source code: https://github.com/jonas-ka/pi-icr-analysis Zenodo code DOI: https://doi.org/10.5281/zenodo.4553515 Programming language: Python Feb. 20, 2021J. Karthein * ,@ , D. Atanasov , K. Blaum , D. Lunney , V. Manea , M. Mougeot CERN, 1211 Geneva 23, Switzerland Max-Planck-Institut für Kernphysik, 69117 Heidelberg, Germany IJCLab, CNRS Université Paris-Saclay, 91406 Orsay, France *: present address:
Massachusetts Institute of Technology, MA-02139 Cambridge, USA @ : corresponding address: [email protected] Analysis m ethods and code for very high-precision m ass m easure m ents of unstable isotopes I-ICR analysis software February 2021 J.Karthein et al.
Groundbreaking new measurement technique for Penning traps was developed by S. Eliseev et al. [El13] at GSI based on preliminary studies by Eitel et al. [Ei09] in Mainz. This so-called phase-imag-ing ion-cyclotron-resonance (PI-ICR) technique allows for the direct mass determination in a non-scanning approach. It thus allows for up to 25 times faster measurements [El14, Ka19b] at the same or better levels of precision. This technique was successfully implemented at the ISOLTRAP setup [Ka19b, MaKa20] for which an analysis code was developed in Python [Ka19c]. Given an optimal beam characteristics and high statistics, a new analysis approach presented be-low allows for a multifold reduction of the statis-tical uncertainty compared to the state of the art without any intervention on the experimental apparatus. It can even be applied to all already-existing datasets and opens the door for very-high-precision mass determination required for stringent tests of the Standard Model with uncer-tainties in the order of a few (tens) eV/ .
PI-ICR THEORY — The mass determination in a Penning trap is based on the direct determina-tion of the cyclotron frequency , In the case of an ideal Penning trap, is also equal to the sum of the frequencies of the radial eigenmotions, and , of the trapped ion [Br86] (also, see Fig. 3). The magnetic field strength can be eliminated with the ratio , between the cyclotron frequency of a well-known reference ion and the cyclotron frequency the ion of interest as long as these consecu-tive measurements are performed within a time-frame shorter than non-linear fluctuations of the magnetic field (typically ~15 min) [Ke03]. In order to reduce the statistical uncertainty, many of such measurement pairs are performed under different measurement conditions to elim-inate eventual systematic effects. The final ratio is calculated by fitting simultaneously the same polynomial function describing the temporal evolution of the cyclotron-frequency measure-ments, and , respectively: The degree of the polynomial function is adjusted to the lowest degree describing the en-vironment-temperature-induced fluctuations of the magnetic field. One can minimize the reduced of the data distribution to the fit, which typi-cally results in the same polynomial degree . This ratio-determination method has proven to be very robust [Na93, Ka19a, Ka19b, Fi12, El15] c ν c ν c ν + ν − B r ν c ,ref ν c ,ioi rf n ( t ) ν c ,ioi ν c ,ref n f n ( t ) χ n Page of 2 14
Time −0.5−0.4−0.3−0.2−0.10.00.10.2 F r e q u e n c y ( H z ) +4.33785e6 2.93.03.13.23.33.43.53.6+4.33706e6poly-6: χ red Ne Napoly-fit
Fig. 1: Frequency-ratio determination based on the tempo-ral cyclotron-frequency distributions of the ion of interest (red) and the reference ion (blue), as well as simultaneous polynomial fits. The frequency scale differs between the two species to allow for a better display of both data sets. The displayed data is a subset of the data published in [Ka19a]. For details, see text. ν c = ν + + ν − = B π × qm . r = ν c ,ref ν c ,ioi , ν c ,ioi ( t ) = f n ( t ) , ν c ,ref ( t ) = r × f n ( t ) . (1) I-ICR analysis software February 2021 J.Karthein et al. while minimizing the statistical uncertainty and is therefore recommended over the commonly used linear neighbor-to-neighbor interpolation. The polynomial technique offers the clear advantage of modeling the nonlinear magnetic field drifts of the whole distribution instead of just interpolat-ing the (assumed) linear trend to the average time of the measurement. This longterm model-ing is depicted in Fig. 1 for an 8 hrs measurement. The atomic mass of interest , can finally be calculated from the ratio of the cyclotron frequencies of the (singly charged) ions in the Penning trap, the atomic mass of the refer-ence ion from literature [Hu17], as well as the electron mass [St14] and the corresponding binding energies of the valence electron of either atom. These atomic binding energies are typically on the order of a few eV and hence can no longer be neglected at the precision level achieved in this work. The PI-ICR measurement technique makes use of the fact that the mass determination can be re-duced to a radial frequency determination in a Penning trap. In PI-ICR, a radial frequency is derived from the full phase accumulated by an ion in a very precisely measured time : This full phase consists of an integer number of full turns as well as an additional phase between the starting point and the end point of the motion. In practical situations, the integer number of turns can be determined based on a prior (much less precise) knowledge of the ion's mass, but this knowledge can also be built up entirely by PI-ICR. This can be achieved by starting with the time as a fraction of one full period and then stepwise increasing . The more is increased, the smaller also the uncertainty on becomes since the sta-tistical uncertainty of in this case only origi-nates from the determination of the additional phase . One full period for the radial eigenmo- m ioi rm e E B ν r φ full t r nn ∈ ℕ φ r ∈ [0, 2 π [ t r T r = ( t n − t )/ nt r t r ν r ν r φ r Page of 3 14 F u ll m o t i o n o f i o n s i n P e nn i n g t r a p d u r i n g g i v e n t i m e t a cc = t e n d - t -16 -8 0 8 16 Detector position X (mm) -16-80816 D e t e c t o r p o s i t i o n Y ( mm ) P ( t ) = P ( t n ) C Integer number n of full turns + -16 -8 0 8 16 Detector position X (mm) -16-80816 D e t e c t o r p o s i t i o n Y ( mm ) ϕ Final (additional) phase P ( t ) P ( t end ) C φ φ r r -16 -8 0 8 16 Detector position X (mm) -16-80816 D e t e c t o r p o s i t i o n Y ( mm ) P ( t ) = P ( t n ) C Integer number n of full turns + -16 -8 0 8 16 Detector position X (mm) -16-80816 D e t e c t o r p o s i t i o n Y ( mm ) ϕ Final (additional) phase P ( t ) P ( t end ) C φ φ r r Fig. 2: Graphical description of the determination of the full phase for a radial motion using the PI-ICR technique. The full phase of a trapped ion consists of an integer number of full turns (top) as well as a final, additional phase (bot-tom). This final phase of the radial motion ex-tends between the spot position (red) of the ions at the starting time , the position of the ions (highlighted in blue) after the evolu-tion time , and a reference point representing the center (black) of the radial mo-tion (dashed line). φ full n ∈ ℕ φ r ∈ [0, 2 π [ φ r P ( t ) t P ( t end ) t r = t end − t C ν r = φ full π t r = 2 π n + φ r π t r . m ioi = r × ( m ref − m e − E B ,ref / c ) + m e + E B ,ioi / c , (2) I-ICR analysis software February 2021 J.Karthein et al. tions in a typical Penning-trap setup is on the or-der of: ~ μ s ; ~ ms . This basic procedure is depicted in Fig. 2 where the additional phase (highlighted in green) is defined by the position (highlighted in red) of the ions at the starting time , the position of the ions (highlighted in blue) after the evolution time inside the trap at frequency , as well as a reference point (highlighted in black) representing the center of the radial motion (dashed line). In case of a Pen-ning trap, a measurement scheme was developed by S. Eliseev et al. [El14] for the direct determina-tion of the cyclotron frequency in only three measurement steps: 1. The ion motion is cooled by radial centering in the Penning trap via resonant RF-excitation [Kr99]. This position is captured by ejecting the ion onto an on-axis 2D-position-sensitive micro-channel-plate (MCP) detector, which is parallel to the radial eigenmotions of the ion inside the trap (see Fig. 3). This step allows for a calibration of the radial motions and only has to be repeated once every couple of hours since this position is very stable over time [El14]. 2. The ion is prepared on a pure motion by resonant RF-excitation [Kr99] and after evolu-tion at for precisely the accumulation time the position is recorded as in step 1. Dif-ferent ions can be prepared over several cy-cles at exactly the same starting point in order to be able to repeat this step for increased statistics and thus reduced sta-tistical error. Before ejection onto the detec-tor, the ion motion is converted into a (slow) motion by RF-excitation to reduce the an-gular spread on the detector. This conversion hereby conserves the absolute value of the accumulated phase but switches the sign of the direction of the phase evolution [El14], hence the plus sign for in Eq. (3). 3. The ion is prepared on a pure motion as in step 2 and immediately converted into a pure motion by resonant RF-excitation [Bl03, Ei09]. After evolution at for precisely the same accumulation time the position is recorded as in step 1 and 2. It is important to also use the same starting point as for step 2 in order for Eq. (3) to be valid. T ν + T ν − φ r P ( t ) t P ( t end ) t r = t end − t ν r C ν c ν + ν + t acc P ( x , y ) ν − φ final = φ + + φ − ν + ν − ν − t acc P ( x , y ) Page of 4 14
Electron cloudMCPDelay linesIon t t X1 t X2 t X1 t X2 t Y2 t Y1 Penning trap B ⃗ Projection & magni fi cation ν + ν — ν z Ion motion inside a Penning trap:= superposition of three eigenmotions ( ν + , ν — , ν z ): Fig. 3: PI-ICR detection principle. The ion motion inside the Penning trap is projected onto the detector while being magnified due to the magnetic field gradient. At the detector, the individual ion is converted into an electron cloud when hitting the micro-channel plate (MCP). The resulting signal triggers the first time event. The electron cloud is then accelerated onto a delay-line array (here: Roentdek DLD40), where the induced signal in the X and Y wires propa-gates towards their extremities, generating in total four additional time stamps. These five time stamps allow for a precise position detection.
I-ICR analysis software February 2021 J.Karthein et al.
Each of these three measurement steps is re-peated multiple times per measurement cycle in order to allow for a determination of the mean position of the so-called spots on the detector (see Fig. 2). These steps reduce the determination of the cyclotron frequency , to the determination of the final phase , de-fined by the spots in steps 1, 2 and 3 instead of requiring one additional spot for each frequency at thanks to the same preparation scheme.
ANALYSIS CODE — RECONSTRUCTION
The analysis is divided into three parts. First, the raw data is read and transformed into position information. Second, the position information is used to de-termine the total phase and the cyclotron fre-quency. Third, the mass of the ion of interest is calculated based on the temporal distribution of cyclotron frequencies. Thanks to the Jupyter notebook environment [Kl16], all parts of the code are easy to read while being fairly customizable. They are independent and can be adjusted to the individual require-ments of different experiments or used for a dif-ferent measurement technique. As an example, the third part concerning the determination of the mass based on the simultaneous fitting tech-nique can also be used for all Penning-trap mass spectrometry techniques based on frequency-ra-tio determination, such as the time-of-flight ion-cyclotron-resonance (ToF-ICR) technique in it’s single-pulse [Gr80] or Ramsey-type application [Ge07] with no adjustments. The first part of the code deals with the event-position reconstruction by importing the five time stamps (MCP, X1, X2, Y1, Y2) from the LabView data acquisition. These timestamps are visualized in Fig. 3. The MCP time is used for a calibration of the four time signals in the X and Y direction of the detector by a simple subtraction. Next, the difference between the two time stamps per di-rection gives the position on the detector along this direction. The detector should be calibrated once by stepwise increasing the excitation ampli-tude of one of the radial ion motions until no ions are detected. This allows to set the proper "time windows" for the individual time stamps marking the edges of the detector. Since it is possible to capture more than one ion inside the trap, the data acquisition at ISOLTRAP is configured such that all events in one mea-surement sequence are saved in their individual channel (MCP, X1, X2, Y1, Y2; see Fig. 3) after being converted to NIM standard by a constant-fraction discriminator. Due to the pulse-height distribution of the sig-nals and limits on minimum discriminator threshold, sometimes one or two of the five re-quired time stamps per position are missing. This leads to a disturbance in the order of incoming events, which is why a computationally demand-ing comparison between all individual time stamps per measurement sequence has to be performed. This calculation ensures finding the real events within all acquired time stamps by performing the sum check described above. The combinatory computation is performed with a cascaded for-loop, which compares each time stamp in all of the five channels with each other. Smart if-conditions are imposed to skip false ν c φ final t Page of 5 14 ν c = ν + + ν − = 2 π n + + φ + π t acc + 2 π n − + φ − π t acc = 2 π ( n + + n − ) + φ final π t acc , (3) I-ICR analysis software February 2021 J.Karthein et al. events and thus reduce the calculation time. This part of the code is CPU-parallelized where appro-priate. A hardware-accelerated (GPU) alternative using 5-dimensional meshgrid tensors (one di-mension each for X1, X2, Y1, Y2, MCP, and thus re-placing the need for the combinatorics based on for-loops), tensor masking (replacing the compar-ison- and if-conditions) and just-in-time compila-tion, all based on the machine learning library JAX from Google [Br18], was tested and found to be about 10 times slower on the available office-PC hardware due to the much-increased number of computations per event. But since the data acquisition is limited in measurement time to en-sure linearity in the magnetic-field drifts (on the order of ~15 min), the position reconstruction typ-ically only takes about one second per file for about tens or hundreds of files in a typical exper-iment. The position information is then saved in a new *.csv-file such that this step only has to be performed once.
ANALYSIS CODE — FREQUENCY EVALUATION — The next part in the analysis is to translate the posi- tion information into a phase information and from that to calculate the cyclotron frequency of the ions. It was purposely chosen to perform the analysis on a file-by-file basis rather than a batch analysis to ensure full control over all cuts and fits. Furthermore, the Pandas library [Mc10] was used extensively to ensure maximal readability and adjustability of the code. In a first step, sev-eral analysis cuts can be performed: -
ToF-cut : The time-of-flight (ToF) distribution from the trap to the detector allows for a limi-tation in the mass of interest to remove cont-aminants of different masses. ToF outliers are reduced automatically within a 10-Sigma range. -
Position cut : The position region of interest can be limited to the desired area of the detector (e.g., to cover only the blue spot in the top part of Fig. 2). Moreover, position outliers are de-tected and reduced automatically within a 10-Sigma range of the mean position. This cut also allows for the removal of contaminating species (such as single isomeric states; see Ref. [MaKa20]).
Page of 6 14
Phase (rad) R a d i u s ( mm )
2D multivariate polar fit
MC-simulated3- σ band (fit)fit X (mm) Y ( mm ) Cartesian transformation
X (mm) Y ( mm ) ×
2D projection Cartesian fit X MC = 7.07 Y MC = 7.07 X = 7.09 ± 0.02 Y = 7.03 ± 0.02 X = 6.85 ± 0.02 Y = 6.80 ± 0.02 a) b) c) Fig. 4: Multivariate Gaussian 2D-fit for a Monte-Carlo (MC) simulated circular distribution in polar coordinates (a) and its trans-formation into Cartesian coordinates (b). For a visual comparison, the corresponding 2 ⨉
1D Gaussian fit is shown in c). Its deviation from the theoretical values compared to the fit uncertainties (without correction for correlation in case of the projection fit in c) is significant while the multivariate fit in polar coordinates (a, b) agrees with the theo-retical values . Spot shapes like the one shown in b) can occur, e.g., due to timing jitters during the ejection procedure or in the direct projection of the (fast) motion. X MC , Y MC X MC , Y MC ν + I-ICR analysis software February 2021 J.Karthein et al. - Z-cut : Limiting the number of ions per capture to reduce space-charge effects [Ke03]. In each step, the number of removed spots are printed transparently into the output console. Next, the positions for the reference, , and spots have to be determined. The approach so far with the SHIPTRAP LabView code [Ka20] was using a 2 ⨉
1D Gaussian model in X and Y directions to fit the unbinned position data (see Fig. 4c). In our approach, we are using a 2D multivariate Gauss-ian model with the maximum likelihood estima-tion for the unbinned position data. The differ-ence between the two models is that the multi-variate version properly accounts for two-dimen-sionality of the distribution with its rotation in space, whereas the regular 2 ⨉
1D Gaussian model merely fits the individual projections onto each dimension. In a perfect experiment with a circular distribution in the X-Y-space there will be no dif- ν + ν − Page of 7 14
X (mm) Y ( mm ) X (mm)
X (mm) C o u n t s Phase (rad)
Phase (rad)
Phase (rad) φ MC = 0.7854 φ = 0.7854 ± 0.0126 φ MC = 0.7854 φ = 0.7855 ± 0.0158 φ MC = 0.7854 φ = 0.7855 ± 0.0154 φ = 0.7854 ± 0.0022 δφ / δφ = 5.72 φ = 0.7855 ± 0.0044 δφ / δφ = 3.59 φ = 0.7780 ± 0.0036 δφ / δφ = 4.28 a) b) c)d) e) f) N data = random data3- σ band (fit)fit Fig. 5: Monte-Carlo (MC) simulation of three-spot scenarios for a symmetric, aligned elliptical, and rotated elliptical (with re-spect to the ion motion along the dashed circle) 2D Gaussian distribution with 1000 ions, as well as their multivariate fit in 2D (top row) and their phase-projected 1D fit (bottom row) in comparison. Elliptical spots as shown in a) and b) can occur due to asymmetries and imperfections in the trap or transport potentials. The true phase is compared to the average phase resulting from the mean X/Y coordinates from the multivariate 2D fit and the average phase from the 1D fit. Both averages are calculated from the fit results of 1000 random MC samples for each of the three dis-tributions. The given uncertainties correspond to the average one-sigma uncertainty for the given quantity over the 1000 simulations. Additionally, the reduction factor in statistical uncertainty from the 2D to the 1D fit is shown for the three spot scenarios. The symmetric and aligned spots agree well with the true phase value, whereas the 1D fit for the rotated elliptical spot shows a >2 sigma deviation from the true value. φ MC φ φ I-ICR analysis software February 2021 J.Karthein et al. ference between the two models. Deviations will however occur in reality. If we assume an exag-gerated case with a highly elliptical or asymmet-ric spot at a rotation of 45 degrees (as depicted in Fig. 4b), a 2 ⨉
1D Gaussian model would result in a circular distribution instead of an elliptical one since it essentially fits only the projections in X and Y. The phase of both fits might agree in case of a symmetric distribution, but the 2D-projection Gaussian model naturally results in highly corre-lated fitting parameters, which increase the off-diagonal elements of the fit’s covariance matrix. These correlations lead mathematically to a sig-nificant additional contribution in the statistical uncertainty, which typically is not accounted for. In the multivariate case on the other hand, corre-lation for PI-ICR spots at any ellipsoidal deforma-tion or rotation at any point of the circle is negli-gible due to the intrinsic rotation matrix in the model. Hence, the spot in Fig. 5b can be properly described as well as any other spot in the simula-tion. This model was chosen over the proposed alternative approaches described below due to its robustness and versatility for any spot shape. It builds on the spatial distribution based on ini-tial buffer-gas cooling, as derived in Eq. (12) ff. in [El14]. The minimization of the negative log-likelihood function for the maximum likelihood estimation is performed with the iMinuit library [On12], which ports the Minuit2 minimizer of CERN’s analysis framework ROOT [Br97] to Python. It was chosen over different SciPy [Vi20] minimizers due to its robustness and direct support at CERN. The library allows in addition the integration of au-tomatic differentiation, i.e., using exact deriva-tives instead of numerical derivation in the min- imization as available, e.g., in the JAX library [Br18]. This integration is planned in a future per-formance update, but is not essential for the analysis code at this stage. From the positions of the reference, , and spots one can determine the phase with the arctangent function "atan2" of the NumPy library [Ha20a], which has the whole circle as image set instead of just . This phase finally allows for the calculation of the cyclotron fre-quency according to Eq. (3). The number of full turns for each eigenfrequency has to be predetermined as described above using the stepwise increase of the accumulation time . All analysis data from the raw events to the fitting information is saved in a pickle-file (*.p) for traceability and transparency. All results are ad-ditionally saved in a comma-separated-value file (*.csv) for the independent calculation of the fre-quency ratio and then the mass.
ANALYSIS CODE — FREQUENCY-RATIO — The calcu-lation of the final mass is performed based on the exact time stamp of the measurement and the calculated and can thus also be used for different Penning-trap mass spectrometry tech-niques like the ToF-ICR technique in its single-pulse or Ramsey-type excitation scheme since they also use the calculated and the exact time stamp of the measurement for the frequency-ra-tio determination. The software uses a polynomi-al function according to Eq. (1) to describe the temporal distribution of all measured -val-ues using the least squares method (see Fig. 1). The statistical uncertainty is calculated from the covariance matrix of the fit, which is performed for all possible polynomial degrees with the number of data points to avoid ν + ν − φ final ] − π /2, π /2[ ν c n +, − t acc ν c ν c f n ( t ) ν c n ∈ [1, ν c − ν c ν c Page of 8 14
I-ICR analysis software February 2021 J.Karthein et al. overdetermining the fit. The software automati-cally chooses the fit with a closest to one for the final ratio , but the degree could also be fixed, e.g., to the known degree (if known) of the temporal fluctuations during the time of the measurement. The results of this analysis step are saved in a separate *.csv-file. From the final ratio , the mass can be calculated according to Eq. (2), using the latest literature mass for the reference ion [Hu17]. It is common to publish the final ratio to update if ever a more pre-cise is available.
ALTERNATIVE APPROACH 1: STRONG RADIAL DE-FORMATION — In case of strong radial deforma-tion of the Gaussian distribution in Cartesian co-ordinates along the ion motion (as depicted in Fig. 4c), a transformation for all spot positions into polar coordinates is highly recommended. As emphasized by Fig. 4a, the same distribution can be described much better in this coordinate sys-tem. Applying a regular multivariate Gaussian model to the data in Cartesian coordinates would result in a deviation to the actual mean position of the spot. In the code, a simple transformation to polar coordinates can be performed at the be-ginning of the "Spot fitting" section as well as back to Cartesian coordinates after the fitting. By this, no further adjustments have to be made to the code. An example is given in the subfolder "alternative-approaches".
ALTERNATIVE APPROACH 2: VERY HIGH-PRECISION — In the following, we present an alternative analysis approach that allows — given the pro-posed preparation — up to 10 times smaller sta-tistical uncertainty on the same dataset com-pared to the robust multivariate approach de-scribed above. This is a game changer for many applications, such as Weak-interaction Standard Model tests based on the determination of the element [Ka19a]. Given the exact same exper-imental operation and high statistics, relative statistical uncertainties of can be reached for the first time in a few hundred mil-liseconds measurement time only. The analysis approach is based on a transforma-tion of the spot positions into polar coordinates. This transformation results in a position informa-tion separated into a phase coordinate and a ra-dial coordinate instead of X and Y. Since in PI-ICR one is interested in the mean position of the dis-tribution along the phase coordinate, the ap-proach only uses the phase information and ne-glects the radial position information. This reduc-tion has the advantage that the number of free parameters in the maximum likelihood estima-tion reduces from five, in case of a multivariate 2D Gaussian model (mean X, mean Y, sigma X, sigma Y, rotation angle used in the covariance matrix) to only two free parameters (mean phase, sigma phase). We performed extensive Monte Carlo (MC) simu-lations to study the effect and accuracy of this new approach for many different, realistic sce-narios by varying the following properties: - spot shapes: ellipticity from 1 (Fig. 5a) to 2 (Fig. 5b/c); larger ellipticities are not recommended due to lack of control over systematics - rotation with respect to the tangent to the cir-cular trajectory of the spot center: 0 (see Fig. 5a/b) to 45 degrees (see Fig. 5c); behavior re-peats afterwards - position on the detector: 0 to 45 degrees with respect to the horizontal axis and counted counterclockwise; behavior repeats afterwards χ rr m ioi m ref r m ioi m ref V ud δν c / ν c ≤ − Page of 9 14
I-ICR analysis software February 2021 J.Karthein et al. - ratio between trajectory radius and spot diam-eter: 10 mm / 2 mm up to 15 mm / 2 mm - number of ions per spot: 10, 100, 1000, 10000 Each MC simulation was repeated 1000 times in order to calculate average fit result values, aver-age fit uncertainties, the variance of the individ-ual fit results, and the skewness of the fit results scattering to ensure the statistical significance of the following statements. An example code is given in the subfolder "alternative-approaches". The results clearly demonstrate the power of this approach as long as the spot is reflection sym-metric with respect to a radial axis passing through the trajectory center. This is depicted in Fig. 5a/d and 5b/e. The top row shows the initial unbinned MC-random distribution (here for 1000 ions per spot) and their unbinned maximum like-lihood estimation using the multivariate 2D Gaussian models described in the "ANALYSIS CODE — FREQUENCY EVALUATION"-section. Their results are displayed in red for the mean X/Y po-sition and in blue for the 3-sigma ellipse. The evaluation agrees with the true value while pro-viding a small statistical uncertainty for the posi-tions parameter. From the average X and Y posi-tion information, a phase can be calculated as described in the same section mentioned above. In the bottom row, the position information dis-played in the top row is transformed point by point into polar coordinates and the unbinned phase information is shown. In case of experi-mental data, the position information would first need to be calibrated by the position of the cen-ter of the circular motion. This average position of the measurement’s step 1 must be determined in any approach first with the 2D multivariate fit. The phase information is also displayed in a binned manner to guide the eye. The results of the unbinned maximum likelihood estimation of the 1D phase data is displayed in red for the mean phase position and in blue for the 3-sigma range. For the mirror-symmetric cases with re-spect to the radial axis, i.e., column one and two, both the new approach and the multivariate ap-proach agree on average in the phase determina-tion with the theoretical value for any number of simulated ions per spot: 10, 100, 1000, 10000. However, the new approach reaches much-re-duced statistical uncertainties as displayed in the uncertainty ratio between the 2D ap-proach versus the 1D approach. These reduction factors for the average statistical uncertainty found in the MC simulations performed for all mentioned spot conditions are presented in Tab. 1 for different numbers of ions per spot. The simulations were performed for 10000, 1000, 100 and 10 ions to cover the range of statistics found in a realistic experiment. For a very low number of ions (10-100) as would be accumulated in a typical experiment at the production limita-tion of present RIB facilities [Ka20] (few ions per δφ / δφ Page of 10 14
Tab. 1: Statistical uncertainty improvement factors for the 1D-phase fitting approach compared to the 2D multivariate Gaussian fitting approach for different number of ions and circular/elliptical mirror-symmetrical shapes (see left and middle columns of Fig. 5). The improvement factors are calculated for the mean statistical uncertainties of 1000 MC simulations. The improvement factor ranges correspond to the evaluation at different angles (0-45 degrees) and realistic ratios between the trajectory radius and spot diameter (2 σ ) of 5-8. Improvement factorscircle ellipse10 1.8 1.2 - 1.3100 3.2 2.0 - 2.31000 5.6 - 5.7 3.6 - 4.010000 9.8 - 10.1 6.4 - 7.1 ion
I-ICR analysis software February 2021 J.Karthein et al. second/minute/hour), the reduction in statistical uncertainty is moderate between a 20% and over 300% improvement. However, it is much more in-teresting for experiments at higher production rates where typical Standard Model tests are per-formed. In these situations where much more data can be accumulated, significant improve-ment factors on the statistical uncertainty of 4-10 can be achieved just by applying the proposed fitting approach. While the robust 2D multivariate fitting approach can be applied for any type of spot deformation introduced by different kinds of systematic ef-fects (e.g., jitter in the timing, fringe fields, or edge effects in the transport to the detector, fluc-tuations in the potentials, …), the 1D phase fitting approach will result in a systematic deviation from the true value as soon as the spot is no longer symmetric to the motion axis. This situa-tion is emphasized in Fig. 5c/f for a spot with 45 degree rotation and an ellipticity of 2 (in case of stronger deformed spots, better ion preparation as described in Ref. [El14] is recommended). This systematic absolute deviation is shown in Fig. 6 for different spot rotation angles and dif-ferent number of ions per spot. Each data point represents the mean deviation for 1000 MC simu-lations and their standard deviation as error bar. Independent of the number of ions in the spot, a maximal absolute deviation of ~0.01 rad was found for the most extreme deformation and ro-tation. The resulting relative shift for the cy-clotron frequency of an ion of mass m ~ 100 u in a Penning trap of 6 T ( → ~1 MHz) measured at typical 100 ms and 1000 ms phase accumulation time would amount to and , respectively. This is on the order of the latest es-timate on the upper limit of the systematic un-certainty of for the ISOLTRAP setup as described in Ref. [Ke03] . This limit will be revisit-ed during planned systematic studies in the fol-lowing months. Thus, in typical mass determina-tions of short-lived isotopes with low production ν c × − × − × − Page of 11 14
Spot rotation angle (°) −0.020.000.020.040.060.08 0 5 10 15 20 25 30 35 40 45
Spot rotation angle (°)
Spot rotation angle (°) max at 32° N spot = 10: Δ = 0.010 ± 0.021 δφ = 0.03 ... 0.06 N spot = 100: Δ = 0.010 ± 0.007 δφ = 0.01 ... 0.02 N spot = 1000: Δ = 0.009 ± 0.002 δφ = 0.003 ... 0.006 Δ = φ M C − φ D ( r a d ) a) b) c) Fig. 6: Systematic absolute deviation of the 1D phase fit for an elliptical spot of ellipticity 2 and for different spot rotation an-gles with respect to the the tangent to the circular trajectory of the spot center. The results are shown for 10 (a), 100 (b) and 1000 (c) ions per spot. Each deviation represents the mean of 1000 MC simulations and is plotted with its standard deviation as error bar. The behavior is compared to a sine function (red) from which the maximal deviation (~0.1 rad) and the rotation angle for the maximal deviation (~32 degrees) were derived. For comparison, the mean individual sta-tistical uncertainty for 1000 repeated MC simulations is shown with its range for the given spot rotation angles. Δ δφ I-ICR analysis software February 2021 J.Karthein et al. yields (few ions per second/minute/hour) and short half-lives (~100s ms), the moderate statisti-cal improvement factors of the 1D phase ap-proach between 20% and over 300% might be at the cost of a possible additional systematic un-certainty from the fit which is on the same order of magnitude. Since these few number of ions per spot make a symmetry judgement difficult even via skewness calculations of the distribution, the robust 2D multivariate approach is recommend-ed. A different conclusion can however be made for so-called Q -value measurements, in which not the actual mass but the mass difference between two isobars is determined. These measurements are typically used for high-precision Standard Model tests since they allow for significant can-cellation of most systematic effects. Hereby, the two ions of interest of same mass number are consecutively prepared at the exact same spot in the trap und thus undergo the exact same sys-tematic shifts. While high statistics would make an optimization for symmetry easy (and is in any situation highly recommended), the fit-induced systematic shift in the individual mass determi-nation would however cancel out even for asym-metric and rotated spots since the shift will be on the same order and direction for both isobars (rigorously true at least for spot distortions aris-ing from transport optics effects). In summary, this new phase fitting approach presents a clear advantage over other existing and the presented robust 2D phase determination method specifi-cally for Q -value measurements. This is true for any number of ions and any spot shape with sta-tistical uncertainties up to 5-10 times smaller than with existing techniques. CONCLUSION/OUTLOOK — The recently developed PI-ICR method has become the state-of-the-art mass spectrometry method for radioactive iso-topes. In this publication, two new analysis schemes were presented for the first time, offer-ing up to 10 times smaller statistical uncertain-ties even on existing datasets. This improvement opens the window for next-generation Standard Model tests and significant contributions to the determination of the neutrino mass via the end point in the electron-capture spectrum using short-lived isotopes at low energy.
ACKNOWLEDGEMENT — We acknowledge financial support by the Max-Planck Society. JK acknowl-edges financial support by a Wolfgang Gentner Ph.D. Scholarship of the BMBF (05E15CHA) and the ISOLDE collaboration for transitional support during the COVID-19-pandemic.
REFERENCES [Bl02]
Blaum, K. et al. "Carbon clusters for ab-solute mass measurements at ISOLTRAP"
Europ. Phys. J. A (2002), 245. doi.org/10.1140/epja/i2001-10262-4 [Bl03] Blaum, K. et al. "Recent developments at ISOLTRAP: towards a relative mass accuracy of exotic nuclei below 10 -8 " J. Phys. B (2003), 921. doi.org/10.1088/0953-4075/36/5/311 [Bl06] Blaum, K. "High-accuracy mass spec-trometry with stored ions."
Phys. Rep. (2006), 1. https://doi.org/10.1016/j.physrep.2005.10.011 [Br18]
Bradbury, J. et al. "JAX: composable transformations of Python+NumPy programs"
GitHub (2018) github.com/google/jax [Br86]
Brown, L., and Gabrielse, G. "Geonium theory: Physics of a single electron or ion in a Penning trap."
Rev. Mod. Phys. (1986), 233. doi.org/10.1103/Rev-ModPhys.58.233 Page of 12 14 I-ICR analysis software February 2021 J.Karthein et al. [Br97]
Brun, R. et al. "ROOT - An Object Orient-ed Data Analysis Framework."
Nucl. Inst. & Meth. in Phys. Res. A
389 (1997), 81-86. doi.org/10.1016/S0168-9002(97)00048-X [Ei09]
Eitel, G. et al. , "Position-sensitive ion detection in precision Penning trap mass spectrometry."
Nucl. Instrum. Meth. Phys. Res. A (2009), 475. doi.org/10.1016/j.nima.2009.04.046 [El13]
Eliseev, S. et al. "Phase-imaging ion-cy-clotron-resonance measurements for short-lived nuclides."
Phys. Rev. Lett. (2013), 082501. doi.org/10.1103/PhysRevLett.110.082501 [El14]
Eliseev, S. et al. "A phase-imaging tech-nique for cyclotron-frequency mea-surements."
Appl. Phys. B (2014), 107–28. doi.org/10.1007/s00340-013-5621-0 [El15]
Eliseev, S. et al. "Direct Measurement of the Mass Difference of
Ho and
Dy Solves the Q -Value Puzzle for the Neutrino Mass Determination." Phys. Rev. Lett. (2015), 062501. doi.org/10.1103/PhysRevLett.115.062501 [Fi12]
Fink, D. et al. " Q -value and half-lives for the double- β -decay nuclide Pd"
Phys. Rev. Lett. (2012), 062502. doi.org/10.1103/PhysRevLett.108.062502 [Ge07]
George, S. et al. "Ramsey method of sep-arated oscillatory fields for high-pre-cision Penning trap mass spectrome-try."
Phys. Rev. Lett. (2007), 162501. doi.org/10.1103/PhysRevLett.98.162501 [Gr80] Gräff, G. et al. "A direct determination of the proton electron mass ratio."
Zeit. Phys. A (1980), 35–39. doi.org/10.1007/BF01414243 [Ha20a]
Harris, C. R. et al. "Array programming with NumPy."
Nature (2020), 357–62. doi.org/10.1038/s41586-020-2649-2 [Ha20b]
Hardy, J. C. & Towner, I. S. "Superallowed 0 + → + nuclear β decays : 2020 criti-cal survey , with implications for V ud and CKM unitarity". Phys. Rev. C (2020), 045501. doi.org/10.1103/Phys-RevC.102.045501 [Hu17]
Huang, W. J. et al. "The AME2016 atomic mass evaluation (I). Evaluation of in-put data; And adjustment procedures."
Chinese Phys. C (2017). doi.org/10.1088/1674-1137/41/3/030002 [Ka19a] Karthein, J. et al. "Q EC -value determina-tion for Na → Ne and Mg → Na mirror-nuclei decays using high-pre-cision mass spectrometry with ISOLTRAP at the CERN ISOLDE facility."
Phys. Rev. C (2019), 15502. doi.org/10.1103/physrevc.100.015502 [Ka19b]
Karthein, J. et al. "Direct decay-energy measurement as a route to the neu-trino mass."
Hyperfine Interact. (2019), 61. doi.org/10.1007/s10751-019-1601-z [Ka19c]
Karthein, J. "PI-ICR analysis software."
Zenodo/Github (2019). doi.org/10.5281/zenodo.3965766 [Ka20]
Kaleja, O. "First direct mass measure-ment of a superheavy element and high-precision mass spectrometry of nobelium and lawrencium isotopes and isomers using the Phase-Imaging Ion-Cyclotron-Resonance technique at SHIPTRAP."
Dissertation, Universität Mainz (2020). [Ke03]
Kellerbauer, A. et al. "From direct to ab-solute mass measurements: A study of the accuracy of ISOLTRAP".
Eur. Phys. J. D (2003), 53–64. doi.org/10.1140/epjd/e2002-00222-0 [Kl13] Kluge, H.-J. "Penning trap mass spec-trometry of radionuclides."
Int. J. Mass Spectr. (2013), 26. https://doi.org/10.1016/j.ijms.2013.04.017 [Kl16]
Kluyver, T. et al. "Jupyter Notebooks – a publishing format for reproducible computational workflows."
Pos. Pow. Ac. Pub. (2016). doi.org/10.3233/978-1-61499-649-1-87 [Kr99] Kretzschmar, M. "A quantum mechanical model of Rabi oscillations between two interacting harmonic oscillator modes and the interconversion of modes in a Penning trap."
AIP Confer-ence Proceedings . (1999), 242–51. doi.org/10.1063/1.57446 [Lu18] Lunney, D. "Extending and refining the nuclear mass surface with ISOLTRAP."
J. Phys. G: Nucl. Part. Phys. (2017), 064008. https://doi.org/10.1088/1361-6471/aa6752 Page of 13 14 I-ICR analysis software February 2021 J.Karthein et al. [MaKa20]
Manea, V. & Karthein, J. et al. "First Glimpse of the N=82 Shell Closure below Z=50 from Masses of Neutron-Rich Cadmium Isotopes and Isomers."
Phys. Rev. Lett. (2020), 092502. doi.org/10.1103/PhysRevLett.124.092502 [Mc10]
McKinney, W. "Data structures for statis-tical computing in Python."
Proc. 9th Python in Science Conf . (2010) 56–61. [Na93]
Natarajan, V. et al. "Precision Penning trap comparison of nondoublets: Atomic masses of H, D, and the neu-tron."
Phys. Rev. Lett. (1993), 1998. https://doi.org/10.1103/PhysRevLett.71.1998 [On12] Ongmongkolkul, P. et al. "iminuit - A Python interface to MINUIT" GitHub (2012) github.com/iminuit/iminuit [St14]
Sturm, S. et al. "High-precision mea-surement of the atomic mass of the electron."
Nature (2014), 467–70. doi.org/10.1038/nature13026 [Vi20]
Virtanen, P. et al. "SciPy 1.0—fundamen-tal algorithms for scientific comput-ing in Python."
Nature Meth.17