High-accuracy determination of the beam divergence error in free-fall absolute gravimeters

Gravimeters are instruments for measuring the acceleration due to gravity, g. As their name indicates, absolute gravimeters measure the absolute value of g—compared to relative gravimeters that only measure gravity differences. For some applications, absolute gravimeters are preferred to relative gravimeters—even if the absolute value of g is not crucial. This is because absolute gravimeters have a high long-term stability. An example of such an application is the observation of the post glacial rebound [1], where measurements over many years are necessary (see [2] for an overview over gravimetry). State-of-the-art absolute gravimeters follow the principle of free-fall. The most common of such instruments is the model FG5, or the later version FG5X (extended drop length), from the company Microg-LaCoste, Lafayette (CO), USA [3]. In this device a macroscopic test body is put into free-fall in a vacuum vessel. This test body contains a retroreflector that acts as a test mirror in a laser interferometer. The laser interferometer thus tracks the free-fall position of the test body—with respect to an inertially isolated reference retroreflector—during a total free-fall length of 20 cm to 30 cm. Taking simultaneously time data during the free-fall provides time-position data, from which the acceleration can be calculated. Those instruments are inherently free of drifts.

For other applications, like in geodesy, free-fall absolute gravimeters are needed, as the absolute value of g is necessary for the description of the figure of the earth, or the definition of the height. More recently, absolute gravimeters become also more and more important in metrology, or to be more precise, in mass metrology. Kibble balance experiments [4] require the absolute value of g for the realization of the SI unit kilogram. This is because Kibble balances are basically force comparators—comparing the gravitational force on a mass piece with an electro-magnetic force, generated by a voice-coil actuator. For the conversion of force into mass, the accurate value of local g is necessary, at the instance when the experiment is performed. Most accurate Kibble balances reach combined relative uncertainties of about $1\cdot10^{-8}$.

A difficulty arises when we talk about accuracy of absolute gravimeters. While agreements between different gravimeters are verified in international comparisons—regional comparisons, or key comparisons—such an event cannot evaluate the instrument's accuracy, since the true gravity value is unknown. In order to evaluate its accuracy, an uncertainty budget has been elaborated for the FG5, mainly based on a publication of Niebauer et al in 1995 [5]. Since then, various research papers have been published, focusing on possible (in the budget of [5]) under- or overestimated effects. One of such investigations treats the so-called beam divergence of the laser beam of the interferometer (in the literature also named diffraction correction). This effect has been evaluated in [5], based on a reference to [6], to be on the order of −2 μGal (1 μGal = 10 nm s⁻²), with an attributed uncertainty of 0.2 μGal. Later this effect has been re-investigated, e.g. by van Westrum and Niebauer [7] and Robertsson [8]. Van Westrum et al measured the beam waist diameter with a knife edge method. They concluded that the effect is on the order of 1 μGal in most gravimeters. The evaluation of the error was based on the formula given by Dorenwendt and Bönsch [9], considering a perfect Gaussian beam shape, and applied to a balanced interferometer. Robertsson did a theoretical study considering a more realistic, namely unbalanced laser interferometer, as is actually the case in the FG5(X) systems. In his study Robertsson already pointed out the importance whether the interferometer is balanced or not, if uncertainty levels on the order of below 1 μGal are aimed for. Moreover, the smaller the beam waist the smaller the uncertainties must be for the determination of the beam waist, as the divergence error increases with decreasing beam waist. Finally, Robertsson stresses that the beam shape may deviate from a perfect Gaussian, and therefore it might be necessary to include higher order structures.

A deviation from a perfect Gaussian can be characterized by a parameter called M². During the 10th International Comparison of Absolute Gravimeters (ICAG) 2017 Wang et al compared the beam waists of different participants [10]. Wang et al measured the beam waists with a commercial beam profiler (Model CinSquare, from CINOGY Technologies GmbH, [11]), which, besides the beam waist radius, also provides an estimate of M². Beam divergence errors have been estimated ranging from 1 μGal to 5 μGal. The correction applied were based on the formula of Dorenwendt and Bönsch, thus for a balanced interferometer. The authors mention that the main reason for the higher beam divergence errors were badly adjusted collimators. After a re-adjustment, the errors were reduced, resulting in a better agreement between the participants. Unfortunately, no link was presented in this study between the beam divergence error and M².

In 2018 and 2022 also Křen and Palinkas reported on beam divergence errors in absolute gravimeters [12, 13]. The authors calculated the beam divergence error for two different collimators, the original one, that came along with the commercial gravimeter FG5(X) (from Microg-LaCoste), and a commercial collimator (from stock) from the company Thorlabs (model TC25APC-633, [14]). It was shown that, although the beam spot diameter of the Microg-LaCoste collimator was larger than that from the Thorlabs collimator, the beam divergence error of the latter was smaller. This fact sounds counter intuitive, because after the formula given by Dorenwendt and Bönsch, the error should be smaller the bigger the beam diameter. The authors explained this paradox with the worse beam shape quality of the Microg-LaCoste collimator, which can be quantified from M². The measurements were based on beam profiling (home made device) at various distances, in order to obtain the necessary parameters. In [13] the same authors presented an improved setup and analysis method. The estimated divergence errors were −3.8 μGal for the Microg-LaCoste collimator, and −2.8 μGal for the one from Thorlabs.

The research on the beam divergence error shows that the measured gravity value can be erroneous by several microgal. The question arises whether the models that were used for the evaluation of this effect are correct. This cannot be verified when only one collimator, with one specific focal length is used. Another problem is, that a beam profiler cannot acquire the wavefront parameters, which are necessary to evaluate whether the beam can be described by a Gaussian distribution, or if other polynomials are needed. In such a case, the formula of Dorenwendt and Bönsch is not valid any longer.

