Analysis of Telescope Wavefront Aberration and Optical Path Stability in Space Gravitational Wave Detection

Abstract

Space-based gravitational wave detection programs, such as the Laser Interferometer Space Antenna (LISA) or Taiji program, obtain gravitational wave signals by measuring the change in the distance between three satellites by laser. The telescope is an important part of the measurement system, and its function is to transmit and receive laser signals. Due to changes in the space environment, the telescope will inevitably introduce additional dynamic aberrations, which will bring optical path errors to the inversion of gravitational wave signals. Taking LISA as an example, to achieve pm-level measurement accuracy at the detection frequency of 0.1 mHz–1 Hz, the stability requirements of the telescope are less than 1 pm/Hz^1/2. This paper theoretically deduces the aberration types that affect the telescope’s stability and conducts simulation analysis according to the actual phase demodulation method, which verifies the theory’s correctness. In addition, using this theory, it can be concluded that under the condition that the total size of the telescope aberration is determined to be stable, reducing the ratio of rotationally symmetric aberrations such as “spherical aberration” and “defocusing” among common aberrations can significantly improve the stability of the telescope. The conclusion guides the optical system design of LISA or Taiji.

1. Introduction

Space gravitational wave detection missions such as the Laser Interferometer Space Antenna (LISA) and Taiji program are all composed of three spacecraft to form an equilateral triangle with an arm length of millions of kilometers [1,2,3], and gravitational waves in the universe are observed by measuring the extremely small phase changes of the laser light during the propagation of the spacecraft. Therefore, measuring the phase change with high precision is one of the keys to detecting gravitational waves. Taking Taiji as an example, to successfully detect gravitational waves, the measurement noise of the detector needs to satisfy 8 pm/Hz^1/2 in the 1 mHz–0.1 Hz frequency band [4]. This ultra-low noise level places stringent demands on every part of a gravitational wave detector. As shown in Figure 1, the telescope system is a subsystem of laser interferometry, which is mainly responsible for local laser emission and reception. The far-field beam of the laser emitted by the spacecraft after propagating for millions of kilometers is intercepted and received at the entrance pupil of the telescope in the adjacent spacecraft, and then condensed by the reflective telescope and finally reaches the detector to interfere with the local Gaussian beam. During this process, the extreme environment of the telescope in orbit causes telescope deformation, lens spacing changes, and material refractive index changes, which will cause the far-field beam to introduce various aberrations of different sizes that will affect phase demodulation after interference with the local Gaussian beam. This effect will ultimately affect the detection of gravitational waves. Therefore, studying the coupling relationship between wave aberration and phase has a guiding significance for the design of the telescope wavefront index and aberration control.

Peter L Bender and CP Sasso et al. analyzed and studied the phase noise of the far-field wavefront caused by primary aberrations [5,6,7,8]. Eugene Waluschka analyzed the effect of the wavefront error on the received beam phase noise of distant spacecraft through numerical and ray tracing [9,10,11]. However, in the optical–mechanical–thermal integration analysis, the aberrations are often complex; high-level aberrations exist, and the dynamic changes of the aberrations will also have a greater impact on the phase. Zhao et al. made a more detailed analysis of the influence of the tilt jitter error in the interference process of a double beam (far-field beam and local Gaussian beam) on phase demodulation. At the same time, the relationship between the aberration and the phase brought by the telescope system was briefly analyzed, and it was pointed out that when the telescope aberration is within λ/20, the demodulation effect on the phase information can be maintained at 25 pm/nrad. However, there was no quantitative analysis of specific aberrations [12,13,14]. Therefore, a more complex and comprehensive analysis of the effects of aberrations and their dynamic changes on the phase is necessary to explore the influencing factors of telescope stability.

