Mahzad Kalantari

Franck Jung

Nicolas Paparoditis

Jean-Pierre Guedon

ROBUST AND AUTOMATIC VANISHING POINTS DETECTION WITH THEIR UNCERTAINTIES FROM A SINGLE UNCALIBRATED IMAGE, BY PLANES EXTRACTION ON THE UNIT SPHERE

Keywords: Photogrammetry, Geometry, Vanishing points extraction, Image Orientation, Error estimation

This paper deals with the retrieval of vanishing points in uncalibrated images. Many authors did work on that subject in the computer vision field because the vanishing point represents a major information. In our case, starting with this information gives the orientation of the images at the time of the acquisition or the classification of the different directions of parallel lines from an unique view. The goal of this paper is to propose a simple and robust geometry embedded into a larger frame of image work starting with an efficient vanishing point extraction without any prior information about the scene and any knowledge of intrinsic parameters of the optics used.After this fully automatic classification of all segments belonging to the same vanishing point, the error analysis of the vanishing points found gives the covariance matrix on the vanishing point and on the orientation angles of the camera, when using the fact that the 3D directions of lines corresponding to the vanishing points are horizontal or vertical. A validation of estimated parameters with the help of the photo-theodolite has been experimented that demonstrate the interest of the method for real case. The algorithm has been tested on the database of a set of 100 images available on line.

INTRODUCTION

In the conical geometry characteristic of human vision, or of photography, the parallel lines in the object space result in the image as pencils that intersect on vanishing points. Therefore to every vanishing point, a 3D direction is associated. The intersection of the lines in the image, corresponding to a given vanishing point, is mathematically an ill-conditioned problem. The various lines intersect with very small angles, which badly conditions the classic model of intersection. Another problem is a problem of classification, namely which segment is associated to which vanishing point. Since the years 80 many authors have tried to provide solutions to these two difficulties. In order to summarize the difference between all these approaches, one could say that it resides in the differences of parameter spaces for modelling and resolving the problem. Several types of spaces for modelling exist. One of the most widely known is the accumulation space on the Gauss sphere, that has been introduced by [START_REF] Barnard | Interpreting perspective images[END_REF] and that since then has been used, modified, and optimized by many different authors (Magee and Aggarwal,1984) [START_REF] Shufelt | Performance evaluation and analysis of vanishing point detection techniques[END_REF]. Other authors used other accumulation spaces using the Hough transform (Lutton et al., 1994) (Quan and Mohr, 1989) [START_REF] Tuytelaars | A cascaded hough transform as an aid in aerial image interpretation[END_REF]. The modelling can also be performed in the space object, which was used by many authors [START_REF] Teller | Automatic recovery of relative camera rotations for urban scenes[END_REF]. Then must be mentioned all authors who have modelled the problem in the 2D space of the image (van den Heuvel, 1998) [START_REF] Brauer Burchardt | Robust vanishing point determination in noisy images[END_REF] [START_REF] Schaffalitzky | Planar grouping for automatic detection of vanishing lines and points[END_REF] [START_REF] Rother | A new approach for vanishing point detection Pattern Anal[END_REF]. In general this kind of approach requires the use of probabilistic models [START_REF] Almansa | Vanishing point detection without any a priori information[END_REF]. Most methods of extraction of vanishing points are composed of 2 steps, one for the vanishing points detection, and another for the association of the different directions of the segments in the image with these vanishing points. The method presented in this paper has for parameter space an unit sphere without accumulation process. It works in only one step, so that the extraction is made simultaneously with the classification. The details of the geometry framework used in the method is described in the section 2. The different steps of the algorithm are described in the section 3. In order to have an assessment of the uncertainty on the location of the vanishing point, a very detailed covariance analysis has been done, and is the topics of the section 4. The section 5 illustrates the results of classification directly on some images taken as examples. Finally, in order to compare with a ground truth, a precise measure instrument that can give the orientations of the image acquisition was necessary, and thus a photo-theodolite (camera mounted with a precise goniometer) has been used.

GEOMETRIC FRAMEWORK

In this section, is presented the geometry on which the algorithm is based. The origin O is chosen at a distance D from the image plane. The axes Y and Z are parallel to the axes of the image, and thus the X axis is perpendicular to the image plane. l1 and l2 are the extremities of a given segment in the image. If this segment belongs to a family of segments who converges on the vanishing point V, one is interested in the corresponding family of planes of the space who pass through O, each containing one of these segments. V being on the line carrying the segment, a given plane contains therefore always V, and therefore the line OV too. For each of this family's planes, the normal vectors passing through O are all coplanar, because they are included in the plane passing by O and perpendicular to OV (Figure 1). And the normal in O to this last plane, (that we call "vanishing plane" in the present paper for this reason), intersects the image plane in the vanishing point V. Thus the extraction of V may be obtained through the research of the best plane containing all the extremities of normal vectors to the previously defined planes, each of them containing O, V a given segment. It is important to note that if O is chosen at the optic centre, and therefore if D is the focal length F, then the plane previously defined is the one that, in the literature of computer vision, is often called the interpretation plane [START_REF] Barnard | Interpreting perspective images[END_REF]) (Weiss et al., 1990)(Antone and[START_REF] Teller | Automatic recovery of relative camera rotations for urban scenes[END_REF]. But the present work doesn't stand in this very particular situation, which requires a previous knowledge of the intrinsic parameters of the camera. Here the point O is taken in an arbitrary way (Figure 2), obviously outside of the image plane, and for simple reasons of numeric stability, at a reasonable distance from the centre of the image. This remark is important, because the algorithms presented here require no prior knowledge of the optics used to extract correctly the points V. It works in the same way as the eye does, when it looks at an image, and rebuilds the vanishing points mentally without knowing the elements of the optics used to get this image.

(1)

N is the vector, passing through O, normal to the plane formed by l1, l2 and O. Each line of the ZY plane is defined by its polar equation (equation 1) as a function of and , polar coordinates of every segment in the image plane (using the origin of the image) (Figure 3). The normal vector to this plane can be calculated directly from and . The extremities of all these normalized normal vectors lie therefore on a unit sphere (equation 2).

That way, a mapping on this unit sphere has been performed (Figure 4). Some authors have used this geometry. But all of them took their origins on the optic centre, therefore at a distance F of the image, which requires the knowledge of the intrinsic parameters. Here, the problem of detection of the different families of directions of the image corresponding to parallel directions of the object space, is converted into the detection of all meaningful planes formed by the extremities of the normal vectors issued from O. In the next parts of this paper the extraction method of the planes is deepened. It is completed by an error analysis on the results obtained for the vanishing point.

ALGORITHM OF EXTRACTION OF THE VANISHING PLANES

Four important points of this algorithm can be mentioned: 1 -there is no need to make any hypotheses on the number of vanishing points present in the image, 2 -the algorithm is robust and use levels that adapt to the precision of detection of the segments, 3 -all vanishing points are detected; there is possibly an overdetection, but no under-detection, 4 -no prior information on the optics intrinsec parameters is requested.

The basic algorithm for the detection of the vanishing points requires two steps, followed if needed by an error analysis of the vanishing point.

The main steps are:

1 -Segments detection in image, 2 -Calculation of the normal vectors and extraction of the vanishing planes, 3 -Analysis and error propagation on the vector of vanishing point.

Detection of segments in the image

As it is the case for every algorithm of extraction of vanishing points, the first step is the detection of the image segments. Many algorithms of detection of segments are available [START_REF] Burns | Extracting straight lines[END_REF][START_REF] Weiss | An error analysis for surface orientation from vanishing points[END_REF])(Deriche et al., 1992). In the proposed method the detection of the segments is based on the Deriche Vaillant algorithm (Deriche et al., 1992). The advantage of this approach is that it can produce for every segment a covariance matrix for the line carrying this segment. The equation of this line is expressed in polar coordinates (½ and µ) and a covariance matrix is provided on these parameters.

But any other segments detector may also be chosen.

Calculation of the normal vectors and extraction of the vanishing planes

Once the segments detected on the image as well as the lines carrying these segments, the normal vectors can be easily calculated (equation 2). Now the problem lies in the extraction of the different planes in a set of clouds of points. For overcoming this problem, a method inspired of the RanSac(Fischler and Bolles,1987) was implemented. RanSac is a robust estimator based on a random sample principle. All necessary parameters for the use of the RanSac are described briefly in (Fischler and Bolles, 1987), this article being adapted for an unique cloud of points. But in the case of the vanishing planes, several clouds of points coexist. It is therefore necessary to bring some changes to the RanSac method.

The main changes concern : -The error tolerance for establishing datum/model compatibility t (Algorithm 1) -The maximum number of attempts to find a consensus set m -The lower bound on the size of an acceptable consensus set T The principle of the algorithm is the following one: 1. -2 normal vectors are randomly selected (N1, N2) 2. -Calculation of the vectorial product between the 2 normal vectors, the vectorial product of N1 with N2 defines the normal to the plane (Vp) 3. -The whole set of normal vectors is browsed to find other normals coplanar with the two first ones, with the condition that Vp is orthogonal to them with a tolerance t (Algorithm1) 4. -After m iterations the plane that aggregated the biggest number of normal vectors is kept as a vanishing plane 5. -The normal vectors classified in a plane are withdrawn from the main set 6. -The whole previous process is repeated until there are only 4 normal vectors left in the set. So that here, T = 4 7. -In each set of normal vectors selected, the best plane is computed by a classical least-squares adjustment.

The number of samples cannot be calculated conventionally with the probabilistic method given by Fischler and R. C. Bolles (Fischler and Bolles, 1987) because no prior knowledge is available. It is impossible to say in advance what percentage of points is good and what percentage is not. A normal vector belonging to a plane is a "inlier" for its plane and will be considered an outlier for another plane, it is therefore impossible to decide on a strict probabilistic basis. In an empiric way, it has been found satisfactory to take for m the half of the number of segments of the image. So the choice of m adapts himself to the image, and therefore the process remains entirely automatic.

ERROR PROPAGATION ON THE VANISHING POINT

After the step of the classification of the planes, it is important to be able to quantify, by a covariance analysis, the uncertainty of the vector of the vanishing point. The impact of the uncertainty of the normal vectors calculated on the vector of the vanishing point is thus computed. The normal vectors are the observations. It is therefore necessary to adjust the best plane fitting these normal vectors while taking in account their uncertainties. The process of least squares used here, using the covariance matrix of the observations, is the method of Helmert-Gauss ((Helmert, 1872) [START_REF] Cooper | Control surveys in civil engineering[END_REF]). It is supposed that all error on the normal vectors (N) have a Gaussian distribution. The cost function here is the square of the distance between the extremities of the normal vectors and the plane that passes through the origin: a, b, c are the components of the extremity of each vector of the normal to the plane. The aim is to estimate and to calculate the covariance matrix on a, b and c. The method of Helmert Gauss is formulated here in the following way:

