## 3 Two-Way Ranging Formulations
\begin{split}\varphi_{1}(t_{1}+\delta t_{1})&= \varphi_{1}^{\text{TX}}(t_{1}+\delta t_{1}-\tau_{21}-\tau_{12})+\varphi_{2}^{ \text{OPLL}}(t_{1}+\delta t_{1}-\tau_{21})\\ &-\varphi_{1}^{\text{TX}}(t_{1}+\delta t_{1})-N_{1}-\varepsilon _{1},\\ \varphi_{2}(t_{2}+\delta t_{2})&=\varphi_{2}^{\text{ OPLL}}(t_{1}+\delta t_{2}-\tau_{21})+N_{2}+\varepsilon_{2}.\end{split}  \tag{6}\]  
The two-way ranging phase of LRI at nominal time \(t_{1}\) is computed as the difference of the transponder and master phase:  
\[\varphi_{\text{LRI}}(t_{1})=\varphi_{2}(t_{2}+\delta t_{2})-\varphi_{1}(t_{1} +\delta t_{1}). \tag{7}\]  
By substituting Equation (6) into Equation (7) and applying the Taylor expansion since the light travel time and time tag error are small, we have  
\[\varphi_{\text{LRI}}(t_{1})\approx f_{1}(\tau_{12}+\tau_{21})+f^{\text{OPLL} }(\delta t_{2}-\delta t_{1})+(N_{1}+N_{2})+(\varepsilon_{1}+\varepsilon_{2}), \tag{8}\]  
Figure 2: Phase observables of the two-way ranging system at actual time and their combinations. Detailed descriptions are in the main text. where \(f_{1}\) and \(f^{\text{OPLL}}\) represent the master laser frequency and the heterodyne frequency in the OPLL, respectively. Each frequency can be divided into a constant part \(\bar{f}\) and a tiny variable part \(\delta f\). Therefore, one can write  
\[f_{1} =\bar{f}_{1}+\delta f_{1}, \tag{9}\] \[f^{\text{OPLL}} =f^{\text{OPLL}}+\delta f^{\text{OPLL}}.\]  
The two-way ranging phase is then written as  
\[ \begin{split}\varphi_{\text{LRI}}(t_{1})&\approx \bar{f}_{1}(\tau_{12}+\tau_{21})+\delta f_{1}(\tau_{12}+\tau_{21})+f^{\text{ OPLL}}(\delta t_{2}-\delta t_{1})+\delta f^{\text{OPLL}}(\delta t_{2}- \delta t_{1})\\ &+(N_{1}+N_{2})+(\varepsilon_{1}+\varepsilon_{2}),\end{split}  \tag{10}\]  
the unit of which can be converted from cycle to meter by using the nominal laser frequency \(f_{0}\). Therefore, the one-way range \(\rho_{\text{LRI}}\) with the unit of meter reads  
\[ \begin{split}\rho_{\text{LRI}}(t_{1})&=\frac{c}{2f _{0}}\varphi_{\text{LRI}}(t_{1})\\ &\approx c\frac{\bar{f}_{1}}{2f_{0}}(\tau_{12}+\tau_{21})+c\frac{ \delta f_{1}}{2f_{0}}(\tau_{12}+\tau_{21})+c\frac{f^{\text{OPLL}}}{2f_{0}}( \delta t_{2}-\delta t_{1})\\ &+c\frac{\delta f^{\text{OPLL}}}{2f_{0}}(\delta t_{2}-\delta t_{ 1})+\frac{c}{2f_{0}}(N_{1}+N_{2})+\frac{c}{2f_{0}}(\varepsilon_{1}+\varepsilon _{2}).\end{split}  \tag{11}\]  
For gravity field recovery, the measured range needs to be converted into the instantaneous range between satellites, which means the distance between two satellite positions at exactly the same time. Here, we use \(\rho_{\text{inst}}\) to represent the instantaneous range and we define an equivalent light time \(\tau\) corresponding to the instantaneous range as \(\tau=\rho_{\text{inst}}/c\), where \(c\) denotes the speed of light in a vacuum. The first two terms in Equation (11) can then be written as  
\[c\frac{\bar{f}_{1}}{2f_{0}}(\tau_{12}+\tau_{21})+c\frac{\delta f_{1}}{2f_{0}} (\tau_{12}+\tau_{21})=\frac{\bar{f}_{1}}{f_{0}}\rho_{\text{inst}}+\frac{\bar {f}_{1}}{\bar{f}_{0}}\delta\rho_{\text{path}}+\frac{\delta f_{1}}{f_{0}}\rho _{\text{inst}}+\frac{\delta f_{1}}{f_{0}}\delta\rho_{\text{path}}, \tag{12}\]  
where the pathlength error \(\delta\rho_{\text{path}}\) due to the non-zero light travel time is defined as  
\[\delta\rho_{\text{path}}=c\left(\frac{(\tau_{12}+\tau_{21})}{2}-\tau\right). \tag{13}\]  
By dropping the small terms in Equations (11) and (12), the intersatellite range measurement can be written as  
\[ \begin{split}\rho_{\text{LRI}}(t_{1})&\approx\frac{ \bar{f}_{1}}{f_{0}}\rho_{\text{inst}}+\frac{\bar{f}_{1}}{f_{0}}\delta\rho_{ \text{path}}+\frac{\delta f_{1}}{f_{0}}\rho_{\text{inst}}+c\frac{f^{\text{OPLL }}}{2f_{0}}(\delta t_{2}-\delta t_{1})\\ &+\frac{c}{2f_{0}}(N_{1}+N_{2})+\frac{c}{2f_{0}}(\varepsilon_{1}+ \varepsilon_{2}),\end{split}