Quantum optical phase estimation is a critical task in many applications such as quantum imaging [1–3], quantum sensing [4–7], and gravitational wave detection [8,9]. To date, most works have focused on the estimation of a constant phase φ, in which a Mach–Zehnder interferometer (MZI) is the most commonly used device [10,11]. The precision of estimation is limited by the shot noise when the classic resources are used. This limit is often called the standard quantum limit (SQL), Δφ∝1/N, where N is the average number of photons in the probe state [12,13]. Many efforts have been taken to improve the precision. Most of them focus on utilizing nonclassical states to reduce the quantum noise, such as the squeezed states and entanglement states [14–16]. It has been shown that the maximally entangled number state (N00N) is the optimum probe state to reach the Heisenberg limit (HL) Δφ∝1/N [17,18]. Moreover, for the constant phase estimation, the variance of estimation (1/vN or 1/vN) will decrease indefinitely as the number of measurement v increases.

However, it is not enough to just estimate the constant phase because many signals of interest in the real world are time varying and stochastic [19–23]. Thus, how to estimate such a time-varying phase with high precision is of practical importance. We assume that φ(t) is the phase to be estimated. It can be treated as a constant in ti≤t<ti+dt if the dt is small enough, i.e., φ(t)=φi, so the phase can be discretized as (φ0,φ1,…,φi,…,φn−1,φn). There is a set of observations (r0,r1,…,ri,…,rn−1,rn) used to estimate the {φi}, where ri is the observation at time i. Compared to the estimation of constant phase, ⟨(φi−φi±s)2⟩ between the phases at two different times (i,i±s) increases as s increases. Therefore, the correlation between ri±s and φi decreases as s increases, and the number of observations that can be used to improve precision is limited. So there is a limited precision for the stochastic phase estimation, even for infinite measurement time. The mean square error (MSE) in the estimation of a stationary Gaussian stochastic phase with a power-law spectrum κp−1/(ωp+λp) using coherent states scales as (κ/N)(p−1)/p, which is called stochastic SQL. Here N is the photon flux [24–26]. Similar to constant phase estimation, there is a stochastic Heisenberg scaling for the estimation of this stochastic phase, which scales as (κ/N)2(p−1)/(p+1) [24–26]. Previous works show that the stochastic SQL can be surpassed with nonclassical resources such as the squeezed light with an adaptive quantum smoothing technique or canonical phase measurement [27–32]. However, practical methods to achieve quantum enhancement in the stochastic phase estimation still remain largely unexplored.

In this paper, we propose a method to estimate the stochastic phase using a new measurement device, the SU(1,1) interferometer. Such device, also known as the nonlinear interferometer (NLI), was original proposed by Yurke and experimentally demonstrated with a signal-to-noise ratio (SNR) surpassing that of the MZI [33–43]. By combining the NLI with adaptive feedback technique, we show that MSE of the estimation is reduced over a range of parametric amplifier gain G compared to that of MZI. For a fixed photon flux, there is an optimal G that minimizes the MSE and maximizes the precision. In particular, with the optimal G, the precision of our scheme surpasses the stochastic SQL and achieves stochastic Heisenberg scaling asymptotically.

The schematic diagram of the estimation of a stochastic phase with NLI is shown in Fig. 1, in which the NLI contains two parametric amplifiers (PAs). Two input modes of the first PA are injected with a coherent state |α⟩ and a vacuum state. The phase φ(t) to be estimated is imposed to one arm of the interferometer. One of the output modes is measured with the homodyne measurement. The measured results after a displaced operation yield photocurrent r(t). The phase Φ(t) in the other arm and the phase θ(t) of the local oscillator are adaptively controlled based on φf(t), which is estimated from r(s) for all the regions s<t. In the NLI, the first PA plays the role of beam splitting. If we define c^in,d^in to be the annihilation operators of the two inputs and C^,D^ to be the annihilation operators of the outputs, the relation of the input/output of the PA can be written as C^=Gc^in+gd^in†, D^=Gd^in+gc^in†, where the gain G of the PA can also be parametrized by hyperbolic functions [44,45] and G2−g2=1 [36]. The second PA, which has the same gain G, plays the role of recombination, so the complete input/output relation of the NLI is

When we perform a homodyne detection at the output mode dout, and the homodyne detection result is added by 2Gg|α|φf(t), the photocurrent can be approximately represented as (Appendix A) Here we have adaptively controlled the feedback phase and the phase of the local oscillator to be Φ(t)=−φf(t)−π, θ(t)=φf(t)+π/2. This adjustment makes each measurement most sensitive and maximizes the phase information obtained. n(t) is the normalized Gaussian white noise from the homodyne measurement, which satisfies ⟨n(t)n(s)⟩=δ(s−t). Moreover, the photocurrent has Gaussian stationary statistics, and σf2=⟨[φ(t)−φf(t)]2⟩ is stationary MSE. We have defined the photon flux in the interferometer to be |β|2=(G2+g2)|α|2.

According to the time span of the observations in the estimation, the time-varying phase estimation can be divided into three cases. Prediction: the future phase φi+m is estimated with observations (r0,r1,…,ri). Tracking: the current and previous observations (r0,r1,…,ri) are used to estimate the current phase φi. Smoothing: the measurement results beyond the time i are also used to estimate phase φi. To analyze these three estimation problems simultaneously, we introduce a general desired signal d(t)=φ(t+ε), which is estimated with photocurrent r(s≤t). Here ε can be any real number, and the three kinds of phase estimation can be defined according to the value of ε. Based on the measurement photocurrent r(t), the desired signal can be estimated as where ho(t,τ) is the impulse response function, which represents the output at time t if the input at time τ is an impulse, and it can be marked as hop(t,τ),hot(t,τ),andhos(t,τ) for prediction, tracking, and smoothing, respectively. The optimum linear processor ho(t,τ) for minimizing the MSE ξ(t)=⟨[d(t)−df(t)]2⟩ satisfies [46] which is the Wiener–Hopf equation and Kdr(t−η)=⟨d(t)r(η)⟩, Kr(ϵ−η)=⟨r(ϵ)r(η)⟩. The correlation function only depends on the time difference due to the fact that received photocurrent and desired signal are jointly stationary and time invariant. We can obtain the optimum linear impulse response function ho by solving the Wiener–Hopf equation, and the MSE of the phase estimation can be calculated at the same time.

In the MZI case, the splitting and recombination of light are accomplished by 50:50 beam splitters (BSs). The relation of the input/output of the beam splitter is C^=1/2(c^in+id^in) and D^=1/2(ic^in+d^in). Similar to NLI, we set the two feedback phases Φ(t)=φf(t) and θ(t)=φf(t)+π for the most sensitive estimation. The photocurrent can be calculated as r(t)≈|β|φ(t)+n(t) (Appendix A). Here we should note that the two input modes of the first BS are injected with a coherent state |β⟩ and a vacuum state, which makes the photon number flux inside both interferometers the same. From the two photocurrents, we can derive the relation SNRNLISNRMZI=4G2(G2−1)(2G2−1)[2G2(G2−1)σf2+1], where SNRNLI and SNRMZI are the SNRs of these two interferometers, respectively. Figure 2 shows the ratio of the two SNRs (blue surface). As G increases from 1, the SNR of the NLI increases with G, which agrees with the previous analysis of NLI [36]. Moreover, the two parameters (G,g) in the relation of the input/output of the PA are asymmetric. If the gain G is close to 1, for the fixed photon flux |β|2, the part of signal of Eq. (2) is close to zero and the SNR of the NLI is smaller than that of the MZI. When G exceeds a certain threshold, the SNR of the NLI surpasses that of the MZI. For finite σf2, the further increase of G will reduce the SNR. This result can be understood from Eq. (2): when G is large, the signal term increases linearly with G while the noise term increases quadratically with G. Since the MSE in the estimation of time-varying phase cannot be arbitrarily small, we expect there is an optimal G for stochastic phase estimation, which is different from the case of measuring a constantphase.

As an example, we consider the situation that the time-varying phase φ(t) to be estimated follows an Ornstein–Uhlenback stochastic process, which can be found in many practical physical processes and is defined by [31] Here λ−1 is the correlation time of φ(t). dV(t) represents the Wiener process, which satisfies ⟨dV(t)dV(s)⟩=δ(s−t)dt. κ is the magnitude of the Wiener noise. The expectation value of φ(t) is 0 and its statistics are stationary, which means the correlation between the phases at two different times only depends on their time difference. Moreover, the spectral density spectrum of φ(t) is Sφ(ω)=κ/(ω2+λ2). In this situation, Eq. (4) can be solved with the Wiener technique, and the Fourier transform of the optimum linear response function ho is (Appendix B) where P=4G2g2|α|2, Λ=Pκ/Nλ2, and N=2G2g2σf2+1. Equation (6) shows that the optimum linear processor is a low-pass filter with a cutoff frequency of λ1+Λ in the case of ε=0. When this optimum linear processor is used to estimate the phase, the minimum MSE of the estimation is (Appendix C) where The stochastic phase is stationary and the optimum filter is time invariant, so the minimum MSE is a constant. To calculate the MSE analytically, we divide the discussion into three cases according to the value of ε: (i) ε=0, (ii) ε>0, and (iii) ε<0. When ε=0, d(t)=φ(t). This is the phase-tracking case, and ξ=σf2. The integral result ξ of Eq. (7) is still implicit because Λ is a function of σf2. After solving the implicit result, the minimum MSE of tracking is (Appendix C) Similarly, the MSEs of the other two cases can be calculated as (Appendix C) where ε>0 stands for the prediction of the future phase with current measurement outcomes and ε<0 is the case of smoothing. Before investigating the enhancement of phase estimation with NLI, we set the precision of the Ornstein–Uhlenback stochastic phase estimation with a coherent state and the MZI as the classical limit. To compare the MSEs of the two types of interferometers, we make the photon number flux inside them equivalent, i.e., NMZI=NNLI=|β|2. In this case the MSE of phase estimation with the MZI can be written as (Appendix C) where Λ1=|β|2κ/λ2. To investigate the effect of the gain G on the MSE σf2, we consider the stochastic phase tracking with a fixed photon number flux |β|2=1.0×107s−1, κ=1.0×104rad/s, and λ=1.0×105rad/s, and we vary G2 from 1.1 to 50. Figure 3 shows that the MSE σf2 can be reduced using the NLI compared to the classical limit with MZI. There is an optimal G2=7.4 that gives the minimum MSE. For the feasible gain (G2=2.28) with current techniques [47], we can already obtain 40.4% improvement in the MSE of phase tracking. This is expected from the analysis of the SNR. Here the optimal degree of gain depends on the photon number flux |β|2, κ, and λ.

Figure 4 shows the MSE ξ as a function of λε according to Eqs. (9)–(11), where the horizontal axis is the ratio between ε and the correlation time of φ(t). Here we have set the parameters κ=1.0×104rad/s, λ=1.0×105rad/s, G2=7.4, and |β|2=1.0×107s−1. We can conclude three key implications from Fig. 4. First, the MSE becomes large with the increase of ε, and the smallest error is achieved with the smoothing. When ε is close to the correlation time, the MSEs tend to be the mean square variation of the stochastic phase κ/2λ, i.e., we cannot predict the phase away from coherence time. Second, the MSE of phase estimation with tracking is nearly two times of smoothing for both kinds of interferometers. Third, the MSEs of all cases are reduced significantly below the classical limit (red line in Fig. 4) when we use the NLI.

So far we have shown that the stochastic phase estimation can be enhanced with NLI. In the following, we will demonstrate that the NLI achieves the Heisenberg scaling asymptotically. For a fixed photon flux |β|2, the SNR of the measurement photocurrent Eq. (2) is SNRNLI=4G2(G2−1)|β|2[(2G4−G2)σf2+1](2G2−1). There is an optimal degree gain Go maximizing the SNR, and the minimum MSE of phase tracking can be calculated as σf2≈1/2Go4. When Go2≫1, |β|22Go2≫λ, |β|2≫κ, and Go4κ≫λ, the optimal gain Go meets the relationship Go2≈(|β|2κ2)1/3/22/3κ and we can obtain the tracking MSE (Appendix D) Substituting this expression into Eq. (10) yields the MSE of smoothing which means the MSE of stochastic phase estimation with NLI can achieve the stochastic HL scaling [24–26]. Figure 5 shows the optimal smoothing MSE in the two kinds of interferometers and canonical measurement for different mean photon flux, which varies from |β|2=109s−1 to |β|2=1010s−1. It can be seen that the phase estimation with NLI has an enhancement on scaling compared with the classical limit ξ≈12(κ|β|2)1/2 using MZI, and the smoothing MSE can reach the stochastic Heisenberg scaling O[(κ|β|2)2/3] [25]. Moreover, using the NLI we can surpass the minimum MSE of canonical measurement, which is 45(κ|β|2)2/3 [24]. It is worth noting that the MSE of tracking is nearly two times of smoothing for both interferometers. Therefore, the phase tracking can also reach the stochastic Heisenberg scaling. However, the MSE of stochastic phase prediction increases with the increase of ε, and we cannot predict the phase away from coherence time. So the prediction cannot reach the stochastic Heisenberg scaling when ε is relatively large.

In summary, we have proposed the stochastic optical phase estimation with an SU(1,1) interferometer. We find that a suitable range of parametric amplification gain can enhance the estimation, and there is an optimal gain minimizing the MSE, which is different from the estimation of constant phase [34,36]. Moreover, compared with the classical limit with MZI, the MSEs have significant reduction for prediction, tracking, and smoothing simultaneously under the same photon number flux inside the interferometers if we optimize the parametric amplifier gain. At last, we can achieve the stochastic Heisenberg scaling and surpass the minimum MSE by using the canonical measurement. These results highlight the advantages of the SU(1,1) interferometer in stochastic optical phase estimation and provide a new avenue for practical quantum metrology.