A is the Jacobian matrix of with respect to the unknown parameters (a, b and c), B is the Jacobian matrix of with respect to the unknown observations (Nx , Ny and Nz). The resolution of this system is performed by the minimization of , v being the vector of the observations. The resolution of the system uses the multiplier of Lagrange in an iterative way. For each iteration the vectors x and v are calculated. The convergence is obtained when the values don't evolve beyond 1 pixel. Once obtained the covariance matrix on the vector of the calculated vanishing point, one propagates it on its intersection with the image plane. Thus an assessment of the uncertainty on the localization of the vanishing point in the image plane is got.

RESULTS AND ASSESSMENT

In order to show some examples of results of the classification of the different directions of lines in an image, 2 examples in the images data base have been chosen. The first image is the Louvre pyramid in Paris (Figure 5), chosen as it contains a lot of families of directions of parallel lines. In this image the vectors of vanishing points have the same colour as the extremities of normal vectors belonging to the same plane. A second image (Figure 6) shows that the algorithm is robust and its parameters tune themselves automatically with regard to the quality of detection of the segments. The possible methods of assessment for the algorithms of automatic vanishing points extraction are to be sorted. One can certainly value their efficiency in terms of numbers of vanishing points found, but no reference method exists to find the correct value except by visual manner. Besides, to find a conclusive statistical value is not that obvious: it is quite simple to bias the results, either by doing the tests on images of neighbouring geometries (which induces a bias on the random character of the sample), either by choosing too simple images (all results are then favourable), or unusually complex, such of unclassical buildings where all lines are curves (and the results will be abnormally bad). In the present case, the assessment was based on about one hundred relatively varied images (Figure 7), for which the number of vanishing points found (ranging from 0 to 5 according to the visual inspection of images) was satisfactory. An important set of images, which served to deepened tests, has been put on line on the site http://mahzad.kalantari.free-.fr/Recherches.htm . One can also try to measure in an independent and more precise way the coordinates that one should find for the vanishing points, and then two options are possible: to work on synthesis images, or to use a photo-theodolite. This last one permits to get images with a calibrated camera, and the orientation of the optical axis is measured with an extreme precision thanks to the theodolite that is included in the mount of the camera. The validation has been led in the following way: eight successive images have been taken of a given scene, presenting a building whose three vanishing points correspond to orthogonal directions, with different angles of view. The angles of the phototheodolite have been recorded every time. The eight images have been processed in order to extract the vanishing points in an automatic way. The angles Á and •, computed from (Patias and Petsa,1993) using the values of focal length and position of the optic center as results of a calibration, are bound directly to the angles H and V measured on the theodolite. Indeed, unlike the situation that prevails for the automatic recovery of vanishing points using the present algorithms, it is necessary, when exploiting the vanishing points in order to compute these angles, to use the intrinsic parameters. The three angles computed for each image follow the classical denomination in photogrammetry (!, Á and •), and their geometrical significance is presented in the Figure 8. The discrepancies of orientation between the successive images have been compared depending on whether one exploits the vanishing points or the values measured by the theodolite (considered as a reference). The standard deviations for the discrepancies on the angles Á and • have been obtained on a set of eight images of the same scene under different angles of view (Table 1). The value obtained is 0.2 mrd for the component Á, that corresponds to the planimetric orientation of the facades observed. It is equivalent, to a distance around 10 m (mean distance from the camera to the facades) to an error of two mm, which is therefore a satisfactory order of magnitude. One deduces that the algorithm workswith a good precision on this type of scene. The discrepancies are higher, around 10 mrd, with the vertical angles extracted from the vanishing points •, as compared with H angles. This is due to the high uncertainty on the intersection of the images of the vertical lines of the object, the consequence being to provide a comparably poor value for .

Figure 1 .

 1 Figure 1.Geometry of the image and the vanishing plane, whose normal in O intersects the image on the vanishing point V

Figure 2 .

 2 Figure 2. Illustration of the geometric features used in the algorithms: the set of planes built on the segments of the image and containing O intersect in V, whatever the position of O

Figure 3 .

 3 Figure 3. Definition of the reference system on the image

Figure 5 .

 5 Figure 5. Typical result on an image of the Louvre pyramide, Paris. The different families of parallel lines are displayed along with the end of the extremities of normal vectors issued from O relative to each corresponding segment, with the same colour.

Figure 6 .

 6 Figure 6. Automatic classification of segments

Deriche, R., Vaillant, R. and Faugeras, O., 1992. From noisy edges points to 3d reconstruction of a scene : A robust approach and its uncertainty analysis. World Scientific Series in Machine Perception and Artificial Intelligence 2, pp. 71-79. Fischler, M. A. and Bolles, R. C., 1987. Random sample consensus:A paradigm for model fitting with applications to image analysis and automated cartography. In: M. A. Fischler

CONCLUSION

A new algorithm, based on very simple geometric considerations and strongly helped by a very efficient RanSactype methodology of classification, allows for a fully automatic extraction of the vanishing points available in a given image. This method and O. Firschein (eds), Readings in Computer Vision: Issues, Problems, Principles, and Paradigms, Kaufmann, Los Altos, CA., pp. 726-740. Helmert, F. R., 1872. Die Ausgleichsrechnung nach der Methode der kleinsten Quadrate. Teubner, Leipzig. is fully independent of any preliminary knowledge of the intrinsic parameters of the optics used, which facilitates its use in a wide set of practical situations. The computation time required for classical urban scenes ranges around one second with a classical PC, which is adapted to most possible applications. On another hand, a new methodology for testing the algorithms has also been developed, based upon the use of a photo-theodolite. This has provided a very efficient way to check their accuracy, which appears as satisfying. Generally speaking, this method may find some vanishing points in excess (artifacts of lines that intersect without being the images of parallel lines), but at least finds all the true vanishing points.

Lutton, E., Maitre, H. and Lopez-Krahe, J., 1994


Barbara Pascal

Vincent Mauduit

Nelly Pustelnik

Patrice Abry

Scale-free Texture Segmentation: Expert Feature-based versus Deep Learning strategies

Keywords: Deep learning, CNN, Texture, Segmentation, Fractal, Total variation, Wavelets

Texture segmentation constitutes a central task in image processing, classically based on two-step procedures consisting first in computing hand-crafted features devised from a priori expert knowledge and second in combining them into clustering algorithms. Deep learning approaches can be seen as merging these two steps into a single one with both discovering features and performing segmentation. Using fractal textures, often seen as relevant models in real-world applications, the present work compares a recently devised texture segmentation algorithm incorporating expert-driven scale-free features estimation into a Joint TV optimization framework against convolutional neural network architectures. From realistic synthetic textures, comparisons are drawn not only for segmentation performance, but also with respect to computational costs, architecture complexities and robustness against departures between training and testing datasets.

I. INTRODUCTION

Context. Automated image segmentation constitutes a crucial task in image processing and computer vision, for many different purposes ranging from medical imaging [START_REF] Trullo | Joint segmentation of multiple thoracic organs in ct images with two collaborative deep architectures[END_REF] to autonomous driving [START_REF] Song | Lateral driving assistance using optical flow and scene analysis[END_REF]. In the last decade, the tremendous increase in computational/storage capabilities has triggered a massive use of deep learning for segmentation, because deep convolutional neural networks have the potential to discover and aggregate relevant information from large scale shapes to fine scale structures. Recently, numerous real-world applications, related to biological tissues [START_REF] Jennane | Fractal analysis of bone x-ray tomographic microscopy projections[END_REF], [START_REF] Marin | Mammographic evidence of microenvironment changes in tumorous breasts[END_REF], geological samples [START_REF] Angelo | Integrated seismic texture segmentation and clustering analysis to improved delineation of reservoir geometry[END_REF], satellite images [START_REF] Roux | A wavelet-based method for multifractal image analysis. III. applications to high-resolution satellite images of cloud structure[END_REF],. . . drove the focus specifically on texture segmentation, at the core of the present work. Texture segmentation however differs drastically from object detection. Indeed, textures are mostly characterized by small scale statistical properties rather than by geometry and large scale structures. Related works. For years, texture segmentation was performed via a classical two-step procedure: First, prior knowledge or expert choice driven features are computed (e.g., Gabor, gradients, differences of oriented Gaussians,. . . ) ; Second, these features are combined via a clustering algorithm (cf. [START_REF] Andrearczyk | Texture segmentation with fully convolutional networks[END_REF] for a state-of-the-art review). Recently, research focus has been on combining these two steps into a single one to improve interface detection and thus segmentation performance. This has been first envisaged by retaining hand-crafted features but modifying classical frameworks (e.g., spectral clustering based on multiple local cues in [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF], nonnegative matrix factorization in [START_REF] Yuan | Factorization-based texture segmentation[END_REF] or combining estimation and detection into a convex optimization formulation [START_REF] Pascal | Joint estimation of local variance and local regularity for texture segmentation. application to multiphase flow characterization[END_REF]). Recently, Deep Learning renewed this topic, jointly performing feature selection [START_REF] Donahue | Caffe: Convolutional architecture for fast feature embedding[END_REF] as well as segmentation, first for semantic segmentation [START_REF] Sermanet | Overfeat: Integrated recognition, localization and detection using convolutional networks[END_REF], [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF], [START_REF] Hariharan | Hypercolumns for object segmentation and fine-grained localization[END_REF], [START_REF] Luc | Semantic segmentation using adversarial networks[END_REF], rapidly followed by texture segmentation [START_REF] Cimpoi | Deep filter banks for texture recognition, description, and segmentation[END_REF], [START_REF] Andrearczyk | Texture segmentation with fully convolutional networks[END_REF], [START_REF] Ustyuzhaninov | One-shot texture segmentation[END_REF]. Amongst others, fractal textures, also referred to as scale-free because their statistics are not controlled by particular scales, are considered Work supported by CBP (Blaise Pascal Center) with the use of SIDUS (Single Instance Distributing Universal System) implemented by E. Quemener. relevant models for and massively used in numerous real-world applications very different in nature (cf. e.g., [START_REF] Pentland | Fractal-based description of natural scenes[END_REF], [START_REF] Hnat | Scale-free texture of the fast solar wind[END_REF], [START_REF] Wendt | Wavelet leader multifractal analysis for texture classification[END_REF], [START_REF] Keller | Texture description and segmentation through fractal geometry[END_REF], [START_REF] Chaudhuri | Texture segmentation using fractal dimension[END_REF], [START_REF] Pascal | Joint estimation of local variance and local regularity for texture segmentation. application to multiphase flow characterization[END_REF]), ranging from physics of rough materials [START_REF] Roux | A wavelet-based method for multifractal image analysis. III. applications to high-resolution satellite images of cloud structure[END_REF] to art investigations [START_REF] Abry | When Van Gogh meets Mandelbrot: Multifractal classification of painting's texture[END_REF]. Interested readers are referred to several other contributions for further examples [START_REF] Keller | Texture description and segmentation through fractal geometry[END_REF], [START_REF] Kaplan | Extended fractal analysis for texture classification and segmentation[END_REF], [START_REF] Riaz | Texture classification using rotation-and scale-invariant gabor texture features[END_REF] and to a recent work [START_REF] Xu | Fractal dimension invariant filtering and its CNN-based implementation[END_REF], implementing a fractal dimension-based contour detector as a Convolutional Neural Network (CNN). Goal, contributions and outline. Making use of mixtures of synthetic fractal textures (defined in Section II-A), the goal of the present work is to compare expert-knowledge driven texture segmentation exploiting a priori chosen scale-free features (such as local regularity) within an advanced unsupervised joint TV based optimization framework [START_REF] Pascal | Joint estimation of local variance and local regularity for texture segmentation. application to multiphase flow characterization[END_REF] (see Section II-B) versus blind or non informed supervised CNN based Deep Learning (DL) approaches (described in Section III) inspired from those in [START_REF] Andrearczyk | Texture segmentation with fully convolutional networks[END_REF]. Comparisons are reported in Section IV not only in terms of absolute segmentation performance but also in terms of architecture complexities (three DL architectures are devised and compared), of learning complexities (computational costs and hyper-parameter tuning issues are discussed) as well as of robustness with respect to both the size of the training set and departures between training and testing sets, a realistic issue in realworld applications, where one cannot always fully control the natural variability encountered in data.