In this research article we present a new approach to investigate the beam divergence error. Firstly, four different commercial collimators are measured with a Shack–Hartmann wavefront sensor, providing the wavefront curvature and the irradiance distribution. Since each collimator has a different beam waist parameter we expect a different beam divergence error. By means of numerical simulations the relative error for the gravimeter is calculated, considering the exact path lengths of both beams, the reference beam, and the test beam within the laser interferometer. Then, we perform gravity measurements with each of the four collimators, installed into PTB's free-fall absolute gravimeter FG5X-263.

In this section an analytical formula for the detected light power of a split Gaussian beam which travels two distinct paths inside an interferometer is derived. Furthermore, the relevant associated length errors are obtained and explained.

Laser beams are commonly described by Gaussian beams which are solutions of the paraxial wave equation [15, 16]. In complex notation a simple stigmatic Gaussian beam on the detector is given by

$$\begin{equation} E\left(r,z\right) = A \cdot V\left(r,z\right)\exp\left(-\mathrm{i} kz\right), \end{equation} \tag{ 1 }$$

where A is an arbitrary real amplitude we do not need to consider here, k is the wavenumber, r is the radial coordinate on the detector, and z is the travelled distance from the waist position of this beam. The complex amplitude $V(r,z)$ can be expressed with the complex beam parameter $q(z) = z+\mathrm{i} z_{\tiny \textrm{R}}$ as

$$\begin{equation} V\left(r,z\right) = \frac{1}{q\left(z\right)}\exp\left[-\mathrm{i} k\frac{r^2}{2q\left(z\right)} \right], \end{equation} \tag{ 2 }$$

where $z_{\tiny \textrm{R}}: = \pi w_0^2/\lambda$ is the characteristic Rayleigh length, w₀ is the waist radius, and λ is the wavelength of the beam [15, 16]. In figure 1 the origin of the characteristic parameters of a Gaussian beam can be deduced.

Zoom In Zoom Out Reset image size

We assume that the full power of the interfering beams is captured by the detector. This requirement is in compliance with the usual demand that no significant parts of a laser beam are cut anywhere in its path. After they are split the two beams may travel different lengths $z_{1,2}$ before they are recombined. The detected light power may be described as

$$\begin{align} P\left(z_1,z_2\right) &= \int_0^{2\pi} \mathrm{d} \theta \int_0^\infty \mathrm{d} r\,r \left[E\left(r,z_1\right)+E\left(r,z_2\right)\right] \nonumber\\ &\quad \cdot \left[E^\ast\left(r,z_1\right)+E^\ast\left(r,z_2\right)\right], \end{align} \tag{ 3 }$$

where the asterisk denotes the complex conjugate. Following [17] the integrals in equation (3) yield the analytical formula for the light power of the interfering Gaussian beams:

$$\begin{equation} P\left(\Delta z\right) = \frac{\lambda}{z_{\tiny \textrm{R}}}\left\{1+\frac{1}{\sqrt{1+\frac{\Delta z^2}{4 z_{\tiny \textrm{R}}^2}}}\cos\left[\varphi\left(\Delta z\right)\right]\right\}, \end{equation} \tag{ 4 }$$

with $\Delta z = z_1 - z_2$ and

$$\begin{equation} \varphi\left(\Delta z\right) = k \Delta z - \arctan\left(\frac{\Delta z}{2z_{\tiny \textrm{R}}}\right). \end{equation} \tag{ 5 }$$

The fraction with the square root term is responsible for the interference contrast which decreases for larger Δz [17]. As can be seen in equation (5), the phase $\varphi(\Delta z)$ consists of two contributions:

• (i)  
  the naive phase $k\Delta z$;
• (ii)  
  the phase error $-\arctan(\Delta z/2 z_{\tiny \textrm{R}})$.

The naive phase is the dominating term one would get from the interference of plane waves (which do not exist in reality). It is the phase an unaware operator of an interferometer would base his/her evaluation on. The phase error is the small contribution which originates from the more rigorous description (still in paraxial approximation) by Gaussian beams.

Now, based on above considerations, a relative error can be defined as

$$\begin{equation} \varepsilon_\varphi : = -\frac{\arctan\left(\frac{\Delta z}{2z_{\tiny \textrm{R}}}\right)}{k\Delta z}. \end{equation} \tag{ 6 }$$

It quantifies the error as a fraction of the measurement value, the latter based on the naive phase. This error accumulates as soon as $\Delta z \neq 0$. Thus, it is the error one would get when starting a measurement in a balanced interferometer with the possibility to scan very far out.

The derivative of equation (5) with respect to Δz is obtained as

$$\begin{equation} \frac{\mathrm{d} \varphi\left(\Delta z\right)}{\mathrm{d} \left(\Delta z\right)} = k - \frac{1}{2z_{\tiny \textrm{R}}\left(1+\frac{\Delta z^2}{4z_{\tiny \textrm{R}}^2}\right)} = k + \delta k\left(\Delta z\right). \end{equation} \tag{ 7 }$$

It depends mainly on k with the small error δk. Similarly to equation (6) an error can be defined as

$$\begin{align} \varepsilon_k\left(\Delta z\right) : = \frac{\delta k\left(\Delta z\right)}{k} = -\frac{1}{2kz_{\tiny \textrm{R}}\left(1+\frac{\Delta z^2}{4z_{\tiny \textrm{R}}^2}\right)} = -\frac{1}{k^2w_0^2+\frac{\Delta z^2}{w_0^2}}. \nonumber\\ \end{align} \tag{ 8 }$$

This is essentially the local error of the wavenumber. For $\Delta z^2 \lt\lt w_0^2$ the well known relative length error (RLE) of Dorenwendt and Bönsch is obtained [9]:

$$\begin{equation} \lim_{\Delta z \rightarrow 0} \varepsilon_k\left(\Delta z\right) = -\frac{1}{k^2w_0^2} =: \varepsilon_{\mathrm{DB}}. \end{equation} \tag{ 9 }$$

Therefore, equation (9) is an approximation for the case when Δz is negligibly small, i.e. for balanced interferometers and only very small displacements. Nevertheless, this is used as the standard equation for diffraction correction in absolute gravimeters. In figure 2 the three errors defined above are compared for a Gaussian beam with $w_0 = 0.5$ mm and λ = 632.8 nm.

Zoom In Zoom Out Reset image size

It is apparent that the use of $\varepsilon_{\mathrm{DB}}$ is not advisable for unbalanced interferometers and also not for longer scan ranges, while ε_k can be used also for unbalanced interferometers, but also only for relatively short scans. If the scans become much longer, it should be integrated along the respective sections of Δz.

In this section the measurement procedure along with the setups is described. For our investigation we use the following commercial triplet collimators (all from Thorlabs): $\textit{TC25APC-633}$, $\textit{TC18APC-633}$, $\textit{TC12APC-633}$ and $\textit{TC06APC-633}$.

The measurement is performed in various steps. The first one is for determining the wavefront and irradiance distribution of the beam as it exits the collimator. This is done with a Shack-Hartmann sensor (SHS). The second step is the calculation of the relative error of the gravimeter measurement. Therefore, the parameters, as obtained from the measurement with the SHS are used. In a third step, the collimators are installed into a free-fall absolute gravimeter of type FG5X. This is to validate the theoretical model, which is used for the numerical simulations. Agreement between theory and experiment can then be taken as an estimate of an absolute value for the divergence error of each collimator.

3.1. Procedure for determination of the beam parameters

The evaluation software (SHSWorks v11.021.4 SVN3802) of the manufacturer of the SHS (Optocraft SHSCam HR2-130-U3, [18]) is used to obtain the corrected wavefront, i.e. after subtraction of piston and tilt, and the normalized power density, i.e. the irradiance (cf figure 3), of the measured beam. The wavefront data is given in terms of the used wavelength which is throughout all measurements $\lambda = {632.8}\,\textrm{nm}$.

Zoom In Zoom Out Reset image size

The coordinates of the data points are also obtained from the manufacturer's evaluation software. The micro-lens array of the SHS contains a deliberate defect, i.e. a missing lens in the middle of the viewing area. The purpose of this defect is to identify a unique position which is present in all data sets as NaN-value (not a number). Ideally it should be enclosed by valid data points as it is shown in figure 3, because invalid data is also represented by NaNs. We ensured during the measurements, that the NaN defect is always surrounded by valid data (but not in the center of it, where the maximal irradiance can be usually found). This is important because the coordinates of wavefront and irradiance data may differ a little. To realign the data coordinates we developed the following algorithm:

• (i)  
  The position of the centroid of the irradiance data is subtracted from the respective coordinates.
• (ii)  
  The position of the NaN defect is found in the coordinates of the just centered irradiance data.
• (iii)  
  The position of the NaN defect in the wavefront data is subtracted from the respective coordinates and the position of the NaN defect in the centered irradiance data is added.

We denote the aligned coordinates $x,y$, the corresponding wavefront $W(x,y)$ and the irradiance $I(x,y)$. A complex amplitude representing the measured beam can be obtained via

$$\begin{equation} E\left(x,y\right) = \sqrt{I\left(x,y\right)}\exp\left[-2 \pi \!\mathrm{i}\,W\left(x,y\right)\right]. \end{equation} \tag{ 10 }$$

A simplex fit [19] of equation (1) to $E(x,y)$ by equation (10) is then used to obtain the characteristic beam parameters w₀ and z, which are also the used fit parameters. In appendix A we show how to obtain fairly good starting values from wavefront and irradiance. Since the measurement distance $z_{\mathrm{m}}$ to the SHS is known, a synthetic beam with arbitrarily high resolution and window size over $x,y$ can be calculated via equation (1) at the position of the collimator under test $z-z_{\mathrm{m}}$. This position is chosen to average unwrapped phase and amplitude of all reconstructed beams of all measurements with varying $z_{\mathrm{m}}$. From this averaged beam a final fit is done to obtain a final w₀ and z. For the later simulation we only need w₀.

Alternatively, the RLE of each measurement are calculated by use of equations (9.36)–(9.39) of [17] and then averaged. The resulting quantity is named $\varepsilon_{\mathrm{N,m}}$. Then, with equation (9) w₀ can be calculated from $\varepsilon_{\mathrm{N,m}}$. Both methods yield comparable values (see $w_{0,\mathrm{m}}$ and $w_{0,\mathrm{f}}$ in table 2).

3.2. Numerical simulation of the relative acceleration error

In figure 4, a simplified sketch of the Mach-Zehnder interferometer is shown as it is realized in a FG5(X). According to the definitions in section 2, z₁ would be the sum of the dotted green and all continuous orange paths, while z₂ would be the sum of the dashed red and all continuous orange paths. For the difference Δz the common orange paths cancel out. The unbalance between measurement and reference path lengths is quite obvious.

Zoom In Zoom Out Reset image size

Instead of only integrating $k+\delta k(\Delta z)$ directly along the relevant section of Δz to estimate the error of our gravimeter, we decided to perform also a more holistic simulation of the measurement process. This has the advantage that hitherto unknown contributions of the phase retrieval by capturing the zero-crossings of the cosine-function of the detected interference signal are implicitly included.

According to equation (4), the contrast is also a function of Δz. Therefore, the amplitude of the interference signal is bound to change during the fall of the test mass. Due to this change the zero-crossing times may be shifted and the retrieved phase may be modified. Additionally, with constant time resolution the sampling resolution must decrease due to the acceleration. This could also have an impact on the retrieved phase.

From technical information we received directly from Microg-LaCoste and from the technical manual [20] we estimate the relevant path information. Furthermore, we neglect all optical elements, assuming that all apertures are sufficiently large to avoid beam clipping and treat an overall propagation medium as vacuum. Our estimated reference path is assumed to be $z_1 = {150}\,\textrm{mm}$. For the test path we arrive at $z_2 = {1575}\,\textrm{mm}$. This path is extended by the doubled drop length. The path difference used for the signal generation by equation (4) is then

$$\begin{equation} \Delta z = z_1-z_2-2h\left(t\right), \end{equation} \tag{ 11 }$$

with the time t, the height $h(t) = h_{\mathrm{max}} - gt^2/2$, $h_{\mathrm{max}} = {457.3}\,\textrm{mm}$, and g = 9810 mm s⁻². The time interval for t was chosen according to a free-fall length of $h_{\mathrm{fall}} = {304}\,\textrm{mm}$ (as in the experiment), i.e. $t = {0}\,\textrm{s}\;{\ldots}\;{0.248\,953}\,\textrm{s}$. The chosen time resolution is 21 ns. For the simulation of the interference signal we simply employ the term inside the curly bracket of equation (4) and subtract the mean over all its sampling points to enable the detection of zero-crossings. In figure 5 an arbitrary section of a typical simulated interference signal with the found zero-crossings is plotted. The respective search algorithm finds the root of a linear interpolation between the points adjacent to a detected sign change. The found zero-crossings $t_0(j)$, where j is a counted integer index, are mapped to the "measured" heights via $h_{\mathrm{m}}(j) = j\lambda/4$. Then, the following linear equation system is solved for all j as

$$\begin{equation} a_2 t_0^2\left(j\right) + a_1 = h_{\mathrm{m}}\left(j\right). \end{equation} \tag{ 12 }$$

Finally, the simulated relative error of g is obtained directly by

$$\begin{equation} \varepsilon_{\mathrm{g}} = \frac{2a_2-g}{g}. \end{equation} \tag{ 13 }$$

The continuous blue curve in figure 6 shows the quite surprising dependence of the so simulated $\varepsilon_{\mathrm{g}}$ on the waist radius w₀ of a Gaussian beam with wavelength $\lambda = {632.8}\,\textrm{nm}$. As can be seen in the figure, the absolute value of the relative error reaches a maximum and for even smaller beam waist radii it becomes smaller again. Therefore, it appears to be possible to reduce the diffraction error by going to even smaller waist radii. However, as these beams get more divergent one has to be cautious to avoid beam clipping.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The analytical formula used to obtain the red dashed curve is derived in appendix C. It follows from the integration of $k+\delta k(\Delta z)$ along the relevant section of Δz and is given here for convenience:

$$\begin{align} \varepsilon_{\mathrm{g,a}} = \frac{\arctan\left(\frac{z_1-z_2-2h_{\mathrm{max}}}{2z_{\tiny \textrm{R}}}\right)-\arctan\left(\frac{z_1-z_2-2h_{\mathrm{max}}+2h_{\mathrm{fall}}}{2z_{\tiny \textrm{R}}}\right)}{2kh_{\mathrm{fall}}}. \nonumber\\ \end{align} \tag{ 14 }$$

For large Rayleigh lengths $z_{\tiny \textrm{R}}$ (or, equivalently, large beam waists w₀) this equation reduces to equation (9), i.e. the correction term due to the unbalance of the interferometer becomes insignificant. Towards smaller beam waists, however, the correction term becomes dominant and avoids that $\varepsilon_{\mathrm{g}}$ tends towards minus infinity, what equation (9) would suggest (black dotted line in figure 6). A noticeable result of this fact is that there is a maximal error which will not be exceeded. The limits can be better seen, when equation (14) is approximated by a Taylor expansion as

$$\begin{equation} \varepsilon_{\mathrm{g,a}} \approx -\frac{\beta w_0^2}{w_0^4+\alpha^2}\,, \end{equation} \tag{ 15 }$$

where we have defined $\alpha: = \frac{(\Delta z-2h_{\mathrm{max}})\lambda}{2\pi}$ and $\beta: = \frac{\lambda^2}{4\pi^2}$. It is instructive to observe, that it reflects basically equation (8), which is the case of considering the local error only, without integrating over the scan range. This solution is plotted as green dash-dotted line in figure 6. Again, for $w_0\gg\alpha$ we get equation (9). For $\lim_{w_0 \rightarrow 0}$ α becomes dominant and the $\varepsilon_{\mathrm{g,a}}$ becomes smaller and tends to zero. The maximum error occurs approximately where $w_0 = \pm\sqrt{\alpha}$. In our example this is at about $w_0 = {0.485}\,\textrm{mm}$. Remarkably also, it does not depend on g nor on t but on the drop length $h_{\mathrm{fall}}$. The holistic simulation and the analytical formula agree very well, as can be seen as the blue continuous line in the lower graph of figure 6. Therefore, the other systematic effects mentioned above (contrast and sampling resolution) are not really relevant (which is also a result).

3.3. Measurement with the free-fall absolute gravimeter

In order to verify theory versus experiment, the measurement data with $TC25APC$ were taken as a reference with respect to the other collimators. As

$$\begin{equation} \frac{\bar{g}_i-\bar{g}_{25}}{\bar{g}_{25}} = \frac{\Delta \bar{g}_i}{\bar{g}_{25}} = \varepsilon_{g,i}-\varepsilon_{g,25}\,, \end{equation} \tag{ 16 }$$

where the subscript 25 refers to collimator TC25APC, and i is the index for the respective other collimators, a direct comparison between theory and experiment is possible. For the data, obtained from gravity measurements, the corresponding beam waists were taken, as obtained by the SHS measurement and the following simulation. Then, two relative gravity errors were calculated, $\varepsilon_{\mathrm{m}}: = \varepsilon_{\mathrm{g}}(w_{\mathrm{0,m}})$ and $\varepsilon_{\mathrm{f}}: = \varepsilon_{\mathrm{g}}(w_{\mathrm{0,f}})$, as described above. Both results are depicted in figure 7. The blue continuous line shows the ε_g versus beam waist radius w₀ as obtained from simulation, considering a free-fall length of 304 mm. Red circles show the experimentally observed relative gravity errors versus beam waists $w_{\mathrm{0,f}}$, whereas green triangles correspond to beam waists $w_{\mathrm{0,m}}$. Vertical error bars are type A uncertainties from the gravity measurements (5 gravity measurements were taken with each collimator). Horizontal error bars show the uncertainty in the beam waist determination, as obtained from the SHS measurements and subsequent simulation. As already stated above, data with $TC25APC$ are assumed to agree with the theory, and thus are the reference. For the other collimators the relative gravity difference is taken with respect to $TC25APC$.

Zoom In Zoom Out Reset image size

The lower graph in figure 7 shows the residuals (measured minus simulated data) for better visibility. It can be seen that the measured values are slightly higher than the simulation suggests. However, even for collimator TC06APC the deviation from simulation is smaller than $1.5\cdot10^{-9}$, what is equivalent to a deviation of about 1.5 μGal in the gravity measurements. This is a good agreement between theory and experiments.

The numerical values of the results are listed in table 1. The second column shows the mean measured absolute gravity values with the different collimators (listed in column 1). The third column presents the relative gravity differences, with respect to collimator $TC25APC$. Those numbers can be directly compared with the differences of the $\varepsilon_{\mathrm{f}}$ or $\varepsilon_{\mathrm{m}}$ of the collimators, with respect to the $\varepsilon_{\mathrm{f}}$ and $\varepsilon_{\mathrm{m}}$, respectively, of $TC25APC$. The $\varepsilon_{\mathrm{f}}$ and $\varepsilon_{\mathrm{m}}$ were obtained from their respective beam waists (shown in table 2), that were measured with the SHS. The numbers show good agreement.

Table 1. Results of the gravity measurement campaign with different collimators. For the mean gravity $\bar{g}$ a constant value of 981 251 000 μGal is subtracted. For the calculation of the relative gravity difference, the absolute difference is divided by 981 251 830 μGal. Subscript 25 refers to the value of collimator TC25APC.

Collimator	$\bar{g}_i/$  μGal	$\frac{\Delta\bar{g}_i}{\bar{g}_{25}}$ /10−9	$\varepsilon_{\mathrm{f}}$ /10−9	$(\varepsilon_{\mathrm{f}}-\varepsilon_{\mathrm{f,25}})$ /10−9	$\varepsilon_{\mathrm{m}}$ /10−9	$(\varepsilon_{\mathrm{m}}-\varepsilon_{\mathrm{m,25}})$ /10−9
TC25APC	829.86	0	−2.77	0	−2.77	0
TC18APC	828.34	−1.55	−4.41	−1.64	−4.41	−1.64
TC12APC	823.84	−6.13	−10.10	−7.33	−10.10	−7.33
TC06APC	810.32	−19.91	−23.58	−20.80	−24.08	−21.31

Table 2. Comparison of the beam waists and the relative length error when calculated with the formula by Dorenwendt and Bönsch. Nomenclature of subscripts: DS—Thorlabs data sheet.

Collimator	$w_{0,\mathrm{DS}}$ /mm	$\varepsilon_{\mathrm{DB,DS}}$ /10−9	$w_{0,\mathrm{m}}$/mm	$\varepsilon_{\mathrm{DB,m}}$ /10−9	$w_{0,\mathrm{f}}$/mm	$\varepsilon_{\mathrm{DB,f}}$ /10−9
TC25APC	2.335	−1.86	$1.9106\pm 0.0029$	−2.7785	$1.9117\pm 0.0003$	−2.7755
TC18APC	1.685	−3.57	$1.5058\pm 0.0167$	−4.4734	$1.5142\pm 0.0169$	−4.4239
TC12APC	1.12	−8.09	$0.9753\pm 0.0093$	−10.664	$0.9806\pm 0.0094$	−10.548
TC06APC	0.56	−32.34	$0.5062\pm 0.0101$	−39.586	$0.5263\pm 0.0113$	−36.619

For the sake of completeness table 2 shows the RLE when the formula by Dorenwendt and Bönsch (DB) is taken, rather than the correct simulation by means of equation (4). Also the beam waists, as provided in the data sheet of Thorlabs, are shown. It can also be noted that for a larger beam waist, as with, e.g. collimator model $TC25APC$, the difference in the correction, when equation (9) is taken or equation (4), is negligible for current gravimeter systems. For smaller beam waists, however, the difference can be significant, as, e.g. in the case of collimator TC06APC the correction would be wrong by more than $-15\cdot10^{-9} (-14.7$ μGal).

As we have used four different collimators, with different divergence angles, and all measurements agree well with the simulation shown in equation (4), there is an evidence that our theoretical model reflects the true divergence error of the collimators. That means that the absolute bias due to the beam divergence can be estimated and corrected. For the commercial collimator Thorlabs TC25APC, which is used by some users as an alternative to the original FG5 collimator, the relative bias in the gravity measurement becomes $-2.77\cdot10^{-9}$, or −2.72 μGal. If we would have investigated only one collimator, or two with nearly same divergence angles, such a conclusion would not be possible. This is because a possible systematic offset could not be ruled out, while the possibility that there is a common offset for all collimators, is very low. Interestingly, the result suggests that for very small beam waists, the divergence error becomes smaller again. Therefore, a strategy to reduce the bias, and thus the uncertainty that comes along with it, could be using a collimator which produces a small beam waist. This can be advantageous if bigger beam diameters are restricted by certain apertures in the setup that cannot be made larger. However, there is a trade-off somewhere, because a small beam waist produces a large divergence angle. Then again, at a certain distance, the beam diameter could also become too large to pass the aperture.

Another way to reduce the bias would be to use a smaller laser wavelength. For example, a blue laser of wavelength 405 nm would yield a reduction by 40% of the maximum error. Towards larger beam waists, it would reduce even more (about 60%).

It is also worth to note that these results are applicable to other interferometer setups. For example, in Kibble balance experiments the force factor is determined by means of a velocity measurement of a moving coil, with respect to a permanent magnet [4]. This velocity measurement is usually done with a laser interferometer. A beam divergence error thus would lead to an underestimation of the coil velocity, resulting in an overestimation of the mass. Also an underestimation of local g would lead to an overestimation of the mass in such an experiment. Therefore, it is recommended to investigate such possible error sources in Kibble balance experiments.

In our experiment the free-fall length was 30.4 cm. In other gravimeters the drop length might be different. For example, the FG5 has a free-fall length of only 20 cm. This fact must be taken into account, when the error is estimated for another gravimeter. The simulation of the divergence error must be adapted to each individual drop length.

Finally, we would like to underline that, although the collimators that we have used are commercial ones, every collimator needs to be investigated individually. A replacement of the collimator by an "identical" model might lead to an offset in gravity measurements. Moreover, it is strictly necessary to use the very same fiber when measuring with the SHS and when measuring with the gravimeter. The reason is that the obtained beam waists depend delicately on the numerical aperture of the used fiber. Experiments that we have done, suggest that these can result in a relative difference from fiber to fiber of a few parts in 10¹⁰ for the collimator Thorlabs TC25APC. For collimators with smaller beam waists this deviation can even be bigger. Unfortunately, we cannot give a general uncertainty estimate even for the same type of collimator and fiber, since we have no knowledge about the manufacturing tolerances. The best solution would be to use a collimator with a beam waist that produces a negligible diffraction error. However, this is not possible in the FG5X because of resulting beam clipping. The effect of beam clipping would require further investigation.

The uncertainty estimation of the measurement results is divided into two parts. Firstly, the uncertainty in the determination of the beam waist radius w₀, as this parameter is used as a link between theory and experiment. Secondly, the uncertainty in the gravity experiments.

5.1. Determination of the beam waist radius

5.2. Gravity measurements

For the gravity measurement camapaign m = 5 non-concatenated days of measurement for each collimater were evaluated. The contributing measurement uncertainty is the type A uncertainty ($\sigma_{\mathrm{g}} / \sqrt{m}$) of the respective five measurements. Table 3 shows the mean values, along with the estimated type A uncertainties. Moreover, the total number of accepted drops, as a sum of all five runs, is presented.

Table 3. Results of the gravity measurement campaign with different collimators. For the mean gravity $\bar{g}$ a constant value of 981 251 830 μGal is subtracted. n is total number of drops.

Collimator	$\bar{g} /\unicode{x03BC}$Gal	$u_{\mathrm{g}} /\unicode{x03BC}$Gal	$\sigma_{\mathrm{g}} / \unicode{x03BC}$Gal	$u_{\mathrm{g,r}}$/10−9	n
TC25APC	829.86	0.24	0.54	0.24	8577
TC18APC	828.34	0.20	0.45	0.20	6785
TC12APC	823.84	0.20	0.44	0.20	6086
TC06APC	810.32	0.56	1.26	0.57	8470

5.3. Combined standard uncertainty

The relative combined standard uncertainties for the divergence errors of the collimators is estimated as $u_{\mathrm{c,r}} = (u_{\varepsilon_{\mathrm{m}}}^2+u_{\mathrm{g,r}}^2)^{1/2}$. As the relative uncertainty due to determination of the waist radius is two orders of magnitude lower, its contribution is negligible. The final uncertainties are then essentially those of table 4. For the sake of completeness, the results are presented in compact form again in table 5.

Table 4. Results of the relative gravity errors as simulated from the collimator measurements, and their respective uncertainties.

Collimator	$\varepsilon_{\mathrm{f}}$/10−9	$u_{\varepsilon_{\mathrm{f}}}$/10−9	$\varepsilon_{\mathrm{m}}$/10−9	$u_{\varepsilon_{\mathrm{m}}}$/10−9
TC25APC	−2.7709	0.0030	−2.7709	0.0029
TC18APC	−4.4120	0.0169	−4.4121	0.0167
TC12APC	−10.099	0.0094	−10.099	0.0093
TC06APC	−23.575	0.0113	−24.076	0.0101

Table 5. Final results of the collimator investigation. Given are the absolute and the relative bias, along with the absolute and the relative combined standard uncertainties.

Collimator	$\Delta g / \unicode{x03BC}$Gal	$u_{\mathrm{c}} / \unicode{x03BC}$Gal	$\frac{\Delta g_{\mathrm{r}}}{g}/10^{-9}$	$u_{\mathrm{c,r}}/10^{-9}$
TC25APC	−2.72	0.24	−2.77	0.24
TC18APC	−4.32	0.20	−4.41	0.20
TC12APC	−9.90	0.20	−10.10	0.20
TC06APC	−23.60	0.56	−24.08	0.57

The influence of beam divergence due to the collimators in free-fall absolute gravimeters of type FG5X has been investigated. Therefore, four different collimators, with different focal lengths, and thus different beam waists, have been installed into the gravimeter's Mach-Zehnder type interferometer, and the gravity has been measured. Differences between the gravity values were observed for different collimators. A measurement of the beam parameter of the collimated beam, by means of a SHS were employed to simulate the relative acceleration measurement error during the free-fall path of the measurement retroreflector of the interferometer, with respect to the reference retroreflector. The results showed a good agreement between measured and simulated gravity errors. It can be concluded that our collimator of type $TC25APC$ (in data sheet specified beam waist diameter of 2.335 mm), from the company Thorlabs, introduces a relative bias in the gravity measurement of $-2.77(24)\cdot10^{-9}$, or an absolute bias of $-2.72(24)$ μGal. For beam waists of about 0.6 mm the effect can easily amount to about 20 μGal. Moreover, for a theoretical estimation of this effect for small beams waists, the path length difference Δz of an unbalanced interferometer needs to be taken into account. The reason is that the path length difference becomes of the order of the Rayleigh length. The contributing additional phase error becomes substantial, and cannot be neglected any longer.

The findings are important for gravity measurements where the absolute value of g needs to be known to high accuracy. Especially for the realization of the SI unit kilogram by means of a Kibble balance experiment, the outcomes become twofold. First, the gravity value needs to be known to the same degree to which the mass needs to be calibrated. This is, because the Kibble balance is a force comparator (compares weight to electro-magnetic force), and for the conversion of weight to mass g is necessary. Assuming that a gravimeter suffers from this divergence error and is used in a Kibble balance experiment, then a relative bias of the same amount would enter into the mass determination, however with different sign. The mass value, hence, would be too high by (taking our values) 2.77(24) μg kg⁻¹.

Second, Kibble balances employ laser interferometers to measure the velocity of the coil in a voice-coil actuator, in order to determine the so-called force (or geometric) factor of the voice-coil actuator. The relative uncertainty in the determination of this velocity enters to the same degree into the uncertainty of the mass to be calibrated. For an evaluation of a possible beam divergence effect in such a velocity measurement, equation (8) above can be applied.

The discrepancy with the more divergent beams in figure 7 reminds us of the demur by Robertsson [8] that Gaussian beams may not give us the full picture and it may be necessary to include higher-order modes for the beam analysis and the derived diffraction correction. This view is supported by Mana et al [22] where it has been shown that neglecting higher-order contributions can lead to an underestimation of the diffraction correction by approximately 43%. Therefore, we are currently planning to investigate the influence of higher-order modes in the near future and may report corresponding findings in a forthcoming paper.

Last but not least, let us note that also the commonly used higher-order Hermite- or Laguerre-Gauss modes are subject to the paraxial approximation [16]. If the discrepancy we see in figure 7 hints at a departure from the paraxial approximation and cannot sufficiently be captured by the higher-order description we must resort to non-paraxial models [17].

B A thanks Giovanni Mana for fruitful discussions. The authors also thank Microg-LaCoste for providing the technical drawings and the path lengths for the interferometer of the FG5X absolute gravimeter. C R thanks Susanne Quabis for stimulating discussions on interferometry, and Qiyu Wang for motivating preliminary experiments during his stay at PTB.

In the following appendix sections we make heavy use of the well known formulas for the characteristic Gaussian beam parameters as they are stated in [16]. For sake of brevity and readability we do not write out the dependency on z for w(z) and R(z).

To help the simplex search algorithm [19] to succeed we aim to obtain a fair set of start values for w₀ and z. We use the procedure described as follows:

• (i)  
  The wavefront data $W(x,y)$ is multiplied by the wavelength λ to express it in length units.
• (ii)  
  From a Zernike fit [23] of $W(x,y)\cdot\lambda$ the defocus term D is obtained.
• (iii)  
  The curvature radius of a Gaussian beam can be related to D via:

  $$\begin{equation} R = -\frac{r_{\mathrm{m}}^2}{4D}, \end{equation} \tag{ A.1 }$$

  where $r_{\mathrm{m}}$ is the actual radius of the unit circle of the Zernike fit in length units.
• (iv)  
  From the irradiance data the row and the column which contains the maximum is used to obtain beam profiles for x and y, respectively.
• (v)  
  Using $w = \sqrt{-x^2/\ln[A(x)/A_{\mathrm{\max}}]}$, where $A(x) = \sqrt{I(x)}$, I(x) being the irradiance profile for x, and $A_{\mathrm{\max}}$ is the maximum of A(x), the actual beam width at the measurement position is obtained. The same is done for y and the mean value of the x- and y-value is used further.
• (vi)  
  A lengthy derivation (see appendix B) yields the following equation for the distance from the waist:

  $$\begin{equation} z\left(R,w\right) = \frac{R}{\frac{R^2\lambda^2}{\pi^2w^4}+1}. \end{equation} \tag{ A.2 }$$
• (vii)  
  If the resulting z is not a finite number then we found the unique undefined position z = 0 and w₀ is the already obtained w. This can only happen if D is exactly zero which practically never happens with measured data.
• (viii)  
  If the resulting z is a finite number we first calculate the Rayleigh length via $z_{\tiny \textrm{R}} = \sqrt{[R-z]z}$.
• (ix)  
  Finally, we obtain $w_0 = \sqrt{z_{\tiny \textrm{R}}\lambda/\pi}$.

We set $w_0^2 = z_{\tiny \textrm{R}}\lambda/\pi$ into

$$\begin{equation} w^2 = w_0^2\left(1+z^2/z_{\tiny \textrm{R}}^2\right) \end{equation} \tag{ B.1 }$$

and obtain

$$\begin{equation} w^2 = \frac{\lambda}{\pi z_{\tiny \textrm{R}}}\left(z_{\textrm{R}}^2+z^2\right). \end{equation} \tag{ B.2 }$$

Solving equation (B.2) for $z_{\tiny \textrm{R}}$ yields two solutions:

$$\begin{equation} {z_{\tiny \textrm{R}}}_{1,2} = \frac{\pi w^2}{2\lambda} \pm \sqrt{\frac{\pi^2w^4}{4\lambda^2}-z^2}. \end{equation} \tag{ B.3 }$$

We set this into

$$\begin{equation} z_{\tiny \textrm{R}}^2 = \left(R-z\right)z. \end{equation} \tag{ B.4 }$$

The intermediate result

$$\begin{equation} Rz-\frac{\pi^2w^4}{2\lambda^2} = \pm\sqrt{\frac{\pi^2w^4}{\lambda^2}\left(\frac{\pi^2w^4}{\lambda^2}-z^2\right)} \end{equation} \tag{ B.5 }$$

is squared and transformed into

$$\begin{equation} \frac{R^2z^2\lambda^2}{\pi^2w^4}+z^2 = Rz. \end{equation} \tag{ B.6 }$$

To proceed further we must divide equation (B.6) by z which means, from here on we exclude the solution z = 0 (see appendix A). After division by z the final result equation (A.2) can be obtained after factoring out z and subsequent division by the bracket term.

From $h(t) = h_{\mathrm{max}} - gt^2/2$ we solve for g and obtain

$$\begin{equation} g = \frac{2}{t^2}\left[h_{\mathrm{max}}-h\left(t\right)\right]. \end{equation} \tag{ C.1 }$$

The apparent value of g, i.e. $\tilde{g}$, is obtained by the evaluation of the optical phase along the doubled path (back and forth):

$$\begin{align} \tilde{g} &= \frac{1}{kt^2}\int_{\Delta z_1}^{\Delta z_2}\mathrm{d}\left(\Delta z\right)\;k+\delta k\left(\Delta z\right) \nonumber\\ &= \frac{1}{kt^2}\int_{\Delta z_1}^{\Delta z_2}\mathrm{d}\left(\Delta z\right)\;\frac{\mathrm{d}\varphi\left(\Delta z\right)}{\mathrm{d} \left(\Delta z\right)}. \end{align} \tag{ C.2 }$$

The integration is straightforward. We get

$$\begin{equation} \tilde{g} = \frac{1}{kt^2}\left[\varphi\left(\Delta z_2\right)-\varphi\left(\Delta z_1\right)\right]. \end{equation} \tag{ C.3 }$$

After setting equation (5) into equation (C.3) and rearranging we arrive at

$$\begin{align} \tilde{g} = \frac{2}{t^2}\left[\frac{\Delta z_2-\Delta z_1}{2}+\frac{\arctan\left(\frac{\Delta z_1}{2z_{\tiny \textrm{R}}}\right)-\arctan\left(\frac{\Delta z_2}{2z_{\tiny \textrm{R}}}\right)}{2k}\right]. \nonumber\\ \end{align} \tag{ C.4 }$$

According to equation (11) the limit values for Δz are $\Delta z_1 = z_1-z_2-2h_{\mathrm{max}}$ and $\Delta z_2 = z_1-z_2-2h(t)$. Setting these into equation (C.4) we obtain

$$\begin{align}\begin{array}{l} \tilde g = \frac{2}{{{t^2}}}\left\{ {\left[ {{h_{{\textrm{max}}}} - h\left( t \right)} \right]} \right.\\ \quad \,\,\,\,\left. { + \frac{{\arctan \left( {\frac{{{z_1} - {z_2} - 2{h_{{\textrm{max}}}}}}{{2{z_{\textrm{R}}}}}} \right) - \arctan \left( {\frac{{{z_1} - {z_2} - 2h\left( t \right)}}{{2{z_{\textrm{R}}}}}} \right)}}{{2k}}} \right\}. \end{array} \end{align} \tag{ C.5 }$$

The acceleration error can now be given by

$$\begin{align} \varepsilon_{\mathrm{g,a}}: = \frac{\tilde{g}-g}{g} = \frac{\arctan\left(\frac{z_1-z_2-2h_{\mathrm{max}}}{2z_{\tiny \textrm{R}}}\right)-\arctan\left(\frac{z_1-z_2-2h\left(t\right)}{2z_{\tiny \textrm{R}}}\right)}{2k\left[h_{\mathrm{max}}-h\left(t\right)\right]} \nonumber\\ \end{align} \tag{ C.6 }$$

Setting $h(t) = h_{\mathrm{max}}-h_{\mathrm{fall}}$ we finally obtain equation (14).