In this paper, the theoretical algorithm of phase demodulation of two-beam interference is firstly constructed, the aberration is fitted by the Zernike polynomial, and the influence of the wavefront aberration on the phase signal is deduced: after ignoring the higher-order terms, it can be concluded that the effect of non-rotationally symmetric aberrations (NRSAs) is negligible on an ideal detector; however, rotationally symmetric aberrations (RSA) such as spherical aberration and defocusing have great influences that cannot be ignored. At the end of the manuscript, the relationship between the dynamic aberration of the gravitational wave telescope and the laser phase demodulation on the quadrant photodiode (QPD) is simulated, which is consistent with the theoretical situation. Therefore, this conclusion shows that different aberrations have different effects on the detection of gravitational waves in practical engineering. When designing telescopes, not only is the wavefront index required, but the consideration of specific aberration allocation is also necessary. For the convenience of reading, the following articles will refer to the rotationally symmetric aberration and the non-rotationally symmetric aberration as RSA and NRAS for short.

2. Optical Model

In the gravitational wave detection interferometer, a far-field beam (according to the properties of Gaussian beams, the far-field beam can be approximated as a flat-top beam) is expanded by the telescope system and emitted by the laser in the adjacent satellite. After the telescope intercepts and converges, it interferes with the local Gaussian laser on the QPD. When properly adjusted (according to the DWS signal of the QPD), the center deviation of the double beams can reach below the micrometer level [14,15]. For the convenience of research, we assume that the center positions of the two beams are consistent with the incident angle. The local Gaussian beam on the detector can be written as (in Cartesian coordinates):

$$E(x,y,z{)}_{gauss}=A_1{e}^{-\frac{x^2+y^2}{w{(z)}^2}}{e}^{-i[k(\frac{x^2+y^2}{2R(z)}+{\phi }_1]+i2\pi f_{1}t}$$

(1)

where $f_1$ is the frequency of the Gaussian beam, $A_1$ is the amplitude. ${\phi }_1$ is the initial phase of the Gaussian beam, $w(z)=w_0\sqrt{1+{(z/z_r)}^2}$ is the beam width of the Gaussian beam, $R(z)=z(1+{(\frac{z_r}{z})}^2)$ is the radius of curvature of the equiphase surface of the Gaussian beam, and $z_r=(\pi w_0^2)/\lambda$ is called the Rayleigh criterion or confocal parameter. In the formula, in addition to the time phase factor $i2\pi f_{1}t$ and the initial phase ${\phi }_1$ due to optical element defects, the phase term $\frac{x^2+y^2}{2R(z)}$ is similar to the defocus. Its magnitude depends on the distance z between the interference position and the beam waist.

It can be seen from Formula (1) that when the propagation distance $z~\infty$, we can get $z>>z_r$ and $R(z)\approx z$. This means that the Gaussian beam can be approximated as a plane wave emanating from the point $z=0$ at the far field. Therefore, the far-field beams, separated by millions of kilometers in space gravitational wave detection, can be regarded as flat-top beams, which can be written on the detector as:

$$E(x,y,z{)}_{flat}=\{\begin{array}{cc}{A}_2{e}^{-ik(x\sin \alpha +z\cos \alpha )+i2\pi f_{2}t+{\phi }_2}\hfill & x^2+y^2\le R\\ 0\hfill & else\end{array}$$

(2)

where ${\phi }_2$ is the phase of the flat-top beam due to defects in the optics. For the convenience of derivation, in the subsequent theoretical derivation, the point to the angle $\alpha \to 0$, $A_2$, is the amplitude. Since the phase information is independent of the amplitude, we normalize the amplitudes of the two beams in the subsequent derivation.

The phase demodulation after two-beam interference can be obtained by simply multiplying and integrating complex amplitudes. In the theoretical analysis, the time factor is irrelevant to the phase and can be ignored, so we assume that in a circular detection area of radius $\rho$, the phase information can be written as:

$$\mathrm{int}={\displaystyle {\int }_s{E}_{flat}{E}^{*}{}_{gauss}}d{r}^2={\displaystyle {\int }_s{e}^{-\frac{x^2+y^2}{w{(z)}^2}}{e}^{ik\Delta \phi (\rho ,\theta )}}d{r}^2$$

(3)

where $\Delta \phi (\rho ,\theta )={\phi }_2-{\phi }_1-\frac{x^2+y^2}{2R(z)}$ is the static wavefront aberration of the interference of the two beams. According to space gravitational wave detection programs LISA and Taiji, the wave aberrations discussed in this paper are all smaller than the laser wavelength of 0.1064 μm, so the Taylor expansion of $e^{ik\Delta \phi (\rho ,\theta )}$ can be obtained and we can ignore the higher-order terms:

$$e^{ik\Delta \phi (\rho ,\theta )}\approx 1+ik\Delta \phi (\rho ,\theta )-\frac{1}{2}{k}^2\Delta \phi {(\rho ,\theta )}^2-\frac{i}{6}{k}^3\Delta \phi {(\rho ,\theta )}^3$$

(4)

Assuming that the detector is a complete circle, we perform a complete integration of the detector circle through the formula:

$$\Delta \phi =\arg (\mathrm{int})=\arg (\frac{\mathrm{Im}(\mathrm{int})}{\mathrm{Re}(\mathrm{int})})$$

(5)

We can obtain the complete expression of the phase signal by substituting Equations (3) and (4) into (5):

$$\begin{array}{cc}\Delta \phi & =\arg [{\displaystyle {\int }_0^1\rho e^{\frac{{\rho }^2}{w{(z)}^2}}}({\displaystyle {\int }_0^{2\pi }{e}^{ik\Delta \phi (\rho ,\theta )}d\theta })d\rho ]\hfill \\ & \approx \arg (s_0+i{s}_1-s_2-i{s}_3)\hfill \\ & \approx \frac{(s_1-s_3)(s_0+s_2)}{s_0^2}\hfill \end{array}$$

(6)

where ${\rho }^2=x^2+y^2$, and

$$s_n=\frac{1}{n!}{\displaystyle {\int }_0^1\rho e^{\frac{{\rho }^2}{w{(z)}^2}}}({\displaystyle {\int }_0^{2\pi }{k}^n\Delta \phi {(\rho ,\theta )}^{n}d\theta })d\rho$$

(7)

3. Aberration Theory Analysis

From the complete Formula (6) of the phase signal, it is not difficult to see that, in theory, the wavefront aberration between the far-field beam and the Gaussian beam will be considered the initial phase difference in gravitational wave detection. However, due to factors such as the space environment and satellite jitter, the final performance of the telescope’s mirror surface deformation and the increase in the distance between the primary and secondary mirrors can be equivalent to the changes in the telescope’s aberrations [16]. The telescope’s aberrations will always have subtle changes. Therefore, unstable aberration changes will be an inevitable source of errors in gravitational wave detection. In 2020, Zhao et al. used the first 25 items of the Zernike polynomial to fit the random aberration size of λ/20~λ/30 using the Monte Carlo algorithm on phase demodulation [12,13,14]. However, it only analyzed static aberrations and did not have a detailed decomposition of specific aberrations. This paper will start from the basic Zernike polynomial and derive the effect of its various parts on phase demodulation.

In actual situations, the situation of wavefront aberrations is extremely complicated, and the Zernike polynomial is a very good basis function for fitting the wavefront; for example, the spherical aberration and coma mentioned above can be easily fitted with it [17]. Its expression can be written as:

$$Z_i(\rho ,\theta )=\{\begin{array}{ll}\sqrt{2(n+1)}{R}_n^m(\rho )\cos (m\theta )& i \mathrm{is} \mathrm{even} \mathrm{and} m\ne 0\\ \sqrt{2(n+1)}{R}_n^m(\rho )\sin (m\theta )& i \mathrm{is} \mathrm{odd} \mathrm{and} m\ne 0\\ \sqrt{n+1}{R}_n^0(\rho )& m=0\end{array}$$

(8)

where $R_n^m(\rho )$ is a polynomial concerning the wavefront curvature, m and n are positive integers, and $n\ge m$. $\rho =r/r_0$ is the normalized radial coordinate. The index i is a polynomial-ordering number and a function of n and m. The specific expression $R_n^m(\rho )$ is:

$$R_n^m(\rho )={\displaystyle \sum _{s=0}^{(n-m)/2}\frac{{(-1)}^s(n-s)!}{s!(\frac{n+m}{2}-s)!(\frac{n+m}{2}-s)!}}{\rho }^{n-2s}$$

(9)

According to the above phase demodulation formula, the Zernike polynomial is divided into three parts for calculation. After substituting the first part $\sqrt{2(n+1)}{R}_n^m(\rho )\cos (m\theta )$ into (6), (7), we obtain (10):

$$\begin{array}{l}{s}_0^1=\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)\pi {\omega }^2\\ s_1^1=\frac{\sqrt{2}C\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)AR{\omega }^2\sin \left[2m\pi \right]}{m}\\ s_2^1=\frac{C^2\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)A^2{R}^2{\omega }^2\left(4m\pi +\sin \left[4m\pi \right]\right)}{2m}\\ s_3^1=\frac{2\sqrt{2}{C}^3\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)A^3{R}^3{\omega }^2\left(5+\cos \left[4m\pi \right]\right)\sin \left[2m\pi \right]}{9m}\end{array}$$

(10)

Substituting the second part, we obtain (11):

$$\begin{array}{l}{s}_0^2=\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)\pi {\omega }^2\\ s_1^2=\frac{2\sqrt{2}C\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)AR{\omega }^2\sin {\left[m\pi \right]}^2}{m}\\ s_2^2=2C^2\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)A^2{R}^2{\omega }^2\left(\pi -\frac{\sin \left[4m\pi \right]}{4m}\right)\\ s_3^2=\frac{16\sqrt{2}{C}^3\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)A^3{R}^3{\omega }^2\left(2+\cos \left[2m\pi \right]\right)\sin {\left[m\pi \right]}^4}{9m}\end{array}$$

(11)

Additionally, substituting the third part, we obtain (12):

$$\begin{array}{l}{s}_0^3=\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)\pi {\omega }^2\\ s_1^3=C\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)A\pi R{\omega }^2\\ s_2^3=\frac{1}{2}{C}^2\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)A^2\pi R^2{\omega }^2\\ s_3^3=\frac{1}{6}{C}^3\left(1-{\mathrm{e}}^{-\frac{1}{{\omega }^2}}\right)A^3\pi R^3{\omega }^2\end{array}$$

(12)

where $A=k\sqrt{1+n}$ and C is the correlation coefficient of each part of the Zernike polynomial. It is worth noting that since the integral interval is 0~2π, the s₁ and s₃ corresponding to the first two parts of the Zernike polynomial becomes 0 after ignoring high-order terms, and the third part has nothing to do with the integration interval because there is no θ. Therefore, theoretically, after the high-order term is ignored, the part of the Zernike polynomial that affects the subsequent phase demodulation is mainly the part of $m=0$; this part is the RSA mentioned above.

4. Experimental Simulation

4.1. Phase Demodulation Algorithm of QPD

According to the nature of space gravitational wave detection, the propagation distance of the beam in the two satellites is in the order of millions of kilometers, and it is difficult to determine the direction of the beam. As the phase detector used in this paper is QPD of which structure is shown in Figure 2, it can acquire the DWS (differential wavefront sensing signal) at the same time as acquiring the phase signal, and by analyzing the DWS signal, a suitable algorithm is used to compensate for the pointing error. This part is described in detail in [18], so it will not be explained in this paper.

In the actual phase demodulation, the detector can only detect the light intensity information, and the phase information cannot be obtained through the inversion of the real part and the imaginary part in the formula according to the complex amplitude integration in the theoretical simulation. At the same time, due to the Doppler effect between satellites [19], the frequency difference between the local Gaussian beam and the far-field beam is about ±15 MHz. Therefore, the method we use is the quadrature demodulation algorithm shown in Figure 3 below. The time-varying optical signal acquired on the phase meter can be written as:

$$I=\overline{I}[1+c\cos ({\omega }_{het}t+\Delta \phi )]$$

(13)

where $\overline{I}=\frac{1}{2}(I_{max}+I_{\min })$ is the average intensity of the resulting light field, $c=\frac{I_{max}-I_{\min }}{I_{max}+I_{\min }}$ is the signal contrast, ${\omega }_{het}={\omega }_{gauss}-{\omega }_{flat}=2\pi (f_1-f_2)$ is Heterodyne frequency, and $\Delta \phi$ is the phase difference between the two beams. The subsequent demodulation steps are as follows:

(1)

The received interferometric light signal is converted into an electrical signal by a phase meter.

(2)

We use the NCO (numerical oscillator) to send out sine and cosine signals and multiply the interference signal. Then, we get Formula (11). It should be noted here that in actual detection, the heterodyne frequency is unknown due to the existence of the Doppler effect, and a subsequent feedback mechanism makes the final NCO send out the same sine and cosine signals as the heterodyne frequency. However, for the convenience of research and because the feedback processing has nothing to do with phase demodulation, we assume that the heterodyne frequency is known.

$$\begin{array}{l}\mathrm{int}1={\displaystyle {\int }_0^{2\pi }\cos ({\omega }_{het}t)}\cdot I({\omega }_{het}t)d({\omega }_{het}t)=\overline{I}c\pi \cos (\Delta \phi )\\ \mathrm{int}2={\displaystyle {\int }_0^{2\pi }\sin ({\omega }_{het}t)}\cdot I({\omega }_{het}t)d({\omega }_{het}t)=\overline{I}c\pi \sin (\Delta \phi )\end{array}$$

(14)

(3)

Finally, the phase information we get can be written as Formula (15):

$$\Delta \phi =\arctan (\frac{\mathrm{int}2}{\mathrm{int}1})$$

(15)

4.2. Aberration Simulation Analysis

After the above demodulation model is established, the simulation analysis of the experiment can be carried out in combination with the beam model of the first part. To facilitate understanding, we will start from the primary phase to verify the theory. Among the common aberrations, the spherical aberration and the coma are typical of RSA and NRSA, their shapes are shown in Figure 4 (spherical aberration is (a), coma is (b)). Therefore, we use Zernike polynomials to fit spherical aberrations and coma that vary from λ/100 to λ, then substitute them into the above model to get the resulting graph as follows:

From the results in Figure 5, it can be seen that as the aberration increases, the influence of spherical aberration on the phase is greater than that of the coma. Under the same aberration change, the phase change corresponding to spherical aberration is about five times that of a coma. This is consistent with the above theory.

In actual telescope systems, aberrations are extremely complex, including not only primary aberrations but also advanced aberrations. Therefore, we simulated a more general aberration composition, and to see the effect more clearly, we used the Monte Carlo algorithm to randomly assign the first 15 terms of the Zernike polynomial and fix its total aberration size at λ/60. Then, we obtained a comparison chart (Figure 6) of removing the non-rotationally symmetrical aberrations and not removing them.

The abscissa in Figure 6 represents the Monte Carlo random number, and the ordinate represents the size of the demodulation phase. It can be seen that the red part (including rotational symmetry aberration) has a greater influence than the blue part (not included). Their variances are V_r = 1.9560 × 10⁻⁰⁶, V_b = 2.4549 × 10⁻¹⁴. Therefore, it once again confirms the above theory in this paper. The effects of RSA and NRSA on phase demodulation in space gravitational wave detection are different, and it is obvious that RSA is more important.

4.3. Telescope Stability Analysis

In the LISA or Taiji space gravitational wave detection projects, the evaluation criteria of the telescope’s stability are different from those of traditional space telescopes. Parameters such as the wavefront index and frame deformation are unified as the PSD (power spectral density) of the optical path noise. Taking LISA as an example, the PSD of the optical path noise of the telescope is required to be in the range of 0.1 mHz~1 Hz:

$$PSD\le 1\frac{\mathrm{pm}}{\sqrt{Hz}}\sqrt{1+{(\frac{f_c}{f})}^2}$$

(16)

To more intuitively feel the impact of dynamic phase changes on the evaluation of telescope stability, we simulated 20,000 random samples using the Monte Carlo algorithm. Therefore, we obtained the power spectral density in the frequency range of 0.1 mHz~1 Hz.

It can be seen from Figure 7 that whether it is a noise signal or its PSD, after removing the rotational symmetry aberration, the size of the noise signal can be reduced by about 5 orders of magnitude, while the power spectral density signal is nearly 6 orders of magnitude lower. Additionally, the power spectral density of the system meets the LISA requirements in the range of 0.1 mHz~1 Hz when the rotational symmetry aberration is eliminated. It is worth noting that the variation trend of noise and power spectral density is consistent in the two random cases, and the maximum value of power spectral density is around 10⁻². The reason may be related to the phase demodulation algorithm.

In addition, to quantitatively analyze the degree of influence, we analyzed the power spectral density of the aberration corresponding to different aberration sizes. As shown in Figure 8, we analyze the mean value of the power spectral density of the phase noise corresponding to the inclusion of RSA (orange) and not the inclusion of the RSA (blue) with an aberration size of λ/80-λ. It can be found that as the aberration decreases, the noise level also decreases rapidly. However, compared with the completely random aberration, the noise level drops by nearly 5 orders of magnitude after removing the rotationally symmetric part, and the ratio does not change significantly with the reduction of the overall aberration. Hence, we can think that rotational symmetry aberration plays a dominant role in phase noise.

In the field of data analysis, box plots are a good way to reflect data stability. Next, as shown in Figure 9 (box plot), we fix the overall size of the aberration and redistribute the weight of the rotationally symmetric part in the random aberration. The abscissa represents the proportion of the rotationally symmetric part of the aberration to the total aberration; the ordinate is the power spectral density. In the figure, we have removed outliers, and we kept the minimum, first quartile (Q1), second quartile (Q2), third quartile (Q3), and maximum, where Q2 is also known as the median. It is obvious that Q1, Q2, Q3, and the maximum of the phase noise shift as the ratio increases. This shows that as the ratio increases, the degree of instability of the system also increases, which is fatal to space gravitational wave telescopes.

5. Discussion

Telescopes are one of the important components in space gravitational wave detection missions. Both LISA and the Taiji mission put forward more stringent requirements for telescope indicators. However, as the space environment changes, the telescope will generate additional dynamic aberration, and the aberration will generate phase noise through the interference of the double beams in the gravitational wave measurement, thereby reducing the stability of the telescope.

In this article, we first built the interference model between the far-field beam and the local beam in space gravitational wave detection and theoretically deduced the inversion method of its phase; then, we used the general form of the Zernike polynomial to deduce the aberration form that plays a dominant role in the phase effect. At the same time, to better simulate the phase demodulation method in the actual situation, we deduced the phase demodulation method on the four-quadrant detector.

In terms of simulation, we first simulated the effects of different sizes of spherical aberration and coma (common aberration) on the phase and preliminarily found that the RSA plays a dominant role in the phase noise. Therefore, a Monte Carlo random algorithm was used to simulate the case where the fixed aberration size is λ/60, and the influence of the RSA on the phase demodulation was compared with and without the exclusion of the RSA. The result was consistent with the theoretical derivation. In demodulated phase noise, the variance aberration is about 8 orders of magnitude.

To further analyze the influence of the rotationally symmetric part of the aberration on the stability of the telescope, we took the LISA requirements for the stability of the telescope as an example, established a correlation model between the aberration and the phase noise power spectral density, and simulated the noise power spectral density corresponding to aberrations of different sizes. Then, it was further analyzed that when the aberration size is fixed at λ/30, the rotationally symmetric part accounts for different proportions of the corresponding power spectral density, which finally shows that the stability of the system can meet the requirements when the proportion is less than 1%.

In general, the above analysis and conclusions point out that for the detection of gravitational waves in space, the additional dynamic aberration caused by the telescope’s environmental influence will greatly affect the stability of the telescope, and the rotational symmetry of the aberration plays an important role. Therefore, when designing a telescope, the design of its structure and the allocation of tolerances should minimize the generation of rotational symmetry aberrations, which has guiding significance for the design of optical systems in LISA or Taiji missions.