II. UNSUPERVISED LEARNINGFOR FRACTAL TEXTURES

A. Fractal textures

Fractal textures. A fractal texture consists of a stationary Gaussian field X(z), whose covariance structure is fully defined by its variance Σ 2 and a fractal parameter H (cf. e.g. [START_REF] Wendt | Wavelet leaders and bootstrap for multifractal analysis of images[END_REF] for detailed definition). Textured images correspond to discretization of the field X on a pixel grid, denoted X = (Xn)n∈Ω ∈ R N 1 ×N 2 . Piecewise fractal textures. For the present study, we consider piecewise fractal textures, consisting of a mixture of Q fractal textures (cf. Fig. 1, top left, for an example with Q = 3), each characterized by parameters Σ 2 q , Hq , q = {1, . . . , Q}. Local multiscale analysis. It has been well-documented that fractal textures can be well analyzed using multiscale transforms, such as wavelet coefficients and wavelet leaders (cf. e.g., [START_REF] Wendt | Wavelet leaders and bootstrap for multifractal analysis of images[END_REF]). Wavelet leaders L j,k at scale j and location 2 j k, are defined as a supremum of all wavelet coefficients located into a small spatial neighborhood at all finer scales [START_REF] Wendt | Wavelet leaders and bootstrap for multifractal analysis of images[END_REF]. It has been shown then that textures can be well-characterized by a local multiscale analysis as

L j,k = 2 j →0 vn2 jhn , n = 2 j k, (1)
with v 2 n ∝ σ 2 n (a local variance) and where hn quantifies the regularity of the texture locally around location n. For homogeneous fractal textures, local variance and local regularities are constant over the texture: hn ≡ H and σ 2 n ≡ Σ 2 . For piecewise fractal textures, hn and σ 2 n are expected to form piecewise constant regions where ∀n ∈ Ωq, σ 2 n ≡ Σ 2 q and hn ≡ Hq.

B. Joint TV segmentation

Goal. The aim is to recover the underlying partition Ω = ∪ Q q=1 Ωq (cf. bottom-right on Fig. 1) of the image domain Ω from estimations of hn and σ 2 n performed locally (i.e. pixel-wise) using wavelet leaders. Segmentation scheme. To achieve that goal, we designed in [START_REF] Pascal | Joint estimation of local variance and local regularity for texture segmentation. application to multiphase flow characterization[END_REF] an objective function combining a least square data fidelity term, fitting the mathematical model (1), and total variation penalization, enforcing piecewise constancy, leading to the optimization problem:

v, h := arg min v,h j 2 j=j 1 v + jh -j (X) 2 2 + λPα(v, h) (2)
where j (X) = log 2 (Lj) ∈ R N 1 ×N 2 , j = {j1, . . . , j2} are the non-decimated log-leaders and Pα = TV(v) + αTV(h) (cf. [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] for TV definition). The regularization parameters, λ > 0 and α > 0, have to be tuned to reach the best segmentation accuracy.

To perform the minimization (2), we have recourse to an accelerated proximal primal-dual algorithm with convergence guarantees to the global minimizer v, h [START_REF] Pascal | Nonsmooth convex joint estimation of local regularity and local variance for fractal texture segmentation[END_REF]. The segmentation M ROF (i.e. a label map of same size as the analyzed image) is then obtained by thresholding h as proposed in [START_REF] Cai | Multiclass segmentation by iterated rof thresholding[END_REF], i.e. M ROF = TQ( h).

III. SUPERVISED FULLY CNN

Fully convolutional neural networks (FCNN). -We followed the work of Andrearczyk et al. [START_REF] Andrearczyk | Texture segmentation with fully convolutional networks[END_REF] in which the state-of-the-art semantic segmentation FCN-8s Network [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF], is tailored to texture segmentation. Mathematically, the FCNN takes entry X ( ) (256 × 256 image) and outputs 256×256 label map M ( ) ≡ R θ (X ( ) ), where θ denotes the weights of the network and R θ encapsulates the successive (convolution, activation, . . .) operations. The two major ingredients tuned to texture segmentation are, first, the skip connections with shallow layers (enabling to keep precise localization), and second, making use of a large number of filters at each layer to capture as much as can be of the rich properties of textures. The state-of-the-art FCNN in [START_REF] Andrearczyk | Texture segmentation with fully convolutional networks[END_REF] reaches a size of 8 • 10 7 weights for textured images of size N1 = N2 = 256. Networks designed here. Starting from the state-of-the-art FCNN described above, we further propose two successively simplified versions, whose sizes reduce respectively to 2 • 10 6 and 4 • 10 5 weights. The global architecture of the three networks is thus similar and they only differ in complexity (number of weights). Hence, in the following, they are referred to by their number of weights to emphasize the complexity hierarchy. The proposed networks are composed of (i) Several successive convolutional blocks, or layers, (BLOCKS 1, 2 on Fig. 1) combining a sequence of convolutions with ReLU activations. Layers are terminated by a polling operation, inducing resolution decrease enabling to deal with more global features as one gets deeper into the network. (ii) These layers are followed by a fully convolutional block composed of one convolutional layer having a large number of filters and wider kernels of size 5, one convolution getting to targeted final depth Q and one transposed convolution, upsampling the feature maps (BLOCK F on Fig. 1).

(iii) This is followed by skip connections between the output of shallow blocks (e.g. BLOCK 1 on Fig. 1) and the output of BLOCK F, performing successive appropriate upsampling to produce a label map with same resolution N1 × N2 as input images. (iv) Finally, a softmax activation is applied to produce decision, that is one class label per pixel.

(a) Two classes

(Q = 2) (b) Three classes (Q = 3) (c) Four classes (Q = 4)
Fig. 2: Texture configurations Σ 2 q , Hq , q = 1, . . . , Q .

These four steps are summarized in Fig. 1 for the 4 • 10 5 w. network. For all networks, the default filter size is p = 3, the default stride is s = 1 and a (2, 2) pooling is applied at the end of each block.

The three explored architectures are • 4 • 10 5 w. net: This proposed network is composed of 2 convolutional blocks. More specifically, BLOCK 1 is composed of two layers, each having 32 filters and BLOCK 2 has two layers, each having 64 filters. The fully convolutional layer in BLOCK F contains 256 filters. One skip connection combine information from BLOCK 1 and BLOCK F using Q filters.

• 2 • 10 6 w. net: This proposed network is composed of 3 blocks, BLOCK 1 with two layers, each having 32 filters, BLOCK 2 with two layers, each having 64 filters and BLOCK 3 with three layers each having 128 filters. The fully convolutional layer in BLOCK F contains 512 filters. Two skip connections connect BLOCK 1 and BLOCK 2 to the final block, each having Q filters.

• 8 • 10 7 w. net: Andrearczyk et al. [START_REF] Andrearczyk | Texture segmentation with fully convolutional networks[END_REF] network is composed of 4 convolutional blocks. BLOCK 1 is composed of two layers, with 64 filters, BLOCK 2 is composed of two layers, with 128 filters, BLOCK 3 is composed of three layers with 256 filters and BLOCK 4 has three layers, each having 512 filters. The fully convolutional layer in BLOCK F contains 4096 filters. Three skip connections connect BLOCK 1, BLOCK 2 and BLOCK 3 to final block. Supervised learning -From a database X ( ) , M ( ) L =0 an optimization strategy estimates the optimal weights θ of the network as:

θ = argmin θ d( M ( ) , R θ (X ( ) ) (3)
where R θ summarized the successive steps (i)-(iv), composing the neural network. Finally the supervised segmentation is obtained as

M ( ) FCN = R θ (X).

IV. EXPERIMENTS ON SYNTHETIC TEXTURES

A. Experimental settings

Learning dataset. As sketched in Fig. 2, we investigate Q-class segmentation with Q ∈ {2, 3, 4}. Each training image X ( ) is built from a partition

Ω ( ) = ∪ Q q=1 Ω ( ) q , with q = q ⇒ Ω ( ) q ∩ Ω ( ) q = ∅ (4)
randomly generated, and from the fractal descriptors of the Q textures Σ 2 q , Hq q=1,...,Q . For each Q, we consider two configurations, i.e., two different sets of characteristic descriptors { Σ 2 q , Hq , q = 1, . . . , Q}. In Config. I (in blue on Fig. 1) textures to be segmented have large ∆H ≡ Hq -H q and small ∆Σ 2 ≡ Σ 2 q -Σ 2 q and conversely in Config. II (in red in Fig. 1). To train neural networks, for Config. I, II and Q = {2, 3, 4}, we generated a data set of 2000 X ( ) , associated with ground truth segmentation M ( ) (label map equivalent to the partition

Ω ( ) = ∪ Q q=1 Ω ( ) q )
. Each texture X ( ) contains samples of each Q textures (i.e. Ω ( ) q = ∅, ∀q, ∀ ). TV-based texture segmentation. Following [START_REF] Pascal | Joint estimation of local variance and local regularity for texture segmentation. application to multiphase flow characterization[END_REF], the wavelet transform of textured images is performed using Daubechies wavelets with two vanishing moments [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. The considered scales range from 2 j 1 = 2 1 to 2 j 2 = 2 3 . For each Q and each configuration, the Joint algorithm is run on one training image X (1) over a 20 × 20 grid of regularization parameters (λ, α), so as to select optimal parameters (in term of segmentation accuracy). Parameters λ and α are then frozen to optimal values and used to segment a set of one hundred test images (built in the same way as training images). FCNN supervised training. The training for supervised learning consists in minimizing the loss d, computed over the entire training dataset (cf. Eq. ( 3)) with respect to the weights θ of the network. Twoclass segmentation is trained minimizing binary cross entropy, while three and four-class use categorical cross entropy. The minimization is performed using ADAM optimizer with AMSGrad strategy [START_REF] Reddi | On the convergence of adam and beyond[END_REF], learning rate 2 10

B. Results and Discussion

Convergence. Fig. 3a reports, for Config. I and with Q = 2, the evolution of accuracy for the three networks, trained on the training dataset of 2000 textures. Fig. 3a clearly shows that FCNN trained on the full dataset of 2000 textures do converge to a stable maximal accuracy after the use of 30 epochs. Equivalent plots and conclusions are obtained for other configurations. Performance. Tab. I reports segmentation performance after convergence, defined as segmentation accuracies (percentage of wellclassified pixels). For Q = 2, Tab. I shows that both FCNN and Joint TV optimization Algorithm achieve comparably high accuracy on Config I and II. For Q = 3, FCNN perform as well as for Q = 2 for both configurations, while Joint TV optimization Algorithm only maintains competitive accuracy on Config II. For Q = 4, FCNN performance suffer only from a mild decrease, while those of Joint TV optimization Algorithm significantly degrade. FCNN with different complexities (different depth and width) show similar performance despite a decrease by 200 of the number of parameters from the largest state-of-the-art FCN-8s to the smallest FCNN tested here. Robustness w.r.t. training set size S. To investigate the impact of the training set size, learning was conducted a second time using a small subset of the learning dataset described in Sec IV-A of only S = 20 textures. Fig. 3a shows that convergence requires a significantly larger number of epochs. Yet, Fig. 3b shows, that convergence occurs for The estimation is performed using the formula proposed in Eq. ( 5).

Interface accuracy. FCNN lead to irregular interfaces (as emphasized in [START_REF] Andrearczyk | Texture segmentation with fully convolutional networks[END_REF] and illustrated on Fig. 4). Yet, for applications requiring precise interface length measurements, such as [START_REF] Pascal | Joint estimation of local variance and local regularity for texture segmentation. application to multiphase flow characterization[END_REF], it is crucial to recover accurate and smooth contours. Tab. III compares relative errors on total interface length (i.e., the number of pixels at the border of two different regions) for Q = 2-class segmentation, and shows that joint TV optimization systematically and significantly outperforms FCNN. Moreover, when trained in one configuration and tested on another, FCNN perform particularly bad at measuring interfaces.

Train V. CONCLUSION This work showed that supervised FCNN compares favorably against unsupervised hand-crafted feature Joint TV optimization for fractal texture segmentation, even with reduced training data sets, and reduced complexities (compared to state-of-the-art architectures), at the cost though of much larger computation and memory resources. However, FCNN did not show robustness in training/testing datasets mismatch (even those with largest complexities) and yield poor estimation of interface lengths, thus still making their use debatable in real-world applications. Following recent works in deep learning literature a natural option would be to combine both approaches so as to exploit the advantages of each, as LISTA [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF] did.

Fig. 1 :

 1 Fig. 1: Fully convolutional neural network for texture segmentation. Sketch of FCNN with skip connexions (simplest architecture, 4 10 5 weights).

Fig. 3 :

 3 Fig. 3: Evolution of accuracy along training phase. Accuracy of different FCNN trained on Config. I, for two-classes segmentation.

• 10 7 Fig. 4 : 7 FCNN 8 •

 7478 Fig. 4: Q = 2-class segmentation (trained and tested on Config. I), showing irregular contours for FCNN segmentation.

TABLE I :

 I Segmentation accuracy (% of well-classified pixels).

TABLE II :

 II Computational cost C for training phase.

  /Test Config. I/I Config. II/II Config. I/II Config. II/I

	Joint TV	14 ± 4	13 ± 2	48 ± 22	16 ± 4
	FCNN 8 • 10 7 33 ± 3	36 ± 2	112 ± 24	68 ± 18
	FCNN 2 • 10 6 37 ± 3	41 ± 3	213 ± 72	77 ± 4
	FCNN 4 • 10 5 24 ± 3	18 ± 2	337 ± 131	62 ± 9

TABLE III :

 III Relative error on interface length for two-class segmentation. FCNN are trained on S = 2000 images. (Similar results are obtained with S = 20.)


S Alighanbari

Florin G S Giri

Lucian Constantin

V I Korobov

S Schiller
email: step.schiller@hhu.de

F L Constantin

& S Schiller

Precise test of quantum electrodynamics and determination of fundamental constants with HD + ions

come   L'archive ouverte pluridisciplinaire

Since the inception of quantum mechanics, the precise understanding of three-body systems has represented a challenging fundamental physics problem. Its detailed study, both theoretical and experimental, is an ongoing effort, with a strong rate of improvement. Different three-body systems (for example, the helium atom, lithium ion, helium-like ions, antiprotonic helium atom and molecular hydrogen ions (MHIs)) provide the opportunity to test our understanding of quantum physics at the highest levels, in particular, the theory of quantum electrodynamics (QED). In doing so, important fundamental constants of physics (such as the Rydberg constant R ∞ , fine-structure constant α, electron mass m e , proton mass m p , deuteron mass m d and antiproton mass) and particular nuclear properties, such as charge radii, electric quadrupole moments and charge-current moments, can be determined.

The MHIs (HD + , and so on) are molecular three-body systems containing two heavy particles and one light particle (electron). The electronic ground state supports hundreds of metastable rotation-vibration levels. A small subset of them have been studied with different experimental techniques and concerning different aspects since the mid-1960s 2 + Η 5- [START_REF] Biesheuvel | Probing QED and fundamental constants through laser spectroscopy of vibrational transitions in HD +[END_REF][START_REF] Jefferts | Hyperfine structure in the molecular ion H 2 +[END_REF] (for an early review, see ref. 3 ). Over the past decade, the MHIs have come into focus because of their relevance for the metrology of the particle masses 4, [START_REF] Schiller | Test of time-dependence of the electron and nuclear masses with ultracold molecules[END_REF][START_REF] Bakalov | The electric quadrupole moment of molecular hydrogen ions and their potential for a molecular ion clock[END_REF][START_REF] Karr | H 2 + and HD + : candidates for a molecular clock[END_REF][START_REF] Schiller | Simplest molecules as candidates for precise optical clocks[END_REF] . These can be determined from rotationvibration spectroscopic data, an approach independent of the established technique of mass spectrometry in ion traps. An additional opportunity is the determination of the Rydberg constant R ∞ and the proton charge radius, independently from the established technique of atomic hydrogen spectroscopy [START_REF] Beyer | The Rydberg constant and proton size from atomic hydrogen[END_REF][START_REF] Fleurbaey | New measurement of the 1S-3S transition frequency of hydrogen: contribution to the proton charge radius puzzle[END_REF][START_REF] Bezginov | A measurement of the atomic hydrogen Lamb shift and the proton charge radius[END_REF] . The precise value of these constants has been called into question in recent years in connection with the 'proton radius puzzle' [START_REF] Antognini | Proton structure from the measurement of 2S-2P Transition Frequencies of Muonic Hydrogen[END_REF] , and therefore alternative and independent approaches for its determination are highly desirable.

The ab initio theory of the MHIs has made enormous progress in precision over the past 20 years [START_REF] Grémaud | Highly accurate calculation of the energy levels of the molecular ion[END_REF][START_REF] Moss | Energies of low-lying vibration-rotation levels of H 2 + and its isotopomers[END_REF][START_REF] Taylor | Variational calculations on the hydrogen molecular ion[END_REF] , reducing the uncertainty by four orders of magnitude. It currently stands at 1.4 × 10 -11 fractionally for the fundamental rotational transition frequency and 7 × 10 -12 for vibrational transitions, ignoring spin-structure effects [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF] . These uncertainties are smaller than the current (CODATA 2018 [START_REF] Tiesinga | Values of fundamental physical constants[END_REF] ) uncertainties of the masses m e , m p and m d , pointing at the potential of MHI spectroscopy for the metrology of fundamental constants. Here we perform precision spectroscopy of the fundamental rotational transition of HD + . Fundamental constants can be derived by comparison of the measured transition frequency f (exp) with the prediction , the uncertainty originating from the CODATA 2018 uncertainties of the fundamental constants is smaller, 4.4 × 10 -12 , and stems from the uncertainties of r p and r d .

Apart from a matching comparison with a 50-year-old radiofrequency (RF) spectroscopy benchmark result on (ref.

2

+ Η 17
), the ab initio theory could not be tested experimentally at a competitive level, owing to lack of suitable experimental methods. Apart from a few exceptions, the spectroscopic resolution in rotational and vibrational molecular ion spectroscopy in general has been limited by Doppler broadening. Although this broadening can be minimized by trapping molecular ions in an RF trap and sympathetically cooling them by atomic ions, their effective temperature remains of the order of 10 mK, leading to Doppler-limited linewidths not lower than 5 × 10 -8 fractionally [START_REF] Bressel | Manipulation of individual hyperfine states in cold trapped molecular ions and application to HD + frequency metrology[END_REF] . Unresolved hyperfine structure increases linewidths again [START_REF] Koelemeij | Vibrational spectroscopy of HD + with 2-ppb accuracy[END_REF][START_REF] Biesheuvel | Probing QED and fundamental constants through laser spectroscopy of vibrational transitions in HD +[END_REF] , posing a roadblock for testing theory at more precise levels.

Only recently, new methods have been introduced that open up the next generation of precision experiments [START_REF] Wolf | Non-destructive state detection for quantum logic spectroscopy of molecular ions[END_REF][START_REF] Chou | Preparation and coherent manipulation of pure quantum states of a single molecular ion[END_REF] . Specifically for rotational spectroscopy, we have shown [START_REF] Alighanbari | Rotational spectroscopy of cold and trapped molecular ions in the Lamb-Dicke regime[END_REF] that sub-Doppler spectroscopy is possible for a radiation propagation direction transverse to the 'long' axis of the molecular ion cluster (trapped ion cluster transverse excitation spectroscopy, TICTES). The small motional amplitude of the ions along the spectroscopy wave propagation direction compared with its wavelength allows reaching the Lamb-Dicke regime. In the first demonstration [START_REF] Alighanbari | Rotational spectroscopy of cold and trapped molecular ions in the Lamb-Dicke regime[END_REF] , a fractional line resolution of 1 × 10 -9 (full-width at half-maximum (FWHM) relative to absolute frequency) was obtained.

Here we improve the resolution of TICTES by more than two orders of magnitude. This enables a detailed direct study of the fundamental rotational transition of HD + , whose hyperfine spectrum and Zeeman splittings are resolved and systematic effects are determined.

Comparison with our improved theory and new analysis method allows us to establish agreement between theory and experiment at the 5 × 10 -11 level (limited by CODATA2018 uncertainties), not only representing the most accurate test of a molecular three-body system so far, but also demonstrating the power of TICTES, a method applicable to a plethora of molecular ions.

The experiment

We performed spectroscopy of the fundamental rotational transition (v, N) = (0, 0) → (v ′ , N ′ ) = (0, 1) at 1.3 THz. v and N are the vibrational and rotational quantum numbers, respectively.

See Extended Data Fig. 1 for the experimental scheme. The fractional population of HD + ions in the lower spectroscopy state (0, 0) is enhanced using rotational laser cooling [START_REF] Schneider | All-optical preparation of molecular ions in the rovibrational ground state[END_REF] . The transition is detected by resonance-enhanced multiphoton dissociation (REMPD) [START_REF] Roth | Ion-neutral chemical reactions between ultracold localized ions and neutral molecules with single-particle resolution[END_REF] . See Extended Data Fig. 2 for typical data. To achieve a spectroscopy wave with narrow linewidth, high frequency stability and high accuracy, a GPS-monitored, hydrogen-maser-referenced terahertz frequency multiplier is used [START_REF] Alighanbari | Rotational spectroscopy of cold and trapped molecular ions in the Lamb-Dicke regime[END_REF][START_REF] Schiller | Ultra-narrow-linewidth continuous-wave THz sources based on multiplier chains[END_REF] . Compared with our previous work [START_REF] Alighanbari | Rotational spectroscopy of cold and trapped molecular ions in the Lamb-Dicke regime[END_REF] , we performed measurements for different magnetic-, electric-and light-field strengths, and minimized the terahertz wave power. These extensive measurements were enabled by improvements in the long-term stability of the apparatus and improved detection schemes.

The HD + molecule has spin structure in both the lower and the upper rotational levels, due to the presence of (1) the intrinsic spins of the electron (s e ), proton (I p ) and deuteron (I d ), and

(2) of the rotational angular momentum N (Fig. 1). For state description, we use the angular [START_REF] Bakalov | High-precision calculation of the hyperfine structure of the HD + ion[END_REF] ), where F is the total angular momentum. The rotational transition encompasses 32 hyperfine components f i in absence of a magnetic field; of these, ten are favoured (strong) (Fig. 1).

momentum coupling scheme G 1 = s e + I p , G 2 = G 1 + I d , F = G 2 + N (ref.

Systematic shifts

For an accurate comparison between theoretical transition frequencies (computed assuming an absence of perturbing fields) and experimental values (measured in presence of such fields), the systematic shifts must be taken into account. We determined them experimentally. The dominant systematic effect is the Zeeman shift. For a nominal RF drive amplitude, we measured the frequency shifts of all considered components as a function of applied magnetic field. The shifts are consistent with the theoretically calculated ones, except for small deviations. We obtained the transition frequencies corresponding to zero magnetic field by extrapolation.

The quadratic Stark shift due to the ion trap's electric field E(t), oscillating at comparatively low (RF) frequency and leading to a mean-square value <E(t) 2 >, is a second shift, of lower magnitude. For a nominal magnetic field, we measured the frequency shifts of all considered components for a set of trap RF amplitudes. All shifts were found to increase with amplitude, with values in the range of 0.5 to 1.2 kHz kV -2 . We determined the frequencies corresponding to zero RF-field amplitude by extrapolation. For additional information, see Methods and Extended Data Fig. 3.

computed [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF] including all relativistic and radiative corrections up to the relative order α 5 and partially including contributions of the order α 6 (Table 2). The value is updated from the value reported in ref.

( theor ) spin-avg f [START_REF] Alighanbari | Rotational spectroscopy of cold and trapped molecular ions in the Lamb-Dicke regime[END_REF] by using CODATA 2018 [START_REF] Tiesinga | Values of fundamental physical constants[END_REF] updates of the Rydberg constant, the particle masses (in atomic mass units, u), the proton charge radius and the deuteron charge radius. The coefficients E 4 , E ′ 4 and E 5 , E ′ 5 describe the dominant s e .I p and s e .I d interaction, respectively, and have been calculated with high theoretical precision, including all corrections of the order α 2 E F /h and the leading corrections of the order α 3 E F /h, where E F ≈ h(1.4 GHz) is the Fermi contact energy for the hyperfine splitting in atomic hydrogen and h is Planck's constant [START_REF] Korobov | Theoretical hyperfine structure of the molecular hydrogen ion at the 1 ppm level[END_REF] .

The fractional theoretical uncertainties of these spin Hamiltonian coefficients are of the order α 3 ; they are estimated as ε F = 1 × 10 -6 . Furthermore, the theory errors are expected to be nearly equal: ΔE 4 (theor) ≈ ΔE ′ 4 (theor) and ΔE 5 (theor) ≈ ΔE ′ 5 (theor) (Methods).

The other spin coefficients,

E ′ 1 , E ′ 2 , E ′ 3 , E ′ 6 , E ′ 7 , E ′ 8
, and E ′ 9 , have been obtained within the Breit-Pauli approximation. We computed them using our most precise non-relativistic nonadiabatic molecular variational wave functions (Methods, Extended Data Table 1). The omitted terms are of the relative order α 2 . References [START_REF] Korobov | Theoretical hyperfine structure of the molecular hydrogen ion at the 1 ppm level[END_REF][START_REF] Menasian | High-resolution study of (1,1/2,1/2)-(1,1/2,3/2) HFS transition in H 2 + . Bull[END_REF] lead us to estimate a common fractional theory uncertainty equal to α 2 = ε 0 ≈ 5 × 10 -5 .

To determine the impact of the theory uncertainty of a particular Hamiltonian coefficient on a particular spin frequency, we consider the quantity ,

( ) theor ' ' , k k i E Δ γ , i k γ ′ is the derivative ( ) ' ' 9 ' 1 ' spin, ' , ,...,
k i k i E E E E ∂ ∂ = γ
, and similarly for the coefficients relevant to the lower spin level.

The γ values are reported in Extended Data Table 1. Assuming equal theory errors for the pairs (E 4 , E ′ 4 ) and (E 5 , E ′ 5 ), we conservatively estimate the total theory uncertainty of the spinfrequency contribution with the following expression

( ) ( ) ∑ ∑ + - = 9 , 8 , 7 , 6 , 3 , 2 , 1 ' ' , 0 5 , 4 , ' ' , F theor , spin k k i k k i k k i i E E E f u γ ε γ γ ε
The minus sign in the first sum embodies the assumption of equal errors and correlation,

Comparison between theory and experiment

Table 1 presents the comparison between the theory and experimental data of the individual hyperfine components of the rotational transition. We find agreement for all lines, within the combined uncertainties of theory and experiment. The agreement is most stringent for line 16, and it is limited by the prediction's total uncertainty , or about 1.5 × 10 (theory) 16

( ) 0.21kHz u f ≈ -10
fractionally. The agreement is far less stringent than the roughly ten times lower experimental uncertainty would allow. The precise experimental value can therefore serve as a benchmark for tests of future improved spin-structure calculations.

Frequencies related to only the spin structure of the molecule can be obtained from rotational frequency differences , where the spin-averaged frequency is cancelled. All deviations between experiment and theory are smaller than 0.42 kHz in magnitude and are well within the theory uncertainties (CODATA 2018 uncertainties are not relevant here). The most stringent theory-experiment agreement is found for Δf , s p i n ,

Δ i j i j i j f f f f f = -= -s p i n ,
21,19 , within the roughly 0.7-kHz theory uncertainty, but ten times less stringent than the experimental uncertainty would allow.

In view of the relatively large uncertainties for above, we introduce a novel way of comparing experiment with theory, using composite frequencies defined as

(theor) spin,i f c i i i f b f = ∑ , with
appropriate weights b i . We aim to find composite frequencies with small theory uncertainty, and therefore must suppress the contribution of the spin energies' uncertainties without suppressing the spin-averaged energies that give rise to f spin-avg . The latter requirement is satisfied by imposing the 'normalization' condition

1 i i b = ∑ , so that f c = f spin-avg + f spin,c , with spin,c spin, i i i f b f = ∑
. The former requirement is implemented by finding the composite frequency that minimizes the theory uncertainty. We use a conservative measure of theory uncertainty that does not assume any relationship between the theory errors of (E 4 , E ′ 4 ,) and of (E 5 , E ′ 5 ,):

( ) ( ) ∑ ∑ ∑ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ γ + γ = k i k k , i i i ' k ' k , i i c , ε E b E b f u k theor spin . The solution {b i } is found numerically (see 'Composite frequencies' in Methods),
, with negligible uncertainty

(theor) spin,c ({ }) 934.635kHz i f b = (theor) spin,c ( ) 0.001kHz u f = .
We note that this approach for eliminating the spin-energy-related uncertainty is complementary to the more general method recently proposed by some of us in ref.

36

, where the composite frequency is equal to f spin-avg .

From the experimental composite frequency, we deduce the experimental spin-averaged

frequency (exp) (exp) (theor) spin-avg c spin,c 11 exp r ({ }) ({ })
1,314,925,752.910( 17) kHz ( 1.3 10

) i i f f b f b u - = - = = × (2)
The theory uncertainty (via ) is negligible and is therefore not indicated.

( theor ) spin,c f

QED test and determination of fundamental constants

A comparison of equations ( 1) and (2) indicates that our experiment and theory achieve a successful test of three-body physics with a combined fractional uncertainty of 4.8 × 10 -11 (0.064 kHz), limited by CODATA 2018 uncertainties. Comparing the total uncertainty of with the QED contributions listed in Table 2, we see that it is close to the QED contribution of highest calculated relative order of α ( theor ) spin-avg f 6 , f (6) ≈ 0.070(18) kHz. Therefore, more specifically, our experiment furnishes a test of QED at the relative order of α 6 . According to theory, the contributions to f (2) stemming from the finite proton root-mean-square charge radius r p and the deuteron charge radius r d with their CODATA2018 uncertainties are -0.644(3) kHz and -4.120(3) kHz, respectively. The sum of these contributions is put in evidence by our experiment-theory comparison, with a fractional uncertainty of 1.4%.

Our experiment-theory agreement is obtained when including in the hyperfine structure calculation the contribution of the deuteron quadrupole moment Q d , quantified by the coefficient . This contribution is observed here in an MHI for the first time. From the measured hyperfine structure we can extract, independently of any QED contributions, a value for

Q d ' Q E ∝ 9 d with 1.5% fractional uncertainty (Methods).
The experiment-theory agreement can also be used to set improved limits to the hypothetical existence of a spin-averaged fifth force between a proton and a deuteron (Fig. 3, Methods). Compared with previous bounds from MHI spectroscopy, the improvement is a factor of 21 or more for force ranges λ > 1 Å.

We can obtain the combination R ∞ m e /μ pd of fundamental constants from any of the measured rotational frequencies ( )

theor i F ( ) ( ) ( ) ( ) ( ) ( ) ( )
- ∞ × = = = + . u . , cF f m m m R (3)
where the third uncertainty is due to the proton and deuteron radius uncertainties. The value is in agreement with the CODATA 2018 value of 8,966.205 150 41(41) m -1 (u r = 4.6 × 10 -11 ) (Fig.

Conclusion

The performance of the recently introduced TICTES technique for rotational spectroscopy has been improved by more than two orders in both resolution and accuracy, reaching a fractional FWHM linewidth of 3 × 10 -12 and a fractional uncertainty of 1.3 × 10 -11 . This vastly higher performance compared with traditional techniques can be of general relevance to the field of precision molecular physics.

Precise measurements of several rotational hyperfine components of HD + and suppression of the impact of the limited accuracy of the ab initio theory of the spin structure allowed us to establish agreement between experiment and theory at the 5 × 10 -11 level, limited by uncertainties of the CODATA 2018 fundamental constants. To the best of our knowledge, this represents the most accurate test of a molecular physics prediction to date and also provides the most accurate experiment-theory comparison for any three-body quantum system 2, [START_REF] Pastor | Absolute frequency measurements of the 2 3 S 1 → 2 3 P 0,1,2 atomic helium transitions around 1083 nm[END_REF][START_REF] Hori | Buffer-gas cooling of antiprotonic helium to 1.5 to 1.7 K, and antiprotonto-electron mass ratio[END_REF][START_REF] Rengelink | Precision spectroscopy of helium in a magic wavelength optical dipole trap[END_REF] . Specifically, we confirmed the combination of the QED contributions of α 5 relative order, of the proton finite size contribution and of the deuteron finite size contribution, with uncertainty equal to 0.7% of the total contribution. A strongly improved upper bound for a new force between a proton and a deuteron was set.

Spin-energy differences were experimentally determined with three orders smaller uncertainty than previously [START_REF] Bressel | Manipulation of individual hyperfine states in cold trapped molecular ions and application to HD + frequency metrology[END_REF] . The best (effective) line resolution for spin energy is one order higher and the accuracy is 30 times higher than the benchmark experiment on the spin structure of , which has stood unchallenged for 50 years. The spin-energy predictions were confirmed within the uncertainties of the theory predictions, the smallest uncertainty being 0.7 kHz. As the experimental uncertainties are much lower, the obtained spin-energy data offer new benchmark values for future improved ab initio theory of the spin structure.

2

+ Η
We deduced the combinations

( 1 e p d ) 1 R m m m - - ∞ +
and m e /m p of fundamental constants with 2.0 × 10 -11 fractional uncertainty, 2.4 and 3.0 times smaller, respectively, than the CODATA 2018 uncertainties. The proton mass in atomic mass units was determined with the same uncertainty as in CODATA 2018. Interestingly, for the first time,fundamental constants have been determined with competitive uncertainty making use of the rotational motion of a physical system.

Our result also provides independent evidence of the correctness of some of the most precise measurements in atomic and particle physics: Rydberg constant determination via hydrogen spectroscopy, electron mass determination via the bound-electron g-factor, and proton mass and deuteron mass determination via cyclotron motion. Our measurement on a three-body quantum system thus provides an independent link between these one-and two-body systems.

The substantial changes introduced in the CODATA 2018 adjustments of the fundamental constants are confirmed. In particular, the predicted HD + transition frequency is shifted by 0.063 kHz when the CODATA 2014 proton root-mean-square charge radius and Rydberg constant are replaced by the values deduced from the muonic hydrogen experiment (as in CODATA 2018). Our experimental frequency is consistent with the prediction based on these most recent values, within the combined experiment-theory uncertainty of 0.025 kHz.

Beyond the present results, our work has important implications for the near future. First, we suppose that in the spectroscopy of vibrational transitions a similar absolute systematic uncertainty can be achieved as in rotational spectroscopy, because the systematic shifts will not increase substantially with transition frequency. Indeed, the shifts depend on the size of the coefficients of appropriate hamiltonians, and these coefficients do not vary substantially between the levels. If an optical spectroscopic technique with spectral resolution at the roughly 10-Hz level becomes available, total experimental uncertainties on the 10 -13 to 10 -14 level could come into reach. Second, our composite frequency approach obviates the need for a more precise spinstructure theory, both for rotational and vibrational transitions. Therefore, more precise QED calculations of the spin-averaged rotational and vibrational frequencies are both sufficient and well worth pursuing. If this challenging programme is successful, the precision of fundamental constants derived from HD + spectroscopy will further improve. Specifically, the combination of rotational and vibrational spectroscopy results and ab initio theory will eventually allow the determination of the fundamental constants R ∞ , m e /μ pd , r p and r d independently rather than in ( , , , , ) -5

E v N G G F ′ ′ ′ ′ ′ spin, spin 1 2 spin 1 2 ( , , , , ) ( , , , , ) i hf
E v N G G F E v N G G F ′ ′ ′ ′ ′ = - 1 = 0, ΔG 2 =
Fig. 3 | Exclusion plot (95% confidence limit) for a Yukawa-type interaction between a proton and a deuteron, deduced from spectroscopy of MHIs.

The parameter space above the lines is excluded. The assumed interaction is

V 5 (R) = βN 1 N 2 exp(-R/λ)/R
, where R is the proton-deuteron distance, λ is the interaction range, N 1 = 1 and N 2 = 2 are the nuclear mass numbers, and β is the interaction strength. Green lines, this work (full green, numerical; dashed green, analytical, equation (4) in Methods); red line, ref.

14

; blue line, ref. [START_REF] Koelemeij | Vibrational spectroscopy of HD + with 2-ppb accuracy[END_REF] ; orange line, ref. [START_REF] Bressel | Manipulation of individual hyperfine states in cold trapped molecular ions and application to HD + frequency metrology[END_REF] . For comparison, the black lines show the limits for the interaction between the antiproton and the helium-4 nucleus, obtained from two different transitions [START_REF] Hori | Two-photon laser spectroscopy of antiprotonic helium and the antiprotonto-electron mass ratio[END_REF] . See Methods for details. The values were calculated using CODATA 2018 values of the fundamental constants. The main contribution f (0) is of order cR ∞ (m e /μ pd ). Recoil corrections (due to finite masses of nuclei) are included fully at the order α 2 ; the leading recoil corrections proportional to m e /m p or m e /m d are included at the order α 3 . Contributions due to the finite size of the nuclei are included in the f (2)

Table 1 | Experimental rotational frequencies, and comparison with theoretical ab initio frequencies

Line i G 1 G 2 F→G ′ 1 G ′ 2 F ′ ( ) exp i f u ( ) ( ) exp i f ( ) theor i f ( )
term [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF] . The one-loop contribution from μ + -μ -vacuum polarization is included in f (3) . The estimated fractional theory uncertainty of the spin-averaged frequency is u r = 1.4 × 10 -11

METHODS

Experimental procedure

We simultaneously trapped Be + and HD + ions in a linear RF trap driven at 14.16 MHz (Extended Data Fig. 1). The distance between the trap centre and the RF electrodes was 4.3 mm. For translational cooling of the molecular ions, we laser-cooled the atomic ions with a laser at 313 nm and the HD + ions were sympathetically cooled via electrostatic interactions with the Be + ions. We estimated the ion secular temperature as about 30 mK. Typically, roughly 10 2 HD + ions were trapped together with about 2 × 10 3 Be + ions. The number of trapped HD + ions affects the spectral resolution of the rotational transitions, since the Lamb-Dicke regime can only be reached when the ions' displacements in the transverse direction are much smaller than the transition wavelength.

Black-body radiation populates the excited rotational levels of the ground vibrational state until a thermal equilibrium population is reached. We counteracted this by pumping the HD + population into the ground rovibrational state using two lasers. They drive the (0, 2) → (1, 1) and (0, 1) → (2, 0) transitions, and the spontaneous decay from the respective excited states eventually transfers a large fraction of the HD + ions in the rovibrational ground state. A quantum cascade laser at 5.48 μm excited the former transition, and a distributed feedback laser at 2.7 μm excited the latter transition.

After rotational cooling, the terahertz radiation was turned on to drive a transition between specific Zeeman components of a specific hyperfine rotational transition. The terahertz wave intensity was controlled with a half-wave plate, a linear polarizer and via the synthesizer output level. A 1.4-μm laser selectively excited molecules from the (0, 1) level to the (4, 0) level.

Molecules in this level were rapidly dissociated by a 266-nm laser.

The spectroscopy scheme relies on the ability to determine the relative decrease of the number of trapped HD + ions. Resonant excitation of the HD + ions' radial secular motion with an auxiliary a.c. electric field couples to the Be + ion ensemble, heating it and causing a change in atomic fluorescence. This fluorescence change is approximately proportional to the number of trapped HD + ions. Applying the secular excitation before and after the REMPD and calculating the ratio of average fluorescence levels provides the fractional decrease of the number of HD + ions. See Extended Data Fig. 2.

As the REMPD process removes HD + ions from the trap, repeated loadings are necessary. With one loading of Be + , approximately 40 loadings of HD + were performed. For each HD + loading, typically five spectroscopy cycles were performed. Each cycle lasted 60 s and provided one data point.

The magnetic field was BB 0 ≈ 45 μT, directed along the trap axis, except during rotational spectroscopy/REMPD, when the field was changed to B ≈ 30 μT or lower, oriented perpendicular to the trap axis and parallel to the terahertz radiation wave vector (Extended Data Fig. 1). The magnitude and direction of the magnetic field were controlled by three pairs of magnetic coils outside the vacuum chamber.

Owing to the complicated statistics of the ion detection process, we assigned one-half of the FWHM of a line as the statistical uncertainty of a measured transition frequency.

Systematic effects

As a guide to and comparison with the experimental work, the ab initio values for various systematic effects were taken from our previous calculations. Explicit values for the Zeeman effect are given in ref. [START_REF] Bakalov | Magnetic field effects in the transitions of the HD + molecular ion and precision spectroscopy[END_REF] . and for the Stark effect in ref. 48 . The ab initio a.c. polarizabilities at the frequency corresponding to the wavelength 266 nm were computed in ref. [START_REF] Alighanbari | Rotational spectroscopy of cold and trapped molecular ions in the Lamb-Dicke regime[END_REF] .

Trap shift

Several systematic shifts are expected to give rise to a quadratic dependence on RF amplitude.

These include the micromotion-induced Stark shift [START_REF] Berkeland | Minimization of ion micromotion in Paul trap[END_REF] , phase-offset-induced Stark shift [START_REF] Berkeland | Minimization of ion micromotion in Paul trap[END_REF] , and a.c.

Zeeman shift due to an alternating magnetic field at the trap frequency correlated with the electric trap drive.

We therefore measured the dependence of the six lines (including three Zeeman components for line 16 and two Zeeman components for line 19) on the trap RF amplitude. The typical values chosen for the RF amplitude were 150 V, 180 V and 245 V. The precise RF amplitude value for each measurement was determined by measuring the radial secular frequency of Be + . See Extended Data Fig. 3 for an example of the frequency shift when varying the trap's RF field amplitude. Fits, assuming quadratic dependence, furnish the correction to be applied for obtaining each line's extrapolated frequency for zero RF amplitude. The theory of the Stark shift 48 predicts shifts of the same sign (positive) and of similar value for all components considered here. The experimental data are consistent with this prediction.

Zeeman shift

Both the linear and quadratic Zeeman shift coefficients vary substantially among Zeeman components and hyperfine components (compare, for example, lines 16 and 19 in Fig. 2). The frequency splitting of the two Zeeman components 16 ± together with the theoretical linear Zeeman splitting coefficient (7.98 kHz μT -1 (ref. [START_REF] Bakalov | Magnetic field effects in the transitions of the HD + molecular ion and precision spectroscopy[END_REF] )) allows the determination of the (time-and ensemble-averaged) magnetic field affecting the molecular ions. For the data shown in Fig. 2, the nominal magnetic field BB nom = 2.98(3) × 10 T is consistent with the value deduced using spectroscopy of the co-trapped beryllium ions . The observed linewidth of the 16

-5

50

± Zeeman components indicates that the magnetic field is homogeneous to at least 1 part in 30 over the molecule sample.

We measured the frequencies at three different values of magnetic field, for RF amplitudes close to the nominal value of 190 V. Since the RF amplitude varied slightly for the individual measurements, each measured frequency was corrected for the trap shift.

To obtain the B → 0 extrapolated frequency, ( )

exp i f
, for each line, we fitted to the measured line frequencies ( ) ( )

exp i f B the sum of ( ) exp i f
plus a quadratic-in-B and/or linear-in-B dependence, depending on the type of Zeeman component. As an accurate measure of the magnetic field, we used the splitting 16 16

f f -- + .
For m F = 0 → 0 Zeeman components, we assumed a quadratic-in-B dependence. For the two components 19 ± and for the two components 16 ± , we allowed for independent linear-in-B shift coefficients α i,+ , α i,-. For 16 16 , f f

+ - +
, we added to the fit functions the quadratic Zeeman shift predicted by theory. From the fits, we found that the 'positive' and 'negative' shift coefficients of a given line are close: and

.

The input data for the magnetic-field dependence fit are the trap-field-extrapolated line frequencies. The reported uncertainty of each ( )

exp i f
contains both the uncertainty of the magnetic-field extrapolation and the uncertainty due to the trap-field extrapolation.

The magnetic field is produced by three solenoids. They were characterized with a magnetic probe before closing the vacuum chamber. We find the field value deduced from the solenoids' currents agrees with the value deduced from the splitting 16 16 f f --+ , within the experimental uncertainty of the former.

Trap-induced a.c. Zeeman shift

This effect would show up as a variation of the splitting between two Zeeman components with the trap RF amplitude. The 19 ± components were measured at 245 V and 154 V, at the nominal magnetic field. Their frequency difference did not change, indicating a negligible a.c. Zeeman shift.

Light shift due to cooling laser

The 313-nm cooling laser permanently irradiates the ion cluster, including when the terahertz wave is on. Its nominal power is 100 μW and the beam radius is 0.25 mm. We measured the effect of a change of the 313-nm laser intensity on 16 f -. No shift was discernible at the 10-Hz level upon increase of the power by a factor of four.

We computed the scalar, tensor and vector polarizabilities of the rovibrational levels at λ = 313 nm using high-precision variational wavefunctions, similar to ref. 48 , obtaining α s (v = 0, L = 1) = 3.5054, α t (v = 0, L = 1) = -0.955, α s (v = 0, L = 0) = 3.4961 and α t (v = 0, L = 0) = 0, in atomic units. The vector polarizabilities are negligible. The computed light shift is of the order of 0.01 Hz. We therefore set the correction due to the 313-nm wave intensity to zero.

Line pulling

We have no observational evidence that Zeeman components, or micromotion-induced sidebands of other Zeeman components, could affect the measured transitions. The small linewidths of the measured transitions are important in this respect. We did not observe any change of 16 f + , 16 f and 0 16

f at the 10-Hz level upon a 500-Hz change of the trap frequency.

d.c. offsets

For every measurement reported in the manuscript, the HD + ions are located along the symmetry axis of the Be + ion cluster. An offset of 10 V was applied to an electrode to displace the beryllium crystal by about 100 μm from the trap axis along the radial direction. We observed that this offset potential does not have an effect on the position of the HD + ions, as also found in molecular dynamics simulations [START_REF] Alighanbari | Rotational spectroscopy of cold and trapped molecular ions in the Lamb-Dicke regime[END_REF] . We measured the frequency shift of 19 f -caused by this offset potential to be 1(10) Hz. Possible day-to-day variations of the trap compensation voltage are a small fraction of the applied offset. Therefore, the size and uncertainty resulting from these variations are negligible.

Light shift due to the two REMPD lasers

The shift due to the 1.4-μm laser and 266-nm laser waves present during spectroscopy has been determined by performing spectroscopy in a different mode, alternating terahertz irradiation and REMPD laser irradiation. It has been measured for all lines and all Zeeman components discussed here. The shifts are smaller than or equal to 0.039 (17) kHz in absolute value. The measured shifts and their uncertainties are used as corrections.

Other shifts

According to theoretical calculations, the black-body radiation shift 48 and the molecular electric quadrupole shift [START_REF] Bakalov | The electric quadrupole moment of molecular hydrogen ions and their potential for a molecular ion clock[END_REF] can be neglected at the present level of accuracy.

Data analysis

Extrapolation of the measured frequencies to zero magnetic field and zero trap amplitude is done by a standard least-squares method. Standard formulae for the propagation of uncertainties are applied.

Spin coefficients, their uncertainties, and sensitivity of the transition frequencies to the spin coefficients

To allow for an accurate comparison between experiment and ab initio theory, we performed a substantially more accurate computation of the spin-structure coefficients of HD + compared with our earlier work [START_REF] Bakalov | High-precision calculation of the hyperfine structure of the HD + ion[END_REF] . We extended the approach developed in ref. [START_REF] Korobov | Theoretical hyperfine structure of the molecular hydrogen ion at the 1 ppm level[END_REF] and the relevant matrix elements were calculated to ten significant digits. Values of the two spin-structure coefficients for the lower level, E 4 and E 5 , and the nine coefficients for the upper level, E ′ 1 , …, E ′ 9 are reported in the Extended Data Table 1. Using these coefficients in the diagonalization of the spin-structure Hamiltonian of ref. [START_REF] Bakalov | High-precision calculation of the hyperfine structure of the HD + ion[END_REF] , we obtain the spin frequencies f spin,i (Extended Data Table 1). . For a given vibrational level, the rotational dependence of the neglected terms in E 4 and E 5 is nearly zero, because these are contact terms determined by the electronic wave function, which depends very weakly on N. This allows us to assume that the neglected terms in (E 4 , E ′ 4 ) and in (E 5 , E ′ 5 ) are essentially equal, respectively.

Under this assumption, the theory uncertainty of a spin frequency due to these

coefficients k = 4, 5 is set to F ε E γ E γ u k k , i ' k ' k , i k - =
, where γ i,k = -∂f spin,i /∂E k is the derivative of the spin energy of the lower quantum state involved in the transition i with respect to the spin coefficient E k , and γ ′ i,k = -∂f spin,i /∂E ′ k is defined analogously for the upper state. The values of the derivatives are presented in Extended Data Table 1.

The spin Hamiltonian coefficients E 4 ≈ E ′

4 and E 5 ≈ E ′ 5 are similar for the two rotational states, and because the transitions studied here are those between similar spin states, for which , , the spin frequencies are small,

1 1 G G ′ = 2 G G ′ = 2 ' ' i E , E , E , E f 5 4 5 4 spin, <<
and the sensitivities are similar, , i k i k , γ γ ′ ≈ . Therefore, we benefit from important reduction of the theory uncertainties u 4 and u 5 contributed by these four coefficients. Even in the least favourable case, line 14, the uncertainty contribution is less than or equal to u 4 + u 5 ≈ 14 Hz (1 × 10 -11 ), that is, negligible compared with the following contributions.

A second set of coefficients, E ′ 1 , E ′ 6 , E ′ 7 , are one to three orders smaller in magnitude, and have estimated fractional uncertainties of ε 1 ≈ ε 6 ≈ ε 7 ≈ α 2 = ε 0 ≈ 5 × 10 -5 . Their absolute uncertainties, 1.5 kHz to 0.06 kHz, are at a relevant level. They enter the spin-structure frequency uncertainty with contributions

0 ε = ' k k E u .
The fractional uncertainties of the coefficients E ′ 2 , E ′ 3 , E ′ 8 , E ′ 9 are similar to ε 0 , but are not relevant at the present experimental accuracy level because the coefficients themselves are much smaller than the others.

As the details of the theory errors are unknown, the total uncertainty of the spin frequencies is set conservatively as the sum over all u k (instead of the root sum of squares).

The sensitivities γ are obtained by first computing the eigenvalues E spin,i and E ′ spin,i of the Hamiltonian analytically and then computing analytically their derivatives with respect to the individual coefficients E k and E ′ k . These derivatives are then evaluated for the set of current theory values for E k and E ′ k .

Fit of the spin Hamiltonian coefficients

From the six measured transitions, we can derive information about the spin Hamiltonian coefficients and about the true spin-averaged frequency. Under the previous assumption of equal theory errors for (E 4 , E ′ 4 ) and for (E 5 , E ′ 5 ), there are six remaining important quantities (E ′ 1 , E ′ 4 , E ′ 5 , E ′ 6 , E ′ 7 and f spin-avg ), and they can be solved for using a set of equations setting the experimental frequencies equal to the corresponding theoretical frequencies allowing for small deviations from the nominal values. We find , having an uncertainty smaller than the theory uncertainty of about . Furthermore, , and .

The from the nominal values cannot be determined precisely (an aspect that is intrinsic to the favoured transitions), but are consistent with zero.

Composite frequencies

The coefficients of the composite frequency given in the main text are: We consider alternative composite frequencies. One alternative ansatz for finding a composite frequency is to impose the 'insensitivity conditions'

( ) ( ) ' k c i k , i i k c β α E f , γ b E f ∂ ∂ = = ∂ ∂ = ∑ α theor theor 0 0
for a suitable subset {k α , k β } of spin Hamiltonian coefficients. As discussed above, if we assume correlated errors for the pair (E 4 , E ′ 4 ) and (E 5 , E ′ 5 ), then the largest theory uncertainties arise from E ′ 1 , E ′ 6 , E ′ 7 . Four experimentally measured transitions are sufficient to satisfy the three insensitivity conditions for these three coefficients.

The normalization condition is easily imposed in addition. Considering, for example, the lines

' k i ' k , i k , i E f γ Γ ∂ ∂ = = theor
, where i = {14, 16, 21} and k = {1, 6, 7}, is close to zero (about 0.008). This implies that these three transitions are nearly linearly dependent and allow for a composite frequency that nearly satisfies the cancellation conditions (and the normalization condition).

If the correlation assumption is not made, the optimum composite frequency based on lines 14, 16 and 21 yields a comparatively large spin-energy uncertainty of 0.22 kHz. For this reason, in the main text, we determined the composite frequency based on six lines.

A third example is the composite frequency based on the five lines 14, 15, 16, 19 and 20: it yields a theory uncertainty (theor) spin,c ( ) 3H u f ≈ z.

Finally, an example of composite frequency for a vibrational transition is the following. . This is only 3 × 10 -14 relative to the vibrational transition frequency f spin-avg ≈ 58.6 THz.

Fifth force bound

Given the present results, the 95% confidence limit to the strength of the fifth force, β max (λ), is approximately given by Here, ΔY(λ) is obtained numerically from perturbation theory as the difference of the expectation value of R -1 exp(-R/λ) in the two rotational states, where R is the internuclear separation in atomic units, and N 1 , N 2 were defined previously.

We have also obtained an analytical approximation for the rotational transition case where R e is the equilibrium distance, and E rot = f rot /2cR ∞ and E vib are the fundamental rotational transition energy and fundamental vibrational transition energy, respectively. They are in atomic units. The previous bounds on β are also discussed in ref. [START_REF] Salumbides | Bounds on fifth forces at the sub-Å length scale[END_REF] .

Electric quadrupole moment of the deuteron

We are also able to deduce a value for the electric quadrupole moment of the deuteron, Q d . The tensor interaction between Q d and the electric field gradient within the HD + molecule [START_REF] Bakalov | High-precision calculation of the hyperfine structure of the HD + ion[END_REF] contributes to the hyperfine structure. It is quantified by the spin Hamiltonian coefficient

E ′ 9 = 5.666 kHz Q ∝ d .
The ratio E ′ 9 /Q d is available from our theory with small fractional uncertainty ε 0 ≈ 5 × 10 -5 . The frequencies of the rotational transition components are sensitive to E ′ 9 to varying degrees, quantified by ,9 i γ ′ (see Extended Data). We therefore consider a composite frequency c i i i f a f ′ = ∑ that suppresses the spin-averaged frequency, and thus all QED

k i k , i E f γ ∂ -∂ = theor spin,
to the lower state. The entries for line 19 are decimal representations of rational values (see eq. ( 6) in ref. [START_REF] Bakalov | Magnetic field effects in the transitions of the HD + molecular ion and precision spectroscopy[END_REF] ). Note that because of the tracelessness of

F

  μ pd = m p m d /(m p + m d ) is the reduced nuclear mass, c is the speed of light, and ≈ 0.244591781951(33) theory(11) CODATA2018 is a dimensionless normalized frequency computed ab initio, neglecting the hyperfine interactions.encompasses, besides the dominant non-relativistic (Schrödinger) part, essential relativistic, nuclear-size-related and radiative contributions. The nuclear charge radius values (r d ) are from the CODATA2018 adjustment that took into account the muonic hydrogen spectroscopy results. Whereas the uncertainty of due to theory is 1.

Figure 2 →

 2 Figure 2 shows the measured transitions, in the presence of a small magnetic field. The different linewidths are due to the different terahertz wave intensities used and due to the different transition dipole moments. Line 19 includes the two transitions between states of maximum total angular momentum F and maximum projection quantum number m F , (F = 2, m F = ±2) → (F ′ = 3, = ± 3), denoted by 19 ' F m ± , whose Zeeman shift is purely linear, according to theory 37 . The two components were observed at lower resolution and with unresolved Zeeman splitting in ref. 16 . One Zeeman component (19 -) measured at particularly low intensity exhibited a full linewidth of 4 Hz, or 3 × 10 -12 fractionally, indicating the potential of the experimental technique in the context of mass determination. For line 16, we measured a Zeeman pair m F = ±1 → = ±2 (denoted by 16 F m ′ ± ), split by a linear Zeeman shift and weakly shifted by a

  rotational levels then leads to a strong suppression of the contributions related to the theory errors of E 4 , E ′ 4 , E 5 and E ′ 5 . This results in the spin-frequency uncertainties shown in

∞

  4). It results from atomic hydrogen spectroscopy (providing R ∞ ), hydrogen-like ion spin resonance spectroscopy (m e ) and Penning trap mass spectrometry (m p , m d ). Our result's total uncertainty is smaller by a factor of 2.4 compared to the CODATA2018 value and ranks among the most precise measurements of a fundamental constant combination. Owing to the comparatively small CODATA 2018 uncertainty of R ∞ , our improved uncertainty impacts mostly the mass ratio sum . Combining equation (3) with the CODATA 2018 values of R , m e /u and m d /u yields the proton mass p exp theor CODATA2018 / 1.007 276 466 605(20) (21) (45) m u = in excellent agreement with the recent most precise direct measurement 40 p s / 1.007 276 466 598(16) (29) m u = Taking into account a recent Penning trap measurement of m d /m p (ref. 41 ), we also obtain the proton-to-electron mass ratio .152 673 449 (24) (25) (13) m m = (u r = 2.0 × 10 -11 ) in excellent agreement but approximately two times more accurate than the most precise value, obtained by combining two published measurements in Penning traps 40,42 : m p /m e = 1,836.152673374(78) exp .

Fig. 1 |

 1 Fig. 1 | Energy diagram of the spin structures and favoured transitions.

  0 and ΔF = 0, ±1. The ten favoured transitions are shown by coloured lines. The rotational transition frequency of a particular hyperfine component is f i = f spin-avg + f spin,i , with f spin-avg ≈ 1.314 THz and, for favoured transitions, f spin,i ≈ Ο(10 MHz). The six components measured in this work are shown by bold numbers in the diagram.

Fig. 2 |

 2 Fig. 2 | Hyperfine components of the fundamental rotational transition of HD + at 1.3 THz. The red and blue points indicate the cases of terahertz radiation on and off (background), respectively. Green lines are Lorentzian fits. The Zeeman components are indicated by the expression m F → . The terahertz wave intensity varied and was less than 10 nW mm F m ′ -2 . The zero of the frequency scales are set to coincide with the fitted line maxima or means. At each frequency setting, the red and blue data points are both shown with an offset equal to the value of the blue point. Each error bar represents the standard deviation of the mean. The nominal magnetic field is BB nom ≈ 3.0 × 10 T and the trap RF amplitude is approximately 190 V.

Fig. 4 |

 4 Fig. 4 | Comparison of results of this work with literature values.

  018 kHz). The impact of the fundamental constants' uncertainties is given in the text. The change in the value of f (theor) spin-avg f (0) from CODATA 2014 to CODATA 2018 has contributions of -0.041 kHz from the Rydberg constant adjustment and 0.213 kHz from the particle masses adjustments. The change in the value of f (2) due to the proton and deuteron charge radii adjustments is 0.104 kHz.

b

  12 = 0.086 372 0…, b 14 = 0.145 634 8…, b 16 = 0.251 611 1…, b 19 = 0.244 279 2…, b 20 = 0.132 807 4…, b 21 = 0.139 295 5…

14, 16 ,k

 16 19 and 21, the resulting uncertainty from hyperfine theory is u( ) ≈ 2 Hz, much smaller than the uncertainty of the spin-averaged frequency . Thus, the composite frequency has a substantially reduced theory uncertainty compared with those of the individual hyperfine transitions. f is then also numerically close to f spin-avg , . With more available transitions we can impose additional conditions. composite frequency is as follows. As in the main text, we consider a composite frequency that minimizes the spin-coefficients-related uncertainty. If we assume correlated E errors, the linear combination of only three lines, f c = b 14 f 14 + b 16 f 16 + (1 -b 14b 16 )f 21 , yields an uncertainty of 3 Hz (2.4 × 10 -12 ). As in the first alternative, this uncertainty is also much smaller than . The coefficients are b 0814…, b 16 = 0.615… and . Such optimal solutions exist independently of the concrete values of the estimated theory uncertainties of the E coefficients: if the assumed fractional uncertainties ε are doubled, a solution is obtained whose theory uncertainty is correspondingly larger, 6 Hz. The relationship between the solution and the cancellation conditions is that the determinant of the sensitivity matrix (theory) c f ( )

For

  the transition (v = 0, N = 0) → (v ′ = 1, N ′ = 1) the six lines 14, 15, 16, 19, 20 and 21 yield a composite frequency with theory uncertainty

  d i = (2F(i) + 1)/36 and = (2F′(i) + 1)/36 are the degeneracies of the respective spin states, and the sum is over the ten favoured transitions i = 12, …, 21.

Table 1

 1

						presents the experimental transition frequencies ( ) exp i f	(corrected for the
	systematic shifts) and their uncertainties. The uncertainties result from the number of frequency
	measurements, which were taken at different RF drive settings and different magnetic-field
	settings, and the statistical uncertainties of the frequency measurements. The lowest
	experimental uncertainty is achieved for line 16, u( ( ) exp 16 f	) = 0.017 kHz (fractional uncertainty u r
	= 1.3 × 10 -11 ). This represents the best performance level of the TICTES technique as currently
	implemented.
	Theory
	For a compelling comparison between theory and the experimental data, highly precise
	theoretical predictions and qualified estimates of their uncertainties are essential. The ab initio
	transition frequency	f	i	(theor	)	of each hyperfine component i is the sum of two contributions,
	f	(theor) spin-avg	+	f	(theor) spin,i	. The dominant contribution is
									f	(theor) spin-avg	=	1,314,925,752.896(18)	theory	(61)	CODATA2018	kHz

Table 1 (

 1 column 6). They dominate the total uncertainty of the transition frequencies .

	f	i	(	theor	)

  because their spin-structure theory uncertainty is suppressed to a negligible level. Furthermore, we note that the ab initio calculation is performed assuming trial values for m e /m p and m e /m d , and naturally yields the rotational frequencies in atomic units (independent of

	exp the highest precision is obtained by instead choosing the composite frequency f i f and the corresponding ab initio value (theor i f c or the spin-. However, ) a.u.. From these, we compute the scaled, dimensionless values ( ) 4 theor 10 998 1 -× ≈ ... . f n , i ( ) ( ) ( ) 1 theor, e pd theor n i i f m μ F = a.u.. These have an important dependence on r p and r d . The dependence on other fundamental constants is weak, compared to their uncertainties, the largest being ( ) ( ) 3 pd e theor 10 4 ln ln -× ≈ ∂ ∂ μ m F i . Because of this smallness, it is consistent to results in averaged frequency, Rydberg constant value), use the CODATA 2018 values of the fundamental constants in the computation of . This ( )

Table 2 | Contributions to the ab initio spin-averaged rotational frequency

 2

	f	( theor ) spin-avg

Extended Data Table 1 | Spin Hamiltonian coefficients, spin-structure frequencies and spin-frequency derivatives E

  ′ k (E k ) are the updated coefficients of the spin Hamiltonian 35 of the upper (lower) rotational level, in MHz. are theoretical spin frequencies in MHz. γ are the dimensionless sensitivities of the spin frequencies to the spin Hamiltonian coefficients.

	f	( theor ) spin,i
		γ	' i	k ,	( = theor i f ∂ spin,	)	E ∂	' k	refer
	to the upper state and

(

)

contributions, by imposing . We determine the weight set {a 0 = ∑ i i a i } that maximizes the sensitivity-to-uncertainty ratio

. We find a 12 = -0.216 516 7, a 14 = 0.650 806 8, a 16 = -0.909 898 9, a 19 = -0.973 830 3, a 20 = -0.115 369 0.

From the comparison of and , we then deduce . It is consistent with the reference value

obtained from RF spectroscopy of neutral D 2 and theory [START_REF] Pavanello | Determination of deuteron quadrupole moment from calculations of the electric field gradient in D 2 and HD[END_REF] . The precision is expected to improve with progress in MHI spin-structure theory.