### Light wave-packets propagating in the Earth’s space-time

In this section, we will give a brief introduction to the propagation of photons under the influence of the Earth’s gravitational field. Considering the rotation of the Earth, the space-time considered in this analysis can be approximately described by the Kerr metric^{26} and our work will be performed on the equatorial plane for simplicity. Now the Kerr line element in the Boyer-Lindquist coordinate ((t,r,phi )) is reduced to^{26}

$$begin{aligned} ds^2&= -,Big (1-frac{2M}{r} Big )dt^2+frac{1}{Delta }dr^2 nonumber \&quad +Big (r^2+a^2+frac{2Ma^2}{r}Big ) dphi ^2 – frac{4Ma}{r} dt , dphi , end{aligned}$$

(1)

$$begin{aligned} Delta&=1-frac{2M}{r}+frac{a^2}{r^2}, end{aligned}$$

(2)

where *M* and *r* respectively denote the mass and radius of the rotating planet. The Kerr parameter (i.e., normalized angular momentum) can be expressed as (a=frac{J}{M}), a combination of the planet’s angular momentum (*J*) and mass *M*.

In order to clearly describe the propagation of wave-packets from a source on Earth to a receiver satellite at a fixed distance, two observers called Alice and Bob are respectively prepared as the reference frames. More specifically, focusing on a photon sent from Alice at the time of (tau _A), it will arrive at Bob at the time of (tau _B=Delta tau +sqrt{f(r_B)/f(r_A)}tau _A), where *f*(*r*) is the gravitational frequency shifting factor and the (Delta tau) represents the propagation time of the light from Alice to Bob. As is well known, photons can be modeled by the wave packet of massless bosonic field with a distribution of (F^{(K)}_{Omega _{K,0}}), where (Omega _{K}) is the mode frequency peaked at (Omega _{K,0})^{27,28} and (K=A, B) denotes the mode in Alice’s or Bob’s reference frame. For an observer infinitely far away from Alice or Bob, the annihilation operator takes the form of

$$begin{aligned} hat{a}_{Omega _{K,0}}(t_K)=int _0^{+infty }dOmega _K e^{-iOmega _K t_K}F^{(K)}_{Omega _{K,0}}(Omega _K)hat{a}_{Omega _K}, end{aligned}$$

(3)

with the frequency distribution of (F^{(K)}(Omega )). Such operator is naturally applied to modeling a wave packet of the electromagnetic field located and propagating in the space-time. Note that when the frequency distribution (F^{(K)}(Omega )) is normalized (i.e., (int _{Omega >0}|F^{(K)}(Omega )|^2=1)), the creation (hat{a}^{dagger }_{Omega _{K,0}}) and annihilation (hat{a}_{Omega _{K,0}}) operators will satisfy the canonical equal time bosonic commutation relations, (([hat{a}_{Omega _{K,0}}(t),hat{a}^{dagger }_{Omega _{K,0}}(t)]=1)).

Now let us consider a realistic case in which Alice (located on the surface of the Earth, (r_A=r_E)) prepares and sends a wave packet (F^{(A)}_{Omega _{A,0}}) to Bob (located on a satellite at the altitude of (r_B)). Given the effects of the Earth’s gravitational field, the wave packet received by Bob ((F^{(B)}_{Omega _{B,0}})) should be modified, following the relation between the frequency distributions of two wave packets^{9,10}. Meanwhile, the gravity of the Earth also changes the mode frequency (Omega _K) with the shift of (Omega _A=sqrt{f(r_A)/f(r_B)}Omega _B), where (f(r_A)/f(r_B)) represents the shifting function of total gravitational frequency (see more details in the text). Therefore, the total modification induced by the Earth’s gravitational field can be parameterized as

$$begin{aligned} F^{(B)}_{Omega _{B,0}}(Omega _B)=root 4 of {frac{f(r_B)}{f(r_A)}}F^{(A)}_{Omega _{A,0}}left( sqrt{frac{f(r_B)}{f(r_A)}}Omega _Bright) . end{aligned}$$

(4)

One may clearly see that the effect induced by the curved space-time of the Earth cannot be simply corrected by a linear shift of frequency, which indicates the difficulty of compensating such transformation in realistic implementations.

Fortunately, following the relation between such nonlinear gravitational effect and the fidelity of the quantum channel^{9,10}, it is possible to decompose the mode (bar{a}^{prime }) received by Bob into

$$begin{aligned} bar{a}^{prime }=Theta hat{a}^{prime }+sqrt{1-Theta ^2}hat{a}_{bot }^{prime }, end{aligned}$$

(5)

in terms of the mode prepared by Alice ((a^{prime })) and its orthogonal mode (hat{a}_{bot }^{prime }) (i.e., ([hat{a}^{prime },hat{a}_{bot }^{prime dagger }]=0))^{29}. Here (Theta) is the wave packet overlap between the distributions (F^{(B)}_{Omega _{B,0}}(Omega _B)) and (F^{(A)}_{Omega _{A,0}}(Omega _B)), which takes the form of

$$begin{aligned} Theta :=int _0^{+infty }dOmega _B,F^{(B)star }_{Omega _{B,0}}(Omega _B)F^{(A)}_{Omega _{A,0}}(Omega _B). end{aligned}$$

(6)

It is easy to see that (Theta =1) corresponds to a perfect channel, while (Theta <1) represents a noisy channel under the influence of the Earth’s curved space-time. In order to better characterize the frequency distribution of the source, we introduce a specific quantity, i.e., fidelity of ({F}=|Theta |^2) in the following analysis and apply a real normalized Gaussian wave packet to Alice’s mode

$$begin{aligned} F_{Omega _0}(Omega )=frac{1}{root 4 of {2pi sigma ^2}}e^{-frac{(Omega -Omega _0)^2}{4sigma ^2}}, end{aligned}$$

(7)

with the wave packet width of (sigma). We remark here that in the expression of the overlap parameter (Theta) (Eq. (6)), the integration will be performed over strictly positive frequencies, which is justified by the fact that the peak frequency is typically much larger than the spreading of the wave packet (i.e.,(Omega _0gg sigma)). The combination of Eqs. (4) and (7) provides us with

$$begin{aligned} Theta =sqrt{frac{2}{1+(1+delta )^2}}frac{1}{1+delta }e^{-frac{delta ^2Omega _{B,0}^2}{4(1+(1+delta )^2)sigma ^2}}, end{aligned}$$

(8)

where the new parameter (delta) is introduced to quantify the shifting effect

$$begin{aligned} delta =root 4 of {frac{f(r_A)}{f(r_B)}}-1=sqrt{frac{Omega _B}{Omega _A}}-1. end{aligned}$$

(9)

Focusing on the equatorial plane of the Earth described by Kerr metric, the parameter of (frac{Omega _B}{Omega _A}) may be rewritten as^{30}

$$begin{aligned} frac{Omega _B}{Omega _A}=frac{1+epsilon frac{a}{r_B}sqrt{frac{M}{r_B}}}{Csqrt{1-3frac{M}{r_B}+ 2epsilon frac{a}{r_B}sqrt{frac{M}{r_B}}}}. end{aligned}$$

(10)

Here the normalization constant takes the form of (C=[1-frac{2M}{r_A}(1+2a {omega })+big (r^2_A+a^2-frac{2Ma^2}{r_A}big ){omega }^2]^{-frac{1}{2}}) (with the Earth’s equatorial angular velocity (omega)) and (epsilon =pm 1) stands for the direction of the corresponding orbit (i.e., the satellite co-rotates with the Earth when (epsilon =+1)).

Now we will expand Eq. (9) by keeping the first order of the perturbation term ((r_A omega )^2), in order to derive the explicit expression of the frequency shift ((delta)) for a photon exchanged between Alice and Bob. Considering the independency between the perturbative result and the state of the Earth and the satellite (i.e., whether they are co-rotating with each other), one can easily obtain the shift parameter with the following expression

$$begin{aligned} delta&= delta _{Sch}+delta _{rot}+delta _h\&= frac{1}{4}frac{r_S}{r_A}big (frac{r_A-2h}{r_A+h} big )-frac{(r_Aomega )^2}{2}-frac{(r_Aomega )^2}{4}big (frac{3}{4}frac{r_S}{r_A}-frac{2r_Sa}{omega r_A^3}big ), end{aligned}$$

where (delta _{Sch}), (delta _{rot}), and (delta _h) respectively denote the first-order Schwarzschild term, rotation term, and higher-order correction term. Note that (r_S=2M) is the Schwarzschild radius of the Earth, while the parameter (h=r_B-r_A) quantifies the height difference between Bob and Alice. It should be pointed out that when (h=frac{r_A}{2}), the overlap parameter (Theta) is no longer equal to one due to the effects induced by the rotation of the Earth. However, when the satellite moves at the height of (hsimeq frac{r_A}{2}), the combined effects of the Earth’s gravity and Special Relativity (i.e., the Doppler effect generated by the motion of the satellite) will compensate each other ((Theta =1)). Therefore, the photon received by Bob at this height will not experience any frequency shift, which implies that the clock rate of Bob will become equal to that of Alice.

### The influence of Earth’s curved space-time on three types of correlations

In this section, we propose a scheme to test the classical and quantum correlations at long distance, and furthermore quantify the effects of the Earth’s curved space-time on such correlations. In the framework of a pair of entangled photons initially prepared in a two-mode squeezed state (with the modes of (b_1) and (b_2) on the Earth’s surface), we firstly send a photon (with the mode (b_1)) to Alice, with the other photon (with the mode (b_2)) propagating from the Earth to the satellite and then received by Bob. Now the wave packet of photons will be deformed by the curved background space-time of the Earth, based on which the total correlation, classical correlation, and quantum correlation (quantum discord) of this system will be investigated in the curved space-time.

Considering the fact that the two modes ((b_1) and (b_2)) will be received by Alice and Bob at different heights (i.e., within different satellite orbits), the effects caused by the curved space-time should be taken into account in the total process. As was extensively discussed in^{9,10}, such space-time effects on the two-mode squeezed state can be modeled by a beam splitter with the orthogonal modes of (b_{1bot }) and (b_{2bot }). The covariance matrix of the initial two-mode squeezed state is given by

$$begin{aligned} Sigma ^{b_1b_2b_{1bot }b_{2bot }}_0=left( begin{array}{cc} tilde{sigma }(s) & 0 \ 0 & I_4 end{array}right) , end{aligned}$$

(11)

where ({I}_4) represents the (4times 4) identity matrix and (tilde{sigma }{(s)}) denotes the covariance matrix of the two-mode squeezed state

$$begin{aligned} tilde{sigma }(s)=left( begin{array}{cc} cosh {(2s)} {I}_2&sinh {(2s)}sigma _z \ sinh (2s)sigma _z &cosh {(2s)} {I} _2 end{array}right) , end{aligned}$$

(12)

with the Pauli matrix (sigma _z) and squeezing parameter *s*. Following the recent analysis of^{9,10}, one can use lossy channels to model such effects and we only consider Bob’s mode sent to the satellite and Alice’s mode sent to the ground, which means that Alice will experience a perfect channel ((Theta _1=1)). Therefore, in our scheme one may naturally obtain

$$begin{aligned} bar{b}_1 &= ,b_1nonumber \ bar{b}_2& = Theta _2,b_2+sqrt{{1-Theta _2^2}}b_{2bot }, end{aligned}$$

(13)

and such process can be described by a mixed beam splitting of two modes, (b_1(b_2)) and (b_{1bot }(b_{2bot })). For the entire process, the symplectic transformation may be encoded into the Bogoloiubov transformation

$$begin{aligned} S=left( begin{array}{cccc} {I}_2 &0& 0&0 \ 0&Theta _2 {I} _2 &0&sqrt{ 1-Theta _2^2} {I} _2\ 0 &0& -{I}_2&0 \ 0&sqrt{ 1-Theta _2^2} {I} _2 &0&-Theta _2 {I} _2 end{array}right) , end{aligned}$$

based on which the final state after transformation is defined as (Sigma ^{b_1b_2b_{1bot }b_{2bot }}=S,Sigma _0^{b_1b_2b_{1bot }b_{2bot }},S^{T}). Then one could trace over the orthogonal modes ((b_{1bot },b_{2bot })) and obtain the covariance matrix (Sigma ^{b_1b_2}) for the modes ((b_1) and (b_2)) after propagation

$$begin{aligned} Sigma ^{b_1b_2}=left( begin{array}{cc} (1+2sinh ^2s ) {I}_2 &sinh {(2s)},Theta _2,sigma _z \ sinh {(2s)},Theta _2,sigma _z &(1+2sinh ^2s,Theta _2^2 ), {I}_2 end{array}right) . end{aligned}$$

(14)

In our analysis, Eq. (14) will be used to characterize the final form of the two-mode squeezed state that suffers from the influence of the Earth’s gravity.

Let’s remark here that a lossy quantum channel determined by the wave packet overlap parameter (Theta) (which contains the parameters of (delta), (sigma) and (Omega _{B,0})), could fully describe the effect of the Earth’s curved space-time on the quantum state. For a typical case that satellites stay within the geosynchronous satellite orbit, one could straightforwardly obtain the shift parameter as (delta sim 2.5times 10^{-10}), based on the value of the Schwarzschild radius of the Earth ((r_S=9) mm). In addition, we also focus on a typical parametric down converter crystal (PDC) source^{42,43} with a wavelength of 598 nm and Gaussian bandwidth of (sigma =1)MHz (which corresponds to the peak frequency of (Omega _{B,0}= 500) THz). Given the condition of (delta ll (frac{Omega _{B,0}}{sigma })^2ll 1), we can expand the wave packet overlap (Theta) by keeping the second-order term of the shift parameter, i.e., (Theta sim 1-frac{delta ^2Omega _{B,0}^2}{8sigma ^2}). In order to ensure the validity of perturbation expansion for the Matrix element, we furthermore estimate the value of (frac{delta ^2Omega ^2_{B,0}}{8sigma ^2}sim 3.2times 10^{-8}), based on which one may guarantee the validity of the perturbative expansion with the squeezing parameter (sll 7.6) (which corresponds to (sinh ^2(s)ll 10^6)). For convenience, the peak frequency and the Gaussian bandwidth are respectively rescaled to

$$begin{aligned} tilde{Omega }equiv frac{Omega }{Omega _{B,0}}, tilde{sigma }equiv frac{sigma }{sigma _0}, end{aligned}$$

(15)

with (Omega _{B,0}=500) THz and (sigma _0=1) MHz. In the following analysis, the dimensionless parameters (tilde{Omega }) and (tilde{sigma }) will be abbreviated as (Omega) and (sigma), respectively.

Summarizing, the effects of the Earth’s curved space-time will quantified by three types of correlations, i.e., the total correlation ({I}_2), classical correlation ({J}_2) and quantum discord ({D}_2) between the two modes of (b_1) and (b_2) (see the Methods for more details). The behaviors of ({I}_2), ({J}_2) and ({D}_2) are explicitly shown in in Fig. 1, which illustrates the three types of correlations as functions of increasing orbit height *h*. The Gaussian bandwidth, squeezing parameter and frequency of mode (b_2) are respectively fixed at (sigma =1), (s=1) and (Omega _B=1). One may clearly see that, compared with the classical correlation, the quantum discord will be more easily effected by the Earth’s curved space-time. Such tendency has been firstly noted and extensively discussed in the previous works^{44}. Moreover, the three type of correlations between mode (b_1) and (b_2) have exhibited very similar behaviors, i.e., they all initially increase for a specific range of height parameter (hsimeq frac{r_A}{2}) and then gradually approach to a finite value with increasing *h*. One possible explanation to such findings lies in the different roles played by gravitational frequency shift and special relativistic effects, with the increase of the satellite’s height. More specifically, the photon’s frequency received by the satellites at the height of (h<frac{r_A}{2}) will experience blueshift (with increasing correlations), while the corresponding frequencies received at the height of (h>frac{r_A}{2}) will experience redshift (with decreasing correlations). Actually, in the framework of three types of correlations between the photon pairs, the peak values of all correlations (i.e., the parameter (delta =0)) have strongly suggested a detectable frequency transformation from blueshift to redshift^{30}. It should be pointed out that the total frequency shift generated by both special and general relativistic effects should be taken into account (Eq. (9)), when the two parties are located at the same height or in the flat space-time ((delta ne 0))^{30}. In addition, when the satellite is moving at the height of (h=frac{r_A}{2}) with vanishing Schwarzschild term ((delta _{Sch})), photons received on the satellite will generate a tiny frequency shift, in the case of which the lowest-order rotation term (delta _{rot}) and higher-order correction term (delta _{h}) should be taken into consideration.

In order to better understand the relations between the three different correlations and the initial squeezing parameter, in Fig. 2 we also show the evolution of ({I}_2), ({J}_2) and ({D}_2) with *s*, with fixed value for the orbit height at geostationary Earth orbits ((h=3.6times 10^4) km), the frequency of mode ((b_2=1)), and the Gaussian bandwidth ((sigma =1)). It is apparent that although all of the three type of correlations will increase with the squeezing parameter, the mutual information ({I}_2) is much more sensitive to the change of squeezing parameter *s*, compared with the classical correlation and the quantum discord. Similarly, considering the degeneracy between the wave packet overlap parameter (Theta) and the frequency parameter (Omega _B), it is also necessary to investigate the change of ({I}_2), ({J}_2) and ({D}_2) with (Omega _B). The results are illustrated in Fig. 3. Besides the general tendency that all types of correlations will decrease with the frequency parameter, one could also find that the quantum discord ({D}_2) is more sensitive to the change of frequency parameter (Omega _B). Such conclusion strongly implies the future possibility of choosing appropriate parameters to realize quantum communications from the Earth to the satellite.

Finally, with the aim of better quantifying the influence of the Earth’s curved space-time, we will define an additional quantity to describe the change rate of three types of correlations

$$begin{aligned} mu _C=frac{C(Sigma ^{b_1b_2})-C_0(Sigma ^{b_1b_2})}{C_0(Sigma ^{b_1b_2})}, end{aligned}$$

(16)

where (C={I}_2, {J}_2,{D}_2) and the subscript 0 denotes its corresponding value at the height of the satellite ((h=0)). In Fig. 4 we plot the change rate of correlation (mu _C) as a function of the height parameter *h*, with fixed Gaussian bandwidth (sigma =1) and frequency parameter (Omega _B=1). The behavior of the (mu) parameter in three types of correlations has clearly demonstrates the effects of blueshift and redshift at the height of (h<frac{r_A}{2}) and (h>frac{r_A}{2}). Therefore, our analysis has revealed a special height of (h=2r_A) (which corresponds to (delta =-,1)) at which the blueshift and the redshift effects might cancel out with each other. More interestingly, we find that the changes of correlations generated by the total gravitational frequency shift could reach the level of (<0.5%) within the satellite’s height at geostationary Earth orbits.