Since the first direct detection of gravitational waves (GWs) from a coalescing binary black hole system in 2015 [1] , the rate of detections by the LIGO 7 [2] and Virgo [3] observatories has increased [4] , yielding insight into properties, formation processes, and populations estimates of GW sources. As the Japanese KAGRA [5] detector comes on-line and is eventually joined by a third LIGO detector located in India [6] , and as the sensitivities of the currently-operating detectors continue to improve, the rate of detections by the global network is expected to increase to several per day [7] . To maximize the scientific benefit of these detections, accurate and precise calibration of the detectors is essential. In 2009, L. Lindblom estimated that calibration accuracy of 0.5 % or better would be required to optimally extract the information encoded in the signals [8] . Subsequent analyses have also highlighted the importance of reducing calibration uncertainties [9, 10] . Accurate determination of the distance to GW sources requires low overall network calibration uncertainty. Also, the relative calibration accuracy between detectors in the global network plays an important role in sky localization of sources, enabling improved sky maps for follow-up observations by electromagnetic observatories [11] .
In 2017, the LIGO and Virgo collaborations, together with over seventy electromagnetic observatories, reported the multi-messenger observation of a binary neutron star inspiral [12] . These observations enabled an independent measurement of the Hubble parameter [13] , albeit with insufficient precision to resolve the tension between the results reported by the Planck [14] and the SHoES [15] collaborations. This result was constrained by the signal-to-noise ratio (SNR) rather than calibration accuracy. As the number and SNR of GW detections increase, measuring the Hubble parameter with the ∼ 1 % accuracy that is needed to resolve the tension in the current estimates is expected within the next decade. It will require detector calibration accuracy of 1 % or better [16] .
Current gravitational-wave interferometers use systems that are referred to as photon calibrators (Pcals) to produce periodic fiducial displacements of suspended interferometer mirrors via photon radiation pressure [17, 18, 19, 20] . A schematic diagram of a Pcal system installed at an end station of one of the LIGO interferometers is shown in figure 1. Power-modulated auxiliary laser beams reflecting from the interferometer mirrors cause differential relative arm length variations that mimic the variations induced by gravitational waves. The magnitudes of the induced length variations are proportional to the amplitudes of the modulated laser power. Thus the interferometer displacement accuracy depends directly on the accuracy of the calibration of the Pcal laser power sensors. The scheme for transferring laser power calibration from the US National Institute of Standards and Technology (NIST) to the Pcal power sensors, and for coordinating the relative calibration of all of the detectors in the global GW network, is also shown schematically in figure 1.
When power-modulated Pcal laser beams reflect from the surface of the mirror suspended at the end of one of the arms of a gravitational wave interferometer, referred to as an end test mass (ETM), the periodic force exerted on the optic is given by where θ is the angle of incidence of the Pcal beams on the ETM surface, c is the speed of light, P(ω) is the amplitude of the modulated laser power reflected from the test mass, and ω = 2πf is the angular frequency of the power modulation. Conversion of the periodic Pcal-induced forces to displacements requires the force-to-length transfer function of the suspended mirror, S(ω) in units of m/N. S(ω) for a LIGO 40 kg ETM is plotted in the upper-left panel of figure 2 together with the transfer function for a free mass, S(ω) = −1/(M ω 2 ), where M is the mass of the suspended optic. As shown in the upper-right panel of figure 2, at frequencies above 20 Hz S(ω) is well approximated (within 0.1 %) by the response of a free-mass. The Pcal forces can also induce unintended rotation of the test mass due to power imbalance between the two Pcal beams 8 or beam positions that are offset from their nominal locations 9 . If the interferometer beam is also offset from its nominal location at the center of the ETM surface, the rotations will be sensed by the interferometer as length variations. Conversion of the unintended rotation to length variation depends on the torque-to-rotation transfer function for the suspended optic, R(ω), in units of 1/(N m), and the dot product between the displacement vector for the center of force 8 LIGO uses a two-beam configuration to minimize local elastic deformation of the ETM surface in the region sensed by the interferometer beam [21, 22] 9 The nominal locations for the two Pcal beams are diametrically opposed and displaced above and below the ETM center by 111.6 mm. These locations are close to the nodal circle of the drumhead natural deformation mode of the ETM and are chosen to minimize sensing of bulk elastic deformation of the ETM by the interferometer [23] . Modeled torque-to-rotation transfer functions for a suspended LIGO end test mass for both pitch and yaw, and for a free mass, -1/(Iω 2 ). Upper-right and Lower-right panels: Discrepancy between the modeled and free-mass transfer functions. Above 20 Hz, the S(ω) discrepancy is less than 0.1 %, and R(ω) discrepancy is less than 0.3 %.
of the Pcal beams, a, and the displacement vector for the interferometer beam, b, as given in (2) . R(ω) is plotted in the lower-left panel of figure 2 for a LIGO ETM for both pitch and yaw rotations. As shown in the lower-right panel of figure 2, at frequencies above 20 Hz, like S(ω), R(ω) is well approximated (within 0.3 %) by the response of a free-mass, R(ω) −1/(Iω 2 ), where I is the moment of inertia of the suspended optic about the relevant axis of rotation.
Taking into account the longitudinal displacement as well as the apparent displacements caused by unintended rotations [21] , the fiducial length modulation induced by Pcal forces, x(ω), is given by
In practice, we do not know the magnitude or direction of the Pcal center of force displacement vector, a; we can only estimate the maximum magnitude, as described in section 3.2. Therefore the second term in the square brackets determines the relative uncertainty in the displacement amplitude introduced by the unintended rotation of the ETM , which we define to be rot . As (2) shows, calculation of the displacements induced by the Pcal systems requires estimates of three parameters: the angle of incidence of the Pcal beams on the ETM surface, the mass of the ETM, and the laser power reflecting from the highly-reflective ETM surface inside the vacuum envelope. Accurate measurement of these parameters, especially P (ω), in order to accurately estimate the Pcal displacements is the focus of this paper. It is organized as follows: in section 2 we discuss methodology for calibrating Pcal power sensors that are located outside the vacuum envelope in terms of the power reflecting from the ETM. Improvements in measurement methods and compensation for temperature variations are also presented. We also discuss a new method for comparing displacement factors from both end stations to calculate combined displacement factors with reduced uncertainty. In section 3 we use measurements made with the LIGO Hanford Observatory (LHO) interferometer during the O3 observing run, from April 2019 to March 2020, to demonstrate the application of the methods described in section 2. We also give detailed descriptions of how uncertainties in the various measured and calculated parameters are estimated. In section 4 we summarize the results presented in this paper and discuss prospects for further improvements in the accuracy and uncertainty of fiducial displacements generated to calibrate gravitational wave detectors. A list of the symbols used in this paper and the parameters they represent is included in table 8 in Appendix A.
The LIGO Pcal systems incorporate a laser power sensor that is located outside the vacuum envelope and receives almost all of the laser power exiting the vacuum chamber after reflecting from the ETM (the Rx sensor in figure 1 ). A second sensor (Tx) samples a small fraction of the laser power directed into the vacuum chamber and is used for power calibration and optical efficiency measurements. Calibrating the Rx sensor in terms of power reflected from the ETM requires compensating for optical losses between the ETM and the sensor.
The power reflectivity of the ETM is greater than 0.9999 but the anti-reflection coated vacuum windows and the relay mirrors located inside the vacuum envelope [17] reduce the optical efficiency, η, between the transmitter and receiver modules to approximately 0.985 -0.990, i.e. the overall optical loss is about 1.0 to 1.5 %. Measurements made inside the vacuum envelope when the system is vented to atmosphere enable apportioning the overall optical efficiency between the input side (between the Tx sensor and the ETM), η T , and the output side (between the ETM and the Rx sensor), η R . The measured efficiency ratio, β = η T /η R , together with the overall optical efficiency that can be measured outside the vacuum envelope, enables calculation of the optical efficiency factors as η R = η/β and η T = η β.
The power reflecting from the ETM is thus estimated by
where P R (ω) is the power measured at the Rx sensor. The power at the ETM can also be calculated as P (ω) = η T P T (ω) with P T (ω) the power measured on the Tx side. However, this sensor is less reliable because it is subject to variations in the beamsplitter that reflects the small sample of the input light and because it is insensitive to changes in the optical efficiency, η.
Calibration of the power sensors at the interferometer end stations is realized by a three-step process as shown schematically in figure 1: i) a transfer standard referred to as the Gold Standard (GS) is calibrated to SI units at NIST in Boulder, Colorado; ii) the calibration of the GS is propagated to additional transfer standards referred to as Working Standards (WS), one for each observatory, by making responsivity ratio measurements in a dedicated laboratory setup at the LHO, and iii) the calibration of a WS is transferred to end station power sensors at each observatory by making responsivity ratio measurements at the end stations.
The GS is sent to NIST annually for calibration. To transfer the GS calibration to the various working standards, a series of responsivity ratio measurements are made in a laboratory at the LHO [24] . This process has been improved by using a spare Pcal transmitter module [17] that incorporates laser power stabilization and delivers two output beams with powers balanced to within 1 %. The GS and one WS are mounted on automated pneumatic slides that alternate the positions of the two detectors between the two output beams as shown schematically in the upper-left panel of figure 3 . Division of the output voltages recorded simultaneously in a given configuration minimizes variations induced by laser power changes; sequential measurements with the detector positions swapped minimizes the impact of changes in the reflectivity of the beamsplitter that separates the two beams.
The transfer standards and the Rx power sensors are comprised of integrating spheres with Spectralon ® interior shells (Labsphere model 3P-LPM-040-SL) and custom-built photodetectors. Although the integrating spheres are largely insensitive to the incident beam position, angle, polarization, and size, they exhibit laser speckle due to the coherence of the laser light that correlates the output time series [25] . This temporal correlation limits the precision of the responsivity ratio measurements and can introduce systematic errors. Figure 4 shows ten, 100 second long time series for both the GS (upper panel) and the WSH (lower panel) standards during responsivity ratio measurements in the LHO laboratory. The temporal correlation of the data in a given time series, together with the lack of correlation with the simultaneously-recorded time series from the other sensor, is the result of laser speckle in the sensor outputs. During responsivity ratio measurements in the LHO laboratory, the impact of the laser speckle is ameliorated by recording shorter time series and swapping sensor positions more frequently.
A typical reported responsivity ratio value, α W G = ρ W /ρ G , where ρ W is the responsivity of the WS and ρ G is that of the GS, is given by the average of 100 measurements. The data for each measurement is comprised of four, five second long time series sampled once per second. The first two are sampled simultaneously with the detectors in the A-B configuration as shown in the upper-left panel in figure 3 , and the last two with the detector positions swapped (B-A). The square root of the product of the ratios of each pair of time series yields an estimate of the responsivity ratio every twenty seconds (five seconds are required to re-position the sensors).
To elucidate this method, individual ratios (A-B and B-A configurations) for 2800 twenty second long measurement suites are plotted in the upper-right panel in figure 3 . Variations in the beamsplitter ratio are evident in the "mirrored" appearance of the two data sets. These variations are minimized when calculating α W 1W 2 , the square root of the product of the ratios from sequential measurements (one point for each data set), as shown in the lower-right panel of figure 3 . The standard deviation of the relative variation of the ratio measurements is 1.0 × 10 −4 as shown in the normalized histogram in the lower-left panel in figure 3 . Note that the Pcal power sensor responsivities have non-negligible dependence on temperature, as discussed in more detail below and in section 3.1. The data in the lower panels of figure 3 were "de-trended" by correcting for the slight variation (1.1 × 10 −4 /K) in relative responsivity due to the changing laboratory ambient temperature (1.3 K) during this fifteen hour long measurement interval.
To propagate the GS calibration to the Pcal power sensors at the end stations, a series of measurements are made with a working standard [26] . They involve placing the WS alternately in the path of one or the other Pcal beam in both the transmitter and receiver modules and recording time series of the WS, Rx, and Tx power sensors using the digital data acquisition system [27] . The end station data acquisition system is used for the WS signal, rather than the digital volt meter that is used for the responsivity ratio measurements in the laboratory setup. This is both for convenience and to synchronize the WS data with the Rx and Tx sensor data. Propagating the WS calibration to the end station thus involves the additional step of measuring the conversion factor, ζ W , between volts registered by the digital volt meter and digital counts reported by the end station data system. This is accomplished via a calibrated voltage source. Synchronizing these time series reduces the impact of laser power variations and yields the Rx/WS responsivity ratio, ρ R /ρ W = α RW . These measurements also yield estimates of the overall optical efficiency, η.
Combining the measurements described above, the responsivity or calibration factor for the Rx end station power sensors is given by
in units of ct/W. Measurements with transfer standards are made in laboratories with differing temperatures. Thus the temperature coefficients of the transfer standards must be taken into account when transferring laser power calibration from NIST to the sensors at the interferometer end stations. For a given power sensor, the temperature dependence of the responsivity can be described by
where ρ 0 is the responsivity of the sensor measured at a reference temperature, T 0 , and κ is the temperature coefficient of the relative responsivity. By including the temperatures at which each measurement in the calibration transfer process is made, (5) can be rewritten as
where T N is the NIST laboratory temperature, T L is the LHO laboratory temperature, and T E is the end station temperature. The factor ξ LN corrects for differences in the GS responsivity measured at the NIST and LHO laboratory temperatures and the factor ξ EL corrects for differences in the WS responsivity measured at the LHO laboratory and end station temperatures. These temperature-related correction factors can be written as
Here κ G and κ W are the temperature coefficients for the GS and WS sensors normalized to their respective responsivities at their reference temperatures, T N and T L .
The temperature-compensated Rx sensor calibration in (7), together with the optical efficiency factor, η R , given by (3) and the digitized output of the Rx sensor, d R , yield the estimated power reflecting from the ETM,
Combining (2), and (9), the digitized output of the Rx sensor can be calibrated in induced ETM displacement (units of m/ct) as
Here we define the displacement factor, X, as
Like LIGO, most GW observatories have implemented, or plan to implement, Pcals at both end stations. The Pcal systems at each end station are calibrated using the same procedures. Laser interferometers are designed and tuned to sense differential arm length variations induced by ETM motion without regard (except for the sign of the relative displacement) to which ETM is moving. Thus, comparing the Pcal fiducials produced at both end stations in the interferometer output signal directly measures the ratio of the Pcal calibrations at each end. This comparison can be used to reduce the uncertainty in the induced displacements due to factors that are not common to both ends.
The X/Y calibration comparison is realized by modulating the two Pcal systems at closely-separated frequencies within the sensitive band of the interferometer. Comparison of the amplitudes of the peaks in fast Fourier transforms (FFTs) of the interferometer output signal with the amplitudes of the associated peaks in the calibrated Pcal end station sensor outputs yields the X/Y Pcal calibration comparison factor, χ XY . If there were no uncertainties in the displacement calibration factors for both end stations, X X and X Y , this factor would be 1. However, errors induced by uncertainties in factors that are not common to both end stations (cos θ, M, η R , ξ EL , α RW , ζ W , and rot ) cause it to deviate from 1.
Using χ XY and these uncertainty estimates, correction factors for each end station can be calculated. The combined displacement factors are defined by
where C X and C Y are the X-arm and Y-arm correction factors and
The superscript c denotes that these displacement factors are calculated from the combination of the X-end and Y-end calibration results. For the special case where the uncertainties are the same at both end stations, the factors are given by
In general, the correction factors are calculated using the weighted geometric mean, µ g , of 1 and χ XY . The weighting factors are given by the inverse of the estimated variances in the end station displacement factors, X X and X Y , due to uncertainty contributions that are not common to both end stations. With the X X weighting factor applied to 1 and the X Y weighting factor applied to χ XY , the correction factors are given by
The uncertainties in these factors are given by the weighted relative standard error on the geometric mean.
As an example of how the methods described in section 2 can be applied to calibrate Pcal displacement fiducials, the results of measurements made before and during the O3 observing run are presented here. Measurements were made for both the LHO and LLO interferometers, but for simplicity only the LHO results are presented. The measured and calculated values of the parameters that contribute to the displacement factors described in section 2, and their associated relative uncertainties, are summarized in the tables below. We follow the convention used by NIST and detailed in [28] , employing Type A, Type B, or Type C evaluations for estimates of relative standard uncertainties for each parameter 10 . For Type A estimates, the number of measurements is noted in parentheses after the relative uncertainties listed in the tables.
Calibration of the on-line power sensors at the end stations is realized by measurement of the six multiplicative parameters on the right-hand side of (7). The GS and WS power sensors were upgraded in 2018 [29] , well before the start of the O3 observing run in April 2019. The modifications included changes in the mounting of the photodetector housing to the integrating sphere port in order to improve mechanical stability and reduce the amplitude of variations due to laser speckle. The transimpedance amplifier was also changed to reduce complexity and minimize dark offsets. The upgraded GS was then calibrated by NIST in December 2018. The reported GS responsivity, ρ G , is −8.0985 V/W. The relative uncertainty is 0.315 % (1 σ) [30] . This GS calibration, reported in volts 11 per watt of laser radiation at 1047 nm, is the foundation upon which the calibration of the LIGO, Virgo, and KAGRA detectors were based during the O3 observing run.
As shown schematically in figure 1 , responsivity ratio measurements performed at the LHO [24] are used to transfer the GS calibration to the array of working standards, WSH, WSL, WSK, WSV (and eventually WSI). Also, responsivity ratio measurements between one or more of the working standards and the GS, both before shipping the GS to NIST and immediately after its return from calibration at NIST, are used to determine whether changes occurred during shipping. Figure 5 shows the WS to GS responsivity ratio, α W G , for WSH over the period from November 2018 to February 2020. The mean of the four measurements made immediately before sending the GS to NIST for calibration in mid-November 2018 differs from the mean of the three measurements made immediately after it was returned in mid-December 2018 by 0.06 %. This indicates 10 Evaluation of uncertainty by statistical analysis of a series of measurements is referred to as Type A. Evaluation of uncertainty using means other than statistical analysis of measurements, based on scientific judgement using all relevant information available, is referred to as Type B. Uncertainty estimates resulting from combinations of Type A and Type B evaluations of uncertainty are referred to as a Type C. 11 Volts reported by a dedicated GS voltmeter, Keithley model 2100. that changes during shipping were minimal. However, as seen in the upper panel of figure 5 , a significant change in α W G (0.36 %) occurred in June 2019 when the responsivity ratio measurement setup was moved to a different building at the LHO. Measurements of the responsivity ratio between WSH and the Rx sensors at the LHO end stations (see figure 6 ) indicate that the change in α W G was caused by a change in the GS responsivity, ρ G , alone. The exact cause of this increase in responsivity has not yet been identified. The lower panel in figure 5 shows the same data as in the upper panel, but with the data plotted with the red points shifted up (multiplied by 1.00359).
The error bars in the upper panel have been magnified by a factor of twenty for better visibility. But, they have not been used for weighting because the variations in the data, each point the mean of a set of between 25 and 1200 measured values, are dominated by systematic variations, not the statistical variations of the measurements within each set. Weighting by the inverse of the square of the standard errors on the means would thus bias the estimate of the overall mean in favor of suites with larger numbers of measurements. The mean and the relative uncertainty of α W G for WSH, from the lower panel of figure 5 , in which the error bars have been omitted, are listed in table 1.
To calibrate the on-line Pcal Rx and Tx power sensors at the end stations, a series of responsivity ratio measurements using a WS are performed [26] . The Rx to WSH responsivity ratios, α RW , for the LHO X-end and Y-end sensors measured between December 2018 and March 2020, are shown in figure 6 . The error bars are estimated from the 240 second long time series recorded for each element of the measurement suite, using the formalism described in Appendix A of [23] . This method uses the standard deviations of the data sets rather than estimating standard errors on the mean values, because the data is correlated due to laser speckle. The measured values of α RW and their relative uncertainties, for both end stations, are listed in table 1.
For measurements of α W G , a digital voltmeter (Keithley Model 2100) is used to record the output of the WS. An analog-to-digital (ADC) converter that is part of the LIGO data acquisition system is used for measurements of α RW at the end stations. The factor that converts the volts measured in the laboratory to the digital counts measured at the end station, ζ W in (5), is measured using a calibrated voltage source (Martel Model IVC-222HPII). Before going to the end station, the ratio of the Martel output to that of the Keithley voltmeter is measured. The Martel voltage source is then used to measure the volt-to-count conversion factor for the ADC of the WS data acquisition system channel. The combination of these two measurements yields ζ W in units of ct/V. These factors are close to, but vary slightly from, the ideal value of 1638.4 ct/V. The measured values of ζ W and their relative uncertainties are listed in table 1.
The remaining two factors on the right-hand side of (7) are ξ LN and ξ EL . These correction factors account for the temperature dependence of the responsivities of the GS and WS power sensors and the differences between the NIST laboratory, the LHO laboratory, and the end station ambient temperatures. To measure the temperature dependence of the WS responsivity, a temperature sensor (Analog Devices, AD590) was bonded to the photodetector transimpedance amplifier circuit board. The WS was then placed in an oven and heated to about 7 K above the ambient laboratory temperature. It was then moved to the responsivity ratio measurement setup and measurements of the WS/GS responsivity ratio were made as the WS cooled to room temperature.
The upper-left panel of figure 7 shows α W G , normalized to the mean value of 1.1172 (see figure 5 ), and the difference between the WS and GS temperatures plotted versus time as the WS cools. The lower-left panel shows a linear, least-squares fit to the data indicating a relative responsivity temperature coefficient for WS, κ W , of 4.38 × 10 −4 /K. The uncertainty in the fit is 0.06 × 10 −4 /K.
To investigate the temperature coefficient of the GS, we recorded data with both the WS and the GS at the laboratory ambient temperature during a period when the ambient temperature varied by about 1.3 K. The laboratory ambient temperature and α W G normalized to 1.1172 are plotted in the upper-right panel in figure 7 . A linear, least-squares fit to the data, yielding a slope of −2.35 × 10 −4 /K with uncertainty of 0.07 × 10 −4 /K, is shown in the lower-right panel. Using κ W , from the data in the left panels of figure 7 , the inferred coefficient for the GS, κ G , is 6.73 × 10 −4 /K.
Temperature differences between the NIST, the LHO, and end station measurement environments were quantified using a set of digital thermometers that were huddled to assess relative offsets, then deployed to each location. The NIST measurement laboratory temperature was set to 20 • C, the mean temperature of the LHO responsivity ratio measurement laboratory was 23 • C, and the mean X-end and Y-end temperatures were 19.8 • C and 19.7 • C respectively. The end station temperatures varied by about ± 0.5 • C over six months. The calculated values for ∆T LN and ∆T EL , together with the measured GS and WS temperature coefficients, κ G and κ W , are used to calculate ξ LN and ξ EL using (8) . Their values and associated relative uncertainties are listed in table 2. Figure 7 .
Temperature-dependent responsivity ratio measurements enabling estimation of the WS and GS temperature coefficients, κ W and κ G . In all panels the responsivity ratio, α W G , is normalized to its mean value of 1.1172 (see figure 5 ). Upper-left panel: Normalized responsivity ratio and temperature difference between the WS and GS versus time as the WS cools after being heated in an oven. Lowerleft panel: Linear least-squares fit to the data in the upper-left panel. The slope and uncertainty are listed in the legend. Upper-right panel: Normalized responsivity ratio and the ambient laboratory temperature plotted versus time when both sensors were at the ambient temperature that varied more than 1 K over almost ten hours. Lowerright panel: Linear least-squares fit to the data in the upper-right panel with slope and uncertainty listed in the legend.
As noted in section 2.1, the end station measurements made with the WS also yield measurements of the overall optical efficiency, η, for propagation of the laser beams between the Tx and Rx modules. Assuming that β is constant, using (3) the WS end station measurements yield estimates of the optical efficiency, η R , required to calibrate the end station power sensors in terms of laser power reflecting from the ETM using (4). The measured values of η for the LHO X-end and Y-end stations from November 2018 to February 2020 are plotted in figure 8 . The error bars, used for calculating weighted values, are generated using the formalism detailed in Appendix B of [23] . The mean values of η R , β, and η, together with their relative uncertainties, are listed in table 3.
Using right hand side of (2), in which suspension displacement and rotation transfer functions have been approximated by the responses of a free mass, only the mass of the ETM is required to convert the force given by (1) to the displacement factor, X, in (10) . Here, we are neglecting the second term in square brackets because its estimated magnitude is treated as an uncertainty, rot . The ETM masses are measured by the vendor that polishes the mirrors and at the end stations before they are suspended. The calibration of the electronic balances used before suspending have been verified using two 20 kg calibrated reference masses. The maximum estimated uncertainty in the ∼ 40 kg suspended mass is 10 g. The angle of incidence, θ, is determined by the last relay mirrors on the periscope structure located inside the vacuum envelope [17] . These mirrors direct the Pcal beams such that they impinge on the ETM at positions above and below center, on the vertical center line of the face of the ETM. The nominal angle of incidence is 8.72 degrees. The maximum deviations of this angle are bounded by the size of the periscope optics (2 inch diameter) that relay the beams to the end test mass.
When the interferometer beam is not centered on the face of the ETM, as is currently the practice at the LIGO observatories to mitigate the impact of point absorbers in the mirror coatings [31] , Pcal-induced rotations can increase or decrease the sensed ETM displacement. Force-to-angle coupling investigations using the ETM orientation actuators have been used to infer the interferometer beam offsets, but for Pcal beam position offsets we currently only have estimates of their maximum magnitudes.
The Pcal beams are positioned on the ETM surface when the vacuum envelope is vented, using targets that are bolted to the suspension structure surrounding the ETM (see left image in figure 9 ). The reflected beams exiting the vacuum envelope are aligned to the center of the entrance aperture of the Rx power sensor (see right image in figure 9 ). Monitoring of the locations of the beams within the aperture, ∼ 12 m from the Tx module output, indicates that the maximum displacement of the beams from their optimal locations on the ETM surface is ± 2 mm. The powers in the two beams are balanced to within 1 %. With their nominal locations 111.6 mm from the center of the ETM surface, this power imbalance causes the center of force to be displaced by an additional ∼ 0.5 mm. The magnitude of a is estimated by adding the Pcal beam position offset and power imbalance contributions in quadrature.
The interferometer beam position offsets from center are 29 mm for the X-end ETM and 22 mm for the Y-end ETM. The uncertainty introduced by unintended rotation of the ETM, rot , is proportional to the dot product of the Pcal and interferometer beam offset vectors, a and b, i.e. rot ∝ | a|| b| cos φ. Because φ, the angle between a and b, is equally probable to be any value between -π and π, we use a sine wave, or U-shaped, probability density function [32] to estimate the variance in cos φ and form a Type B estimate for rot using (2), rot = M ab/( √ 2I). The values of the relative uncertainty estimates for both the X-end and Y-end stations at LHO are listed in table 4. Because of these large interferometer beam position offsets, rot is currently one of the largest sources of uncertainty for the LIGO Pcal systems. Combining the factors discussed above, the displacement factors for each end station are calculated using (11) . The values of X for both end stations, together with their relative uncertainties are listed in table 5. The overall uncertainties are dominated by the uncertainties due to the unintended rotation of the ETM, rot , and the uncertainty in the calibration of the end station power sensors, ρ R . The relative uncertainties of 0.53 % for the X-end and 0.45 % for the Y-end are smaller than the lowest values previously reported, 0.75 % [17] .
To compare the Pcal calibrations at the two end stations, large-amplitude displacements at frequencies separated by 0.1 Hz were induced by the X-and Y-arm Pcals as shown in the right panel in figure 10 . The calibrated Rx power sensors at each end station provide estimates of the induced displacement amplitudes (see (10) ). The induced periodic displacements appear with high SNR in the interferometer output signal and the Rx power sensor signals as shown in the right panel of figure 10 . The Spectral Line Monitoring tool [33] was used to perform long-duration (100 s) FFT measurements of the line amplitudes in the interferometer and Rx sensor output time series. The left panel of figure 10 shows χ XY , the ratio (X/Y) of the amplitudes of the calibrated X-end and Y-end Pcal Rx sensor output signals, each normalized to the amplitude of the respective peak in the interferometer output signal. To ensure that small variations in the interferometer response function over the 0.1 Hz frequency separation between the two excitations were not impacting the comparison, the excitation and analysis were repeated with the frequencies of the X-end and Y-end Table 6 . Calculated X-end and Y-end combined displacement correction factors, C X and C Y , together with χ XY , µ g , and the non-common factors contributing to end station displacement factor uncertainty (indented), and their uncertainties, for the LHO interferometer during the O3 observing run. For Type A uncertainties, the number of measurements is noted in parentheses.
Param figure 10 ) is 1.00460 with a relative standard uncertainty of 0.014 % and for the other data set (red points in the left panel of figure 10 ) the mean value is 1.00463 with a relative standard uncertainty of 0.046 %. The contributions to the relative uncertainties in X from factors that are not common to both end stations (see section 2.2) are 0.42 % for X-end and 0.32 % for Y-end. Summing these uncertainties in quadrature, the calculated relative uncertainty for the quotient of the displacement factors, X X /X Y , is 0.52 %. Thus the measured value of 1.0046 for χ XY is a ∼ 0.9 σ result.
To calculate combined end station displacement factors, we follow the formalism detailed in section 2.2. We calculate the weighted geometric mean of 1 and 1.0046, and its relative uncertainty, using weighting based on the 0.42 % and 0.32 % non-common uncertainty contributions for X X and X Y . The combined displacement correction factors, C X and C Y , calculated using (14) together with χ XY , µ g and the non-common factors contributing to end station displacement factors, are listed in table 6 with their relative uncertainty estimates. The combined displacement factors, X c X and X c Y , calculated using (12) and their relative uncertainties are listed in table 7. The uncertainties for these combined displacement factors have been estimated by summing the correction factor uncertainties, 0.25 %, in quadrature with the uncertainty of 0.32 % from factors that are common to both end stations, ρ G , α W G and ξ LN . The resulting overall uncertainties of 0.41 % for the combined displacement factors are smaller than those for the displacement factors for each end station, 0.53 % for X X and 0.45 % for X Y . This reduction results from combining the calibrations from both end stations using the measured X/Y comparison factor, χ XY . This lowers the uncertainty contributions from sources that are not common to both end stations from 0.42 % for X-end and 0.32 % for Y-end to 0.25 % for C X and C Y .
The estimated overall 1 σ relative standard uncertainties of 0.41 % for both the X c X and X c Y displacement factors are significantly lower than previously reported values [23] . The reduced uncertainty in the calibration of the GS at NIST (from 0.42 % to 0.32 %) is a key contributor to this improvement. The updated factors are corrected for the temperature dependence of the GS and WS sensors and the ambient temperature differences between laboratories where the measurements were performed (ξ LN and ξ EL ). They also incorporate in-chamber measurements of the ratio between the input-side and output-side optical efficiencies, β. These enable apportioning the overall optical efficiency (η) measured outside the vacuum envelope with the WS to calculate the optical efficiency factors, η R and η T .
Unintended rotation of the ETM due to Pcal forces is one of the biggest contributors to the estimated overall calibration uncertainty. The large displacements of the interferometer beams from the centers of the ETMs (29 mm and 22 mm) are the main cause. The interferometers have been designed to operate with the beams located close to the center of the ETM surfaces. However, point defects in the ETM high-reflectivity coatings [31] have required large beam displacements to optimize interferometer performance, even while suffering the deleterious impacts of operating with mis-centered beams. If the coatings are improved to eliminate these defects and the interferometer beams can be moved to within a few millimeters of center, the uncertainty due to unintended rotation could be reduced by a factor of ten. Finally, using the interferometer to compare the X-end and Y-end Pcal displacement factors (by measuring χ XY ) has enabled generating combined displacement factors, X c X and X c Y , with improved accuracy and lower estimated uncertainty. Realizing the important role laser power sensor calibration plays in the scientific impact of gravitational wave detections, there is increased interest and activity within the national metrology institute (NMI) community to improve power sensor calibration accuracy and precision. A comparative study executed in 2009 under the auspices of the Consultative Committee for Photometry and Radiometry (CCPR) of the Bureau International des Poids et Mesures (BIPM) [34] indicated that there might be systematic differences that are as large as a few percent between power sensor calibrations performed at various NMIs. In March 2019 the GW Metrology Workshop was held at NIST in Boulder, CO to communicate the requirements of the GW community and investigate possibilities within the NMI community for addressing them. A bi-lateral comparison between NIST and the Physikalisch-Technische Bundesanstalt (PTB) in Braunschweig, Germany using a LIGO working standard is currently underway. Additional bi-lateral comparisons are being considered. Furthermore, scientists at NIST are developing a new generation of primary calibration standards using bolometers and nanotube absorbers [35] that are expected to have significantly lower uncertainties than the standards currently being used. These sensors may prove suitable for locating inside the vacuum envelope where several sources of uncertainty would be mitigated, providing laser power calibration directly traceable to SI units in-situ [36] .
There is another significant ongoing effort within the GW community to establish a displacement reference for the Pcals that is similar to what is provided by the NIST calibrations of the GS. It involves deploying rotating mass quadrupoles, Newtonian (Ncal) or Gravity (Gcal) calibrators, near, but outside the ETM vacuum chambers at the observatory end stations [19, 37] . Modulation of the position of the ETM is achieved via variation of the Newtonian gravitational field produced by the rotating quadrupole. The possibility of using similar devices has been considered for a long time [38] . Prototypes are currently installed and being commissioned at the both the Virgo and the LHO observatories.
Establishing continuous, on-line displacement fiducials with sub-percent accuracy is a key element of interferometer calibration. But achieving sub-percent calibration of the interferometer output signals, over the entire sensitive frequency band, and over long observing runs during which changes to improve interferometer performance may have been implemented, is a daunting hurdle. To better characterize the interferometer response function, swept-sine measurements are performed regularly using the Pcal systems. These, together with tracking of time-varying interferometer parameters using continuously-injected Pcal modulations [39] , have contributed to dramatically improved overall interferometer calibration uncertainties for the LIGO detectors. In some frequency bands, the estimated systematic errors during the recently-completed O3 observing run are below 2 % [40] .
Currently, interferometer calibration uncertainty does not limit the extraction of source parameters from detected GW signals. However, as the SNRs of GW signals increase, optimally extracting the encoded information they carry will require reducing interferometer calibration uncertainty. The improvements in the accuracy of the Pcal fiducial displacements described herein should enable realizing the sub-percent interferometer calibration accuracy required to continue extending the scientific reach of the global network of GW detectors. 