© The Authors 2023

1 Introduction

Since the Michelson stellar interferometer was mounted on the Hooker telescope, followed by its successful measurement of the diameter of Betelgeuse in 1920 (Michelson & Pease 1921), ground-based interferometry has been widely used for optical, infrared, and radio astronomy (Beckers et al. 1990; Colavita & Wizinowich 2000; ten Brummelaar et al. 2005; ALMA Partnership 2015). An image of the sky was also reconstructed based on the Van Cittert-Zernike theorem (Born & Wolf 1999). The Fourier transform of the filled U-V plane offers a two-dimensional (2D) image by rotating the baseline and changing its length.

Bracewell (1978) introduced the concept of nulling inter-ferometry to the search for exoplanets around nearby stars by introducing a π phase shift to one of the beams of a two-beam interferometer. When the observed sky consists of a host star and multiple planets, this concept generates a sin² ($\frac{\pi }{\lambda }B\cdot \theta$) fringe pattern that can nullify the host star at the centre of the field-of-view and transmit light from an off-axis point source, such as a planet, where À is the observing wavelength, B is the baseline vector, and θ is the position vector on the sky. Rotating the baseline of the interferometer also modulates the signal of the off-axis source as a function of time, which can be leveraged for signal extraction purposes. Followed by the proposal of mid-infrared nulling interferometry, Angel et al. (1986) noticed that the mid-infrared wavelength range is useful for the characterization of temperate Earth-like planets. This is because of the relatively low contrast between the planet and its host star compared to that observed in the visible wavelength range. In addition, CH₄ and O₃ are atmospheric biosignatures that have strong absorption bands in the same wavelength range (e.g., Des Marais et al. 2002; Fujii et al. 2018). A combination of these two studies led to the construction of a concept for remotely measuring the activity of primitive life on distant planets through detecting a variety of CH₄ and O₃, named Darwin (Léger et al. 1996). Darwin was an ESA-led concept and, in parallel, similar activities were being undertaken on the US side, in the context of Terrestrail Planet Finder-Interferometer (TPF-I; Lawson et al. 2008).

Angel & Woolf (1997) used a cross-correlation technique to efficiently find the modulated signal while rotating the baseline, leveraging Bracewell’s idea. A nuller consisting of four apertures was also introduced to obtain a fourth-order null of the host star. Furthermore, Mennesson & Mariotti (1997) proposed five collectors to suppress modulation of the exozodiacal light during baseline rotation, keeping the fourth-order null. Instead, Velusamy et al. (2003) mentioned the advantage of a dual Bracewell interferometer consisting of two equivalent second-order nullers. Its use allows for the ambiguity of the planetary positions to be overcome with a phase chop. We note that a phase shifter (π/2 for nulling interferometers) is introduced to one of the two beams, and the two states are formed by inserting or removing the phase shifter. In addition, the subtraction of the two chopped states can separate symmetric components, including stellar leakage and background light, from the off-axis point sources. Finally, both the TPF-I and Darwin mission concepts have favored a dual Bracewell interferometer (Cockell et al. 2009). Although TPF-I and Darwin were anticipated to detect and characterize the thermal emissions from Earth-like planets for the first time, they were postponed indefinitely due to the technical difficulties involved.

However, stellar leak is more sensitive to the optical path difference (OPD) and low-order aberrations in the second-order null than in the fourth-order null (e.g., Hansen et al. 2023a). Initially, Lay (2004) quantified the systematic errors generated from the fluctuation of the null depth, which could obscure the modulated planetary signal over the baseline rotation. The ideal phase chop technique can successfully remove most of the systematic noise errors, and only the first-order phase error and the cross-term of phase and amplitude errors remain in the demodulated signal. To identify an Earth-like planet around a Sun-like star at 10 µm, the OPD and amplitude perturbations need to be stabilized to 1.5 nm and 0.1%, respectively (Lay 2004). Because observing such a planet with a signal-to-noise ratio (S/N) of 7 requires an integration time of a few days, according to preliminary analyses for the LIFE mission presented in Dannert et al. (2022), the 1.5-nm OPD stability requirement holds for the same period. This imposes strict requirements on the formation flight and the optical beam transport and combination system. Lay (2006) later proposed stretching the aspect ratio of the rectangular four collector array of the Double Bracewell interferometer to remove instrumental noise induced by the systematic effects from the data to mitigate the requirements by a factor of 10. The method utilized different behaviors between the planet and instability noise achieved by stretching the ratio between the nulling and imaging baseline to 1:6. While this relaxes the requirements on the null stability, the stretching of the baselines requires more fuel consumption compared to baseline rotation.

Instead of using the modulated signals of off-axis objects during baseline rotation, Matsuo et al. (2011) proposed a method for estimating the positions of off-axis objects and obtaining their spectra from a few baselines, focusing on the modulation of the signal along the wavelength domain. This method requires only relative stability among the wavelength channel across the observed spectrum, rather than the stability of the null depth during the baseline rotation. In addition, when the number of baselines is larger than that of the detectable objects in the field of view of the interferometer, we can effectively separate planetary signals from long-term fluctuations caused by systematic effects. The extended method optimizes for continuously rotating the baseline instead of fixing baselines, which could further mitigate the stability requirements of a space-based nulling interferometer. The present study is complementary to developments of the formation flying interferometry (e.g., Hansen et al. 2023b; Matsuo et al. 2022) and ground-based nulling experiments (e.g., Ertel et al. 2020; Ranganathan et al. 2022). These efforts provide a support for the Large Interferometer for Exoplanets (LIFE), which represents a science theme that was recognized as one out of three potential science themes for a future L-class mission in the Voyage 2050 of the European Space Agency. Based on the heritage of Darwin and TPF-I, while leveraging the most recent scientific and technological developments, LIFE will directly detect and characterize the thermal light from habitable planet candidates. LIFE is poised to detect 25–45 terrestrial planets in the habitable zone around nearby F-, G-, K-, and M-type stars under the most conservative assumptions (Quanz et al. 2022; Kammerer et al. 2022).

In the following, we propose the phase-space synthesis decomposition (PSSD) method for extracting the planetary signal. This method could mitigate the rigorous requirements imposed on the nulling space interferometer, which is complementary to ongoing technological demonstrations for LIFE. Section 2 provides a brief overview of PSSD and its mathematical explanation. We perform a numerical simulation to investigate the feasibility of the method using the LIFE simulator (Dannert et al. 2022) in Sect. 3. Section 4 gives a discussion of the limitation of PSSD, along with its advantages and disadvantages. We present our main findings in Sect. 5.

2 Concept

This section provides an overview of PSSD and then introduces its procedures from planet detection to spectral characterization. The analytical equations constructed for PSSD offer a detailed explanation of the processes of planet detection and characterization.

2.1 Overview

The PSSD method is split into two processes: (1) a search for the planetary signal and (2) measurement of the planetary spectrum. Both steps require a continuous baseline rotation and the same operation as the previous method to extract the modulated signal through the cross-correlation or the maximum likelihood of the data obtained while rotating the baseline (e.g., Angel & Woolf 1997; Dannert et al. 2022). However, PSSD differs from the previous method in how the planetary positions are reconstructed. The earlier process transforms interferometric signals collected by spinning the baseline into the planet position by fitting the modulation in both the azimuth and wavelength domains simultaneously. However, focusing on the fact that the spectrum of a G-type star is smoothly distributed at mid-infrared wavelengths (e.g., Husser et al. 2013), a one-dimensional (1D) image parallel to the baseline is first formed by correlation of the interferometric signal only along the wavelength domain, which basically has the same characteristics as a Fourier transform of the signal in the same direction (Matsuo et al. 2011). Because the wavelength dependence of the stellar leak is largely different from that of the planetary signal, the partial correlation only along the wavelength domain can decompose the stellar leak and that signal. Instead of a cross-correlation of the signal in the azimuth domain after spinning the baseline and summing over the 1D images, PSSD transforms a set of 1D positional information into 2D positional information. Thanks to the partial correlation only in the wavelength domain, the effects of a long-term systematic error on image reconstruction can be avoided.

The method displays the advantages of the two previous methods, namely, the cross-correlation method (Angel & Woolf 1997) and Fourier transform of the signal along the wavelength domain (Matsuo et al. 2011). While the cross-correlation method increases the S/N as much as possible, the latter efficiently decomposes the stellar leak and planetary signal. The PSSD method combines the two previous approaches by employing a local cross-correlation of the signal only along the wavelength domain, instead of the full cross-correlation in both the wavelength- and time-domains. Generally, local (or segmented) cross-correlation is used when cross-correlating data in smaller segments. We can directly compare small sections of two arrays of data by cross-correlating corresponding segments, allowing for a more localised analysis. This method is very useful when analyzing complex astronomical phenomena with fluctuations or patterns in various portions of compared signals (e.g., Kovačević et al. 2018).

Thanks to the combination of the two previous methods, PSSD provides three key advantages in terms of planet detections. First, PSSD could shorten the required stability duration from a few days to a few minutes. Second, PSSD could also mitigate the impact of the limited number of baselines on search for the planetary signal. Third, PSSD could have robustness against a larger OPD fluctuation. Utilizing the advantages of the planet detection process, PSSD also offers a method for extracting the planetary spectrum embedded in the stellar leak.

Regarding the first advantage, the required stability duration could be shortened from a few days to a few minutes because we do not use the correlation of the planetary signals collected while rotating the baseline. This is equivalent to the period for obtaining the two-phase chop states. We note that the period of switching between the two-phase chop states is determined such that the slow change of the background can be fully sampled (Absil et al. 2003). PSSD only requires relative stability along the wavelength (i.e., among the spectral data), rather than the stability of signals received while turning the baseline. The continuous and wide wavelength range obtained from space, such as 4–18.5 µm for the LIFE observatory, realizes this alternative approach. We note that the particular type of object (e.g., Jovian planet or terrestrial planet) orbiting the host star is unclear in the planet detection phase because the light is integrated over the entire wavelength range (4–18 µm for the LIFE mission) in this phase. Regarding the second advantage, because PSSD reconstructs a 1D image from one imaging baseline, 2D positional information can be extracted from fewer baselines. In other words, PSSD offers more resistance against a limited number of data collected during baseline rotation than the previous cross-correlation technique. In terms of the last advantage, a fluctuation of OPD during baseline rotation is not correlated with a modulation of the planetary signal along the wavelength domain. Instead, the OPD error contributes to the observing data as a noise because the stellar leak is inversely proportional to approximately the fourth power of wavelength. The PSSD method is able to detect the planet light unless the modulation of the planetary signal is embedded in the stellar leak. Thus, it is more robust against a large OPD error in terms of planet detection.

Next, the planetary spectra are derived based on the positional information of the planets. Because we typically calculate the modulation of planet light while rotating the baseline based on the information of the estimated planet position, the planetary light for each spectral channel can be extracted by fitting the data through the singular value decomposition (SVD) method. However, because the modulation of the planetary signal during the baseline rotation is used for the reconstruction of the planetary spectrum, the reconstructed spectrum is more affected by a long-term systematic error compared to the planet’s detection. As a result, the large contrast between the stellar leak and a temperate planet at short wavelengths prevents us from precisely reconstructing the planetary spectrum in the same wavelength range (Sect. 4.1). Before applying the data to the SVD method, if we fnd that the stellar leak is much brighter than the planetary light, the stellar leak has to be subtracted from the data. On the other hand, it is difficult to measure the OPD change because the number of available photons is very limited at the nulled output.

Here, we find that the stellar leak induced by the systematic OPD error could be measured from the data at short wavelengths. Because warm and temperate planets are much fainter than the stellar leak in the short wavelength range, only the stellar leak mainly contributes to the data. Since there is no strong chromatic aberration in the optical system, thanks to reflectors constructing the optical system of LIFE, the stellar leak may be expressed as a function of wavelength. For example, when the stellar leak is induced by the OPD error, the stellar leak is inversely proportional to approximately the fourth-power of wavelength. Since there are a large volume of data collected during baseline rotation, we could estimate the wavelength dependence of the stellar leak using a simple model, such as an exponential function. We note that modeling the wavelength dependence of the stellar leak was already performed in the data reduction pipeline of GRAVITY (e.g., GRAVITY Collaboration 2020). Once the wavelength dependence is derived, the modeled stellar leak could be extrapolated to a longer wavelength range. The modeled stellar leak is subtracted from the data and is not, in principle, contaminated in the reconstructed spectra. Because the stellar leak is much weaker at the wavelengths longer than 10 µm, the planetary spectrum is less affected even if there exists a chromatic aberration.

We note, however, that if a Jovian planet close to its host star exists in the observing object, the planet is brighter than the stellar leak even at short wavelengths. Subtracting the bright planet from the data is required for modeling the wavelength dependence of the stellar leak at short wavelengths. Thanks to the long imaging baseline, an inner planet can be spatially resolved from the host star. Because both the position and spectrum of the inner planet are obtained through PSSD, we can estimate how the planetary signal is modulated during baseline rotation and subtract it from the data.

We also need to emphasize that PSSD is validated only for objects having smoothed spectra such as a Planck function. If the spectra contain sine components, this indicates that the positions of the objects are shifted from their true positions. The reconstructed spectra are systematically affected by the erroneously estimated positions. Since the low-dispersion spectra of the atmospheres in exoplanets are close to a Planck function, the systematic shifts are not considered in this study. We note that the reflected light of Europa from the Sun in optical does not affect the position on the reconstructed image in spite of the spectrum including a number of lines (Matsuo et al. 2022).

Here, we provide an overview the concrete data reduction of this method. The process of the planet detection consists of the following five steps. First, we subtract the two chop states for each baseline to obtain the sine component of the complex visibility (step 1). Second, we extract the modulated signals of the off-axis point sources along the wavelength (step 2; see Sect. 2.2). Then, we repeat both of these steps during the rotation of the baseline for each set of the two-phase chop states (step 3). Next, we transform a set of reconstructed 1D images into 2D images (step 4; i.e., phase-space synthesis). Finally, we carry out the search for planetary light in the reconstructed 2D image and measure the planet position – if it indeed exists (step 5). We then proceed with the characterization process as follows:

We repeated step 1, but setting a longer integration time (step 6). Next, we model the wavelength dependence of the stellar leak from the collected data at short wavelengths (e.g., 4 to 6 µm) when the stellar leak is much brighter than the planet light due to a large OPD error (step 7). Then, we subtract the stellar leak from the collected data after extrapolating the stellar leak model constructed in the previous step to the long wavelength range (step 8). Next, we construct a matrix equation of the following form from the set of collected data: O = RI, where O is the observable vector, R is the response function, and I is the vector of the input sky (step 9). As the final step, we solve the matrix equation using the SVD method to extract the planetary spectrum (step 10; i.e., phase-space decomposition).

The phase information is summed for planet detection throughout the first five steps. In contrast, the phase information is decomposed in the planet characterization phase throughout the latter five steps. If the stellar leak is not much brighter than the planet light at short wavelengths, we can skip the modeling of the wavelength dependence of the stellar leak and its subtraction (i.e., steps 7 and 8). As shown in Sect. 3.2, when the systematic OPD rms error is 0.75 nm, corresponding to the standard requirement of LIFE (Dannert et al. 2022), the planetary spectra can be precisely extracted without those two steps. Although the first six steps are equivalent in terms of their data reduction, the required integration time is different. Because the planet light is integrated over the entire wavelength range in the reconstructed image through the first set of steps, the required integration time for planet detections is much shorter than that of the planetary spectrum obtained in the latter five steps (related to the characterization). We explain the planet detection and characterization in Sects. 2.2 and 2.3, respectively, and construct the equations required for PSSD.

2.2 Search for planetary signals (phase-space synthesis)

When a dual-Bracewell nulling interferometer with a $\frac{\pi }{2}$ phase chop observes the sky, the observed two-chop states in the unit of photoelectrons are as follows (e.g., Beichman & Velusamy 1999; Matsuo et al. 2011):

$$\begin{array}{l}{O}_{\pm }\left(\lambda \right)=\frac{1}{2}{\displaystyle \iint d^2\theta I\left(\lambda ,\theta \right){\sin }^2\left(\frac{\pi }{\lambda }b\cdot \theta +\delta l_n\right)}\hfill \\ \times \left\{1\pm \sin \left(\frac{2\pi }{\lambda }B\cdot \theta +\delta l_i\right)\right\},\hfill \end{array}$$(1)

where 6 is the position vector in the sky, b and B are the nulling and imaging baseline vectors, I(λ, θ) is the signal without the effect of sky transmission caused by the dual-Bracewell nulling interferometer at wavelength of λ, while δl_i and δl_n are the optical path differences of the imaging and nulling baselines, respectively. The (+) and (−) notations indicate the two chop states. In Eq. (1), we assumed that the two nulling baselines for the dual-Bracewell nulling interferometer have the same optical path difference error, δl_n, for simplicity.

The planetary system is the sum of the host star I_*,(λ, θ), multiple planets N_p, ${\displaystyle {\sum }_k^n{I}_{p,k}}$, local zodiacal light, $I_{lz}\left(\lambda \right)$, and the exozodiacal light I_ez(λ, θ). The spectrally resolved signal for the planetary system is expressed as:

$$\begin{array}{l}{O}_{\pm }\left(\lambda \right)=\frac{1}{2}{\displaystyle {\iint }_{{\text{Ω}}_{*}}{d}^2\theta I_{*}\left(\lambda ,\theta \right){\sin }^2\left(\frac{\pi }{\lambda }b\cdot \theta +\delta l_n\right)}\hfill \\ \times \left\{1\pm \sin \left(\frac{2\pi }{\lambda }B\cdot \theta +\delta l_i\right)\right\}\hfill \\ +\frac{1}{2}{\displaystyle \sum {}_k^{N_p}}{I}_{p,k}\left(\lambda ,{\theta }_{p,k}\right){\text{Ω}}_{p,k}{\sin }^2\left(\frac{\pi }{\lambda }b\cdot {\theta }_{p,k}+\delta l_n\right)\hfill \\ \times \left\{1\pm \sin \left(\frac{2\pi }{\lambda }B\cdot {\theta }_{p,k}\right)+\delta l_i\right\}\hfill \\ +\frac{1}{2}{\displaystyle {\iint }_{{\text{Ω}}_{fo\upsilon }}{d}^2\theta \left(I_{lz}\left(\lambda \right)+I_{ez}\left(\lambda ,\theta \right)\right){\sin }^2\left(\frac{\pi }{\lambda }b\cdot \theta +\delta l_n\right)}\hfill \\ \times \left\{1\pm \sin \left(\frac{2\pi }{\lambda }B\cdot \theta +\delta l_i\right)\right\},\hfill \end{array}$$(2)

where Ω_* and Ω_p,k are the solid angles of the host star and the kth planet, respectively, Ω_fov is the field of view of the interferometer, and θ_p,k is the position vector of the kth planet. As shown above (in step 1), the demolutated signal is given by:

$$\begin{array}{l}O\left(\lambda \right)=O_{+}\left(\lambda \right)-O_{-}\left(\lambda \right)\hfill \\ =\frac{1}{2}{\displaystyle {\iint }_{{\text{Ω}}_{*}}{d}^2\theta I_{*}\left(\lambda ,\theta \right){\sin }^2\left(\frac{\pi }{\lambda }b\cdot \theta +\delta l_n\right)}\hfill \\ \times \sin \left(\frac{2\pi }{\lambda }B\cdot \theta +\delta l_n\right)\hfill \\ +\frac{1}{2}{\displaystyle \sum {}_k^{N_p}{I}_{p,k}}\left(\lambda ,{\theta }_{p,k}\right){\Omega }_{p,k}{\sin }^2\left(\frac{\pi }{\lambda }b\cdot {\theta }_{p,k}+\delta l_n\right)\hfill \\ \times \sin \left(\frac{2\pi }{\lambda }B\cdot {\theta }_{p,k}+\delta l_i\right)\hfill \\ +\frac{1}{2}{\displaystyle {\iint }_{{\text{Ω}}_{fo\upsilon }}{d}^2\theta I_{ez}}\left(\lambda ,\theta \right){\sin }^2\left(\frac{\pi }{\lambda }b\cdot \theta +\delta l_n\right)\hfill \\ \times \sin \left(\frac{2\pi }{\lambda }B\cdot \theta +\delta l_i\right),\hfill \end{array}$$(3)

where the local zodiacal light was assumed to be removed from the demodulated signal because of its symmetrical structure. If the host star is perfectly positioned at the center of the field of view (FOV), the stellar leak disappears in Eq. (3) and contributes only as shot noise. Next, we move to step 2.

There are two approaches that can be applied in step 2: the cross-correlation method (Angel & Woolf 1997) or the Fourier transform (Matsuo et al. 2011). While the former focuses on the modulated signal while rotating the baseline, the latter uses the modulated one along the wavelength domain for each baseline. Here, we combine the advantages of the two methods. We used the correlation method (Angel & Woolf 1997) to extract the signal correlated to the modulation of the planet along the wavelength and derive the planet position for each baseline. After rotation of the baseline, the positions of the planets are obtained. We employed a rectangular array configuration with a baseline ratio of 6:1 based on the baseline of the LIFE mission concept (see Fig. 1). The parameters of the configuration are the same as those used for the numerical simulations in Sect. 3. The configuration is optimized to maximize the throughput of the habitable zone around a G-type star at 10 pc.

Given that the position vector of the correlated signal is α_corr (i.e., the 2D position of the signal), the positional information reconstructed from the jth imaging baseline vector, B_j, is

$$M_{\text{corr,}j}\left({\alpha }_j\right)={\text{Σ}}_i^{N_i}O\left({\lambda }_i\right)\sin \left(\frac{2\pi }{{\lambda }_i}{B}_j\cdot {\alpha }_{\text{corr}}\right){\sin }^2\left(\frac{\pi }{{\lambda }_i}{b}_j\cdot {\alpha }_{\text{corr}}\right),$$(4)

where λ_i is the ith spectral element, and N_i is the number of elements. Also, M_corr,j· tells us about the position of the correlated signal, α_j, projected to the baseline vector, B_j. As shown in Eq. (3), O(λ_i)is equal to the sum of the stellar leak, planetary signals, and background components. Each component has a different spectrum energy distribution.

In order to explain how PSSD works, we perform a simulation under a simple condition that an Earth-like planet is positioned at 1 AU parallel to the x-axis (panel a of Fig. 2). When the azimuths of the imaging baseline are 0 and 45° (panel b of Fig. 2), two 2D images parallel to the imaging baseline vector are generated (see panels c and d of Fig. 2). When the azimuth of the imaging baseline is 0°, the planet light is nulled for the rectangular array because the nulling baseline is parallel to the x-axis (panel c of Fig. 2). In fact, the peak value of the planet is almost 0. In contrast, for the azimuth angle of 45°, the planet light is extracted at 1 AU of the x-axis (panel d in Fig. 2). The peak value is a much larger than that of the nulled planet. The data obtained by the rectangular array also has information on the planet position in the direction perpendicular to the imaging baseline, thanks to the nulling baseline. The planet’s position along the nulling baseline is weakly constrained for each angle of the imaging baseline. Thus, both the imaging baseline and the nulling baseline can be utilized for planet detection.

Because the stellar leak caused by the OPD error is inversely proportional to approximately λ⁻⁴, the first term of O(λ_i) in Eq. (3) is not correlated with ($\frac{2\pi }{{\lambda }_i}{B}_j\cdot {\alpha }_{\text{corr}}$) sin² ($\frac{\pi }{{\lambda }_i}{b}_j\cdot {\alpha }_{\text{corr}}$) along the wavelength domain. The OPD error does not have less influence on the reconstruction of 1D images, compared to the previous method that extracts the modulated signal through the cross-correlation of the data collected while rotaing the baseline. In other words, PSSD has robustness against a large OPD fluctuation, which is the third advantage of PSSD, as discussed in Sect. 2.1.

The 2D positional information is obtained by summing the 1D images collected while rotating the baseline:

$$M_{\text{corr}}\left({\alpha }_{\text{corr}}\right)={\text{Σ}}_j^{N_j}{M}_{\text{corr,}j}\left({\alpha }_j\right),$$(5)

where N_j is the number of collected baselines. Panel (e) shows the 2D image reconstructed through Eq. (5). The pixel value at the planet position on the 2D image corresponds to the sum of the pixel values at the same position on the 1D images. Because PSSD focuses on the signal modulated by the wavelength for each baseline instead of one modulated by the rotation of the baseline, it is less affected by the long-term systematic noise than the purely rotation-based signal extraction, which points to the first advantage of PSSD discussed in Sect. 2.1. The required stability duration could be shortened to a few minutes, corresponding to the period for obtaining the two-phase chop states.

The PSSD method is not less impacted by a sparse U-V sampling. We discuss how the limited number of collected baselines affects the reconstructed 2D image in Sect. 4.2.

We also note the relationship between the correlation and Fourier tranform methods. The Fourier transformation of the demodulated signal along the wavelength gives a 1D image of the sky in parallel to the jth imaging baseline vector:

$$M_{\text{FT},j}\left({\alpha }_{\text{FT},j}\right)={\displaystyle \int d\left(\frac{1}{\lambda }\right)O\left(\lambda \right)\sin \left(\frac{2\pi }{\lambda }\left|B_j\right|{\alpha }_{\text{FT},j}\right){\sin }^2\left(\frac{\pi }{\lambda }\left|b_j\right|{\alpha }_{\text{FT},j}\right)},$$(6)

where α_FT is the 1D coordinate system in parallel to the jth imaging baseline vector. The origin of the coordinate system is the center of the field of view. Comparing Eq. (6) with Eq. (5), we found that both approaches are analytically equal if the sky consists of multiple point sources. A continuous source can be reconstructed only through Fourier transform of the interferometric signal (i.e., complex visibility). However, focusing on the fact that the transmission pattern of the sky induced by the nulling baseline can be utilized to increase the S/N of the planet’s detection, the Fourier transform method requires the condition that the imaging and nulling baselines are aligned for better planet detection. In contrast, the correlation method can be applied to any telescope configuration. Thus, the correlation method would be more utilized for planet detection, compared to the Fourier transform method (Matsuo et al. 2011).

Fig. 1 Rectangular array used for this study. Each filled circle represents a 2 m-diameter telescope. The ratio of the imaging to nulling baselines is 6:1.

Fig. 2 Procedure of image reconstruction. Panel a shows the planet position for this simulation. The planet is positioned at 1 AU in parallel to the x-axis, where (x, y) is the coordinate system of the object plane in the unit of AU. The planetary system at 10 pc consists of only one Earth-radius planet with an effective temperature of 285 K (without its host star). Panel b shows the directions of the imaging baseline for the two reconstructed 1D images (panels c and d). Panel e shows a 2D image converted from the set of the 1D images. The vertical axis of panels c and d and the color bar of panel e represent the number of photoelectrons. The parameters of the telescope and instrument are compiled in Table 2.

2.3 Extraction of planetary spectrum (phase-space decomposition)

Once the planetary light is successfully detected in the reconstructed 2D image, we can estimate the planetary spectrum. We calculate how the planet light is modulated while rotating the baseline based on the 2D planet positional information. The planetary signal for each spectral channel can be extracted from the data collected during baseline rotation. However, because the long-term fluctuation of OPD correlates with the modulation of the planetary signal during baseline rotation, the reconstruction of the planet spectrum is more easily affected by a long-term OPD error. In other words, the characterization of the planet light is more challenging than planet detection. Therefore, if the stellar leak is much brighter than the planetary signal, the bright stellar leak has to be subtracted before extracting the planetary signal from the data.

Here (and as introduced in Sect. 2.1), we consider the question of how the stellar leak changes while rotating the baseline could be measured from the data at short wavelengths. This is because the signals of warm and temperate planets except for a hot Jupiter are negligible compared to the stellar leak at short wavelengths. If we confirm from the reconstructed image that the stellar leak mainly contributes to the data at short wavelengths, the wavelength dependence of the stellar leak can be modeled in the same wavelength range. After the stellar leak model is extrapolated to the longer wavelength range, the stellar leak is subtracted from the demodulated signal (as shown in Eq. (3)).

There are two main systematic error terms: the first-order phase error and the cross-term of phase and amplitude, in the demodulated signal. Because the demodulated signal of the stellar leak in Eq. (3) is characterized by sin² ($\frac{\pi }{\lambda }b\cdot \theta +\delta l_n$) sin ($\frac{2\pi }{\lambda }B\cdot \theta +\delta l_i$), the wavelength dependence of the stellar leak could be simply expressed by a power law of wavelength under the condition that the systematic errors are much smaller than wavelength:

$$I_{\text{leak}}\left(\lambda ,t\right)=\left\{a{\left(\frac{\lambda }{{\lambda }_0}\right)}^{\alpha }+b{\left(\frac{\lambda }{{\lambda }_0}\right)}^{\beta }\right\}{I}_{\text{leak}}\left({\lambda }_0,t\right),$$(7)

where λ₀ is the reference wavelength, and I_leak(λ, t) is the stellar leak model in the demodulated signal at the wavelength λ as a function of time, t. However, we note that different wavelength dependencies may exist in the long wavelength range because a coating dispersion error or a pupil shear error could impact at long wavelengths. Because the coating dispersion and pupil shear errors drastically decreases in the longer wavelength regime, the subtraction of the estimated stellar leak model from the demodulated signal would less impact the reconstructed planetary spectrum at long wavelengths. We also note that thanks to the bright stellar leak at short wavelengths, the systematic aberration could be modeled from the intensity and wavelength dependence of the stellar leak, which is a similar work as measuring coronagraphic low-order aberrations (Guyon et al. 2009).

After subtracting the stellar leak from the demodulated signal shown in Eq. (3), the planetary spectra are reconstructed from the residual data through the SVD method. The matrix equation for the ith spectral element can be written as:

$$O_i=R_i{I}_i,$$(8)

where O is the vector composed of the observed data, R is the matrix of the response function of objects (i.e., sky transmission of object), and I is the vector of the input sky. In order to reconstruct the spectra of the planets, we solve the matrix equation for each spectral channel.

After the observed data is subtracted from the averaged value over the baseline rotation, the O vector of the ith spectral element is

$$O_i=\left(\begin{array}{c}{O}_{i,1}\\ O_{i,2}\\ ⋮\\ O_{i,N_j}\end{array}\right),$$(9)

where O_{i, j} is the observation data of the jth azimuth angle for the ith spectral element. The number of the elements for the O vector is N_j, corresponding to the number of the collected data during baseline rotation. The response matrix for the ith spectral element, R_i, is written as:

$$R_i=\left(\begin{array}{ccc}{R}_{1,1}& \cdots & R_{1,N_p}\\ ⋮& & ⋮\\ R_{N_j,1}& \cdots & R_{N_j,N_p}\end{array}\right),$$(10)

where we assume that the continuum component is removed from the observed vector, O_i, by the subtraction of the two chop states. The R matrix is a N_j × N_p matrix. Each component is expressed as follows:

$$\begin{array}{l}{R}_{1,1}={\sin }^2\left(\frac{\pi }{{\lambda }_i}{b}_1\cdot {\theta }_{p,1}\right)\sin \left(\frac{2\pi }{{\lambda }_i}{B}_1\cdot {\theta }_{p,1}\right),\hfill \\ R_{1,N_p}={\sin }^2\left(\frac{\pi }{{\lambda }_i}{b}_1\cdot {\theta }_{p,N_p}\right)\sin \left(\frac{2\pi }{{\lambda }_i}{B}_1\cdot {\theta }_{p,N_p}\right),\hfill \\ R_{N_j,1}={\sin }^2\left(\frac{\pi }{{\lambda }_i}{b}_{N_j}\cdot {\theta }_{p,1}\right)\sin \left(\frac{2\pi }{{\lambda }_i}{B}_{N_i}\cdot {\theta }_{p,1}\right).\hfill \end{array}$$(11)

The vector of the input sky for the ith spectral element is

$$I_i=\left(\begin{array}{c}{I}_{p,1}\left({\lambda }_i,{\theta }_{p,1}\right)\\ I_{p,2}\left({\lambda }_i,{\theta }_{p,2}\right)\\ ⋮\\ I_{p,N_p}\left({\lambda }_i,{\theta }_{p,N_p}\right)\end{array}\right).$$(12)

The number of elements for the I vector is N_p. When the number of the observed data is much larger than that of the elements for the input matrix, the N_p planetary signals for each spectral element can be decomposed via the SVD method.

Finally, we need to emphasize that there are several ways to decompose the planetary signals and unknown stellar leaks under the condition that the planets are successfully detected. This decomposition could be solved using modified orthogonal projections, or kernels, such as those used in (Laugier et al. 2020). They preserve important properties of the covariance matrix of errors, and this decomposition therefore is well suited for further data whitening approaches.

We could also combine PSSD with wavelet-based signal reconstruction methods (e.g., del Ser et al. 2018), either by suppressing wavelet coefficients at different levels, corresponding to unwanted signals, or convolving them with specially design convolution kernels. It is expected that while the low-frequency signal is mainly caused by the systematic OPD residual, the mid-to high-frequency signals are caused by off-axis point sources and stochastic noises. Knowing the models of the systematic errors and stochastic noises makes it possible to suppress sections of the wavelet decomposition connected to these signals and then, using the inverse wavelet transform, we can rebuild the planetary signal. Reconstruction of the signal will depend on the cadence and its level with respect to noise. With advanced techniques, we could expect to see signal recovery, even when its contribution is up to ~10−30%.

3 Simulations

We performed numerical simulations to check the feasibility of PSSD under the LIFE baseline scenario. First, we briefly explain the simulation setup regarding the target system and instrument. Next, we show the results generated by PSSD under the ideal condition where only the astronomical noise contributes to the data as shot noise. Finally, we include a long-term systematic OPD error in the simulations and show its impact on PSSD.

3.1 Setup

The distance of the considered target is 10 pc. The target system consists of three Earth-sized (R_⊕) planets, a Sun-like star with a solar radius (R_⊙), and an effective temperature of 5778 K, and an exozodiacal dust disk. The semi-major axes of the three planets are 1, 0.73, and 1.5 AU, that is, the same as those of Earth, Venus, and Mars. Given that the effective temperatures of the three planets are simply proportional to the inverse square root of the semi-major axis, the temperatures of planets P1, P2, and P3 were set to 285, 330, and 232 K, respectively. All of the target objects were assumed to emit blackbody radiation. The orbital phases of the three planets were set to 0, −45, and 90°. The phase angle of 0° points along the positive x-axis, where (x, y) is the coordinate system of the sky in the unit of AU. The arrangement of the three planets is shown in Fig. 3a. The exozodiacal light is equal to three times that of the solar system (Ertel et al. 2020). The surface brightness of the exozodiacal light is generated based on the previous model (Kennedy et al. 2015), which is applied to the software tool, LIFEsim (Dannert et al. 2022). Table 1 compiles all the parameters of the target system.

We employed a dual-Bracewell nulling interferometer with imaging and nulling baselines of 87.3 m and 14.55 m (see Fig. 1) so that the maximum of transmission is achieved for the center of the habitable zone around a Sun-like star at 10 pc at a wavelength of 15 µm (Quanz et al. 2022). The diameter of each telescope is 2 m, and the imaging and nulling baselines are perpendicular to each other. Although the observing wavelength ranges from 4 to 18.5 µm, the same as the LIFE baseline, we limited the wavelength range to larger than 8 µm in the planet detection phase. The reason is that the bright stellar leak is more than 100 times brighter than the light of the planets we consider at short wavelengths and deteriorates the performance of PSSD. We note that the shorter wavelengths are effective in looking for inner planets because of the combination of higher spatial resolution and brighter planets in that wavelength range. Thus, PSSD must optimize the wavelength range used for planet detection based on what type of planets we find.

We set the resolving power of the spectrum to 50 for both planet detection and their characterization. The minimum resolving power for planet detection is determined by the required field of view. When the resolving power is 50, the field of view is 1.14 arcsec at 10 µm. For the characterization, a spectral resolution of 30–50 was suggested in Konrad et al. (2022) to detect the various molecules in an Earth-twin atmosphere. The total throughput was set to 0.035, given that the instrument throughput and quantum efficiency are 0.05 and 0.7, respectively. The integration time was set to 55 h for planet detection and 75 days for planet characterization, respectively.

We do not assume a continuous rotation but a discrete rotation in steps of one degree, where the spacecraft come to a halt before rotating again by one degree because of computational cost. We note that the continuous rotation provides a better reconstruction thanks to the continuous U-V coverage compared to the discrete rotation. Assuming that the baseline rotates by 360° at a one-degree interval, the integration time for each baseline is 550 and 180 000 seconds for planet detection and characterization phases, respectively. In addition to the ideal observing case, we studied the feasibility of PSSD under a systematic OPD error. Although there are mainly the first-order phase and the phase-amplitude cross-term in the demodulated signal (Lay 2004), only the former was considered in this simulation. We note, however, that this simulation on the feasibility of PSSD is not largely impacted by the phase-amplitude cross-term because the two systematic components have a similar frequency dependence of $\frac{1}{f}$, called “pink noise,” where f is the frequency. The root mean square (rms) of the OPD error was set to 0.75 nm, corresponding to the LIFE baseline scenario for the case of only phase error (Dannert et al. 2022). The baseline value is larger than that for the case of both phase and amplitude errors. We also consider 5, 10, and 15 times the baseline values, 3.8, 7.5, and 11.3 nm rms errors, to investigate the limitation of PSSD in Sect. 4. Table 2 compiles all the instrumental parameters.

Figure 3b shows the astronomical signals obtained by the Bracewell nulling interferometer under the above-observed conditions. The stellar leak and background, such as the local zodiacal and exozodiacal light, cover the modulations of the three planets. When there is no OPD error, all astronomical signals contribute to the data as the shot noise. We perform the numerical simulations under the ideal condition in Sect. 3.2.1 and then consider the fluctuation of the stellar leak due to the systematic OPD error in Sect. 3.2.2.

Fig. 3 Target planetary system. Left: Configuration of three planets. Right: Signals of planet P1 (black), planet P2 (blue), planet P3 (green), the nulled host star (red), local zodiacal light (brown), and exozodiacal light (purple) with a resolving power of 50 per a unit of time. Target system and instrument parameters are compiled in Tables 1 and 2.

3.2 Results

3.2.1 Ideal condition

We first show the feasibility of PSSD under the ideal condition, in which only shot noise exists due to astronomical sources. After collecting data while rotating the baseline by 360° in steps of 1 degree, we generated a 2D image through steps one through four. Figure 4a shows the reconstructed 2D images for an integration time of 55 h. The S/N for the detection of planets under the ideal condition are compiled in case 1 of Table 3. We successfully detected signals of planets P1 and P2 with S/N of 10.8 and 14.6, respectively. The higher temperature of planet P2 allows us to obtain higher S/N compared to that of planet P1. In contrast, the S/N of planet P3 is only 3.6 because of its lower temperature. We require a longer integration time to achieve the S/N of five for planet P3.

The S/N is defined as the ratio of the planetary signal to the starndard deviation at the same angular distance as its planet. In order to calculate the standard deviation, a 2D image without the planetary signals is generated and divided into annular rings. The noise floor is calculated as the standard deviation for each annular ring. We note that the absolute value of the S/N cannot be directly compared with that calculated by Dannert et al. (2022) because PSSD reconstructs both the signal and the noise in a different way.

Next, we reconstructed spectra of the three planets through steps of 6, 9, and 10, assuming that the planet positions are correctly obtained. Figures 4b, c, and d show the reconstructed spectra of the three planets. The spectra of planet P1 and planet P2 are consistent with the input spectra (solid gray lines). We also derived for each data-point the average and standard deviation by performing the numerical simulations 100 times. The S/N of the spectrum for planet P1 agrees with that of the previous study (Konrad et al. 2022). LIFE could detect the methane and ozone absorption bands at 7.6 and 9.6 µm with S/N of approximately 5 and 15, respectively. This combinationofsimultaneously detected absorption features is thought to be a good indicator of a non-equilibrium atmosphere caused by biological activity on the planet (Kasting et al. 2014).

On the other hand, the S/N worsens at shorter wavelengths than 7.5 µm. This is because the stellar leak drastically increases due to the narrower null pattern on the sky in the shorter wavelength range. The planetary signals rapidly decrease at the same time (see Fig. 3b). In addition, PSSD also obtained the entire spectrum of planet P3 but at a S/N lower than three except for wavelengths longer than 10 µm due to its faintness compared to the other planets. Table 4 compiles the S/N of the reconstructed spectra at wavelengths of5, 7.5, 10, and 15 µm.

Although planet P3 was not detected for an integration time of 55 h, the signal of planet P3 could be obtained while integrating the data in the characterization phase, and its position would be well determined. Planets P1 and P2 were also detected with higher S/N compared to those in the planet detection phase, which can reduce the systematic errors of the reconstructed spectra due to the estimation errors of the planet positions. Thus, we confirmed that PSSD could detect the planetary signals and characterize their atmospheres under the ideal condition, in which the data is affected only by the shot noise due to the stellar leak, background, and planets.

Fig. 4 Image and spectrum reconstruction under photon-noise limited condition. (a) Reconstructed 2D image. The white arrows denote positions of planet P1, planet P2, planet P3, respectively. The integration time was set to 55 h. The unit of the color bar is the number of photo-electrons. Reconstructed spectra of three planets (b) P1, (c) P2, and (d) P3 for an integration time of 75 days. The grey line and grey vertical bar of each panel show the input model and the standard deviation of each data point derived through 100 numerical simulations, respectively.

3.2.2 Systematic error

We consider instrumental noise and investigate its negative impact on planet detection and characterization. In order to evaluate it, we included an OPD error in the numerical simulations as the instrumental noise. However, we note that the phase-amplitude cross-term also contributes to the long-term stellar fluctuation in the demodulated signal under the existence of the amplitude error (Lay 2004). Because both systematic components have the same spectrum in terms of the time domain (e.g., Dannert et al. 2022), the phase-amplitude cross-term would not significantly impact the results. We assumed that the systematic noises have a dependency of $\frac{1}{f}$, where f shows the frequency.

According to Dannert et al. (2022), when the OPD rms error is larger than 0.75 nm, the instrumental noise is dominant over the statistical noise (i.e., fundamental noise) from the astronomical objects at the shortest wavelength. Because the systematic OPD error is much smaller than the observing wavelength, the amount of stellar leak is proportional to the OPD error for each spectral element. The OPD error impacts the null depth and leaves the stellar leak in the subtraction of the two chop states (Eq. (3)). While the former contributes to the data as the Poisson noise, the latter affects the planetary signal during baseline rotation, which prevents us from reconstructing the planetary spectrum. In our simulations, we added the same OPD error to the imaging and nulling baselines to reduce the calculation cost.

Figure 5a compares the spectra of the three planets with the stellar leak left after subtracting the two chop states in the entire wavelength range (i.e., 4–18.5 µm). The stellar leak drastically increases at the shorter wavelengths because of the shallower null depth and the planetary signal drops instead. It is more challenging to perform planet detection and characterization at shorter wavelengths than longer ones. Panels b, c, and d of Fig. 5 compare the modulated signal of planet P1 with the stellar leak at 4, 8, and 12 µm as a function of the azimuth of the imaging baseline. Although the systematic OPD rms error of 0.75 nm does not affect the planetary signal at 12 µm, the stellar leak covers the planetary signal at 4 µm, which is consistent with Fig. 5a.

Figure 6a shows a reconstructed 2D image under the condition that the systematic OPD rms error is 0.75 nm. Thanks to robustness of PSSD against a long-term OPD error, we obtained the same signal-to-noise ratios for detecting the three planets as those for the ideal state.

Finally, we reconstructed the planetary spectra over the entire wavelength based on steps 6, 9, and 10 without subtraction of the stellar leak from the demodulated signal in steps 7 and 8. Panels b, c, and d of Fig. 6 compare the reconstructed spectra of the three planets with the input models. The reconstructed spectra are consistent with the models over the entire wavelength range, except for the shorter range. Comparing the reconstructed spectra with the ideal case, we found that the shot noise limits the performance of PSSD in the wavelength range longer than 7.5 µm. In addition, even though the systematic stellar leak is a few times brighter than the planetary signals (Fig. 5a), the S/N are almost the same as those for the reconstructed spectra under the ideal condition (case 2 of Table 4). This is because the SVD method efficiently separates the mid- to high-frequency components induced by planets from the low-frequency component due to the systematic stellar leak. Thus, the SVD method could reconstruct the spectra of the three planets from the modulations of the planetary signals during baseline rotation while efficiently separating them from the long-term systematic stellar leak.

Fig. 5 Systematic stellar noise. Panel a shows comparison of signals of the three planets P1 (black), P2 (blue), and P3 (green) with the nulled stellar leaks left in the subtraction of the two chop states over the entire wavelength range. The red, brown, gray, and light gray lines represent the stellar leaks for systematic OPD rms errors of 0.75, 3.8, 7.5, and 11.3 nm. Panels b, c, and d show fluctuations in the stellar leak due to an OPD rms error of 0.75 nm (red) and the demodulated signal of planet P1 (black) at wavelengths of 4, 8, and 12 µm, respectively. The integration time of each data point is 550 s.

Fig. 6 Same as Fig. 4 except for a systematic OPD rms error of 0.75 nm.

Fig. 7 Image reconstruction under large systematic noises. Panels a, b, and c show reconstructed images for systematic OPD rms errors of 3.8, 7.5, and 11.3 nm, respectively. Panel d shows the noise floors for systematic OPD errors of 3.8 (black), 7.5 (yellow), and 11.3 nm (brown) were compared with the planetary signals (star symbol). The noise floor for each OPD error was calculated for a reconstructed 2D image without the planetary signals. The gray line shows the noise floor for the ideal case (i.e., only the shot noise) as a reference. Because the reconstructed planetary signals are slightly affected by the OPD error, the black, yellow, and brown star symbols represent the planetary signals for OPD rms errors of 3.8, 7.5, and 11.3 nm, respectively.

4 Discussion

Thus far, we have confirmed that PSSD detects planet light and extracts planetary spectra under the existence of OPD fluctuation, which follows $\frac{1}{f}$. In this section we investigate how much noise amplitude can be withstood by PSSD (Sect. 4.1) and we compare PSSD with the previous method (Sect. 4.2).

4.1 Robustness against a large OPD error

We set the systematic OPD rms error to 3.8, 7.5, and 11.3 nm, values that are equal to the baseline values multiplied by factors of 5, 10, and 15, respectively. Figures 7b, c, and d show the reconstructed images for these three different OPD errors. As the OPD error increases, systematic patterns are brighter in the centre region of each image. In contrast, the noise floor is limited by the shot noise at semi-major axes larger than 1.0 AU (panel Fig. 7d). As discussed in Sect. 3.2.1, the noise floor was calculated for a reconstructed 2D image without the planetary signals.

We derived the S/N of the three planets for each OPD rms error. Although the S/N is gradually worsened as the OPD error increases, the S/N of the planets are higher than 5 except for planet 3 (cases 2–4 of Table 3). Thus, PSSD could successfully detect the inner two planets even under systematic OPD errors 15 times larger than the LIFE baseline requirement. In addition, because the shot noise limits the detection of planet 3, PSSD could also detect planet 3 with a longer integration time.

Next, we reconstructed the spectra of the three planets for an OPD rms error of 7.5 nm by solving the matrix equation with the SVD method (step 9; see left panels of Fig. 8). However, the OPD rms error of 7.5 nm deforms the spectra of planet P1 and planet P2 at the shorter wavelengths than 10 µm and the spectrum of planet P2 over the entire wavelength range. This is because the planetary spectra are reconstructed from the long-term modulation of the signal while rotating the data, which correlates with the systematic OPD error. Thus, the characterization of the planetary atmosphere is much more affected by the long-term OPD error than the planet detection.

Here, we utilize the advantages of PSSD in terms of planet detection. Because a 2D image can be reconstructed even for a large OPD error, we are able to investigate what kind of objects orbits the host star. Once we confirm that the stellar leak is dominant over the planetary signals at short wavelengths, we can then measure the long-term fluctuation of the stellar leak from the data at short wavelengths (step 7) and subtract the stellar leak from the demodulated signal (step 8). There are several steps required to reconstruct the planetary spectrum. First, we apply a low pass filter to the demodulated signal to decrease the impact of statistical noise on the data. Second, the wavelength dependence of the stellar leak is estimated from a large number of the data points collected during baseline rotation. Finally, based on the estimated wavelength dependence, we extrapolate the stellar leak model to the longer wavelength range and subtract it from the demodulated signal shown in Eq. (3) over the entire wavelength. The planetary spectra are reconstructed through applying the subtracted data to the SVD process. Panels b, d, and f of Fig. 8 show the reconstructed spectra of the three planets through the above process for a large OPD rms error of 7.5 nm, corresponding to ten times larger than the baseline requirement of LIFE. In this simulation, the wavelength range applied to estimation of the wavelength dependence was 4–5.5 µm, in which the stellar leak is more than 100 times greater than planetary signals for an OPD rms error of 7.5 nm (Fig. 5a). The S/N for the reconstructed spectra (case 4 of Table 4) are almost the same as those under the ideal condition (case 1 of Table 4). Thus, if we confirm, based on the reconstructed image, that only the stellar leak mainly contributes to the demodulated signal at short wavelengths, the planetary spectra could be reconstructed.

Fig. 8 Spectrum reconstruction under large systematic noise. Panels a, c, and e show the reconstructed spectra of planet P1, planet P2, and planet P3 only through the SVD process. Panels b, d, and ſ show the spectra of planet P1, planet P2, and planet P3 reconstructed through subtracting the stellar leak from the demodulated signal before the SVD method. The OPD rms error was set to 7.5 nm for both the two cases.

Fig. 9 Same as Fig. 7 except for using the cross-correlation method. The noise floor derived under the LIFE baseline requirement (dashed line) is added as a reference.

4.2 Comparison with previous method

We compare PSSD with the previous signal extraction through the cross-correlation or the maximum likelihood of the data obtained while rotating the baseline (e.g., Angel & Woolf 1997; Dannert et al. 2022). The main difference between the two methods is whether a ID image from each baseline is first reconstructed or a 2D image is reconstructed at one time. PSSD performs the correlation of planetary signal among the obtained spectrum and converts the ID information into a 2D image. In contrast, the previous method simultaneously finds the modulation of the planetary signal in both the wavelength- and time-domains. As discussed in Sect. 2.1, the advantages of PSSD are robustness against a large OPD error and a limited number of baselines.

Panels a, b, and c of Fig. 9 show the reconstructed 2D images through the previous cross-correlation method. The images were formed based on the following equation:

$$M_{\text{corr}}\left({\alpha }_{\text{corr}}\right)={\displaystyle {\sum }_j^{N_j}{\displaystyle {\sum }_i^{N_i}O\left({\lambda }_i\right)\sin \left(\frac{2\pi }{{\lambda }_i}{B}_j\cdot {\alpha }_{\text{corr}}\right){\sin }^2\left(\frac{\pi }{{\lambda }_i}{b}_j\cdot {\alpha }_{\text{corr}}\right)}}.$$(13)

Compared with Fig. 7, the OPD error induces brighter systematic patterns at semi-major axes smaller than 1.0 AU, which prevent us from detecting the two inner planets. As shown in panel d of Fig. 9, the intensity of the systematic pattern is roughly proportional to the OPD rms error in the inner region. We note that the noise floor under the baseline requirement of LIFE is almost equal to that for a systematic OPD rms error of 0.75 nm, which is consistent to the previous study (Dannert et al. 2022). The S/N for the detection of planets with PSSD under an OPD rms error of 7.5 nm (case 2 of Table 3) are almost the same as those for the cross-correlation method under the ideal case (case 5 of Table 3). Therefore, there exists a large difference between the robustness against a large OPD error.

We also compare PSSD with the cross-correlation method in terms of the impact of a limited number of baselines on planet detection. The left panels of Fig. 10 show the PSSD reconstructed images with a limited number of baselines and under the existence of only long-term systematic error (without shot noise). In other words, the reconstructed image is not influenced by the integration time. We randomly selected available baselines, which are more sparsely distributed over 360° as the number of baselines decreases. The systematic pattern is slightly brighter as the number of baselines decreases. However, the three planets could be detected even though the fraction of the available baselines is limited to only 8%. In contrast, images reconstructed using the previous method were largely affected by the limited number of baselines because the modulation in the time domain was lost (the right panels of Fig. 10). In addition, the artificial pattern fully covers the three planets if the number of baselines is limited to 8%, and the limited azimuth coverage elongates the point sources, which means that the long-term systematic error modulates the data collected while rotating the baseline.

Thus, PSSD is much less impacted by both a large systematic error and a limited number of available baselines compared to the previous method, which could relax the requirement of LIFE.

Fig. 10 Images reconstructed through PSSD (left panels) and the previous method (right panels) under the condition that the number of baselines is limited by a, b 100%, c, d 28%, and e, f 8%. Only a long-term systematic error was included in the simulations; shot noise was not considered while investigating how the limited number of the baselines affects the image reconstruction. The OPD rms error was set to 3.8 nm.

5 Conclusion

We proposed a method for planet detection and characterization with future nulling space interferometers, such as large interferometer for exoplanets (LIFE). The proposed method is named “phase-space synthesis decomposition” (PSSD). It focuses on the correlation of the planetary signal over the entire wavelength range instead of that along the baseline rotation. Because a 1D image parallel to the baseline can be derived for each baseline, a large number of 1D images are collected after rotating the baseline. A 2D image can be reconstructed by summing over the 1D images. Once the 2D image is obtained, a continuous equation is constructed based on the planet position information, and its solution through singular value decomposition (SVD) allows us to extract the planetary spectra embedded in the stellar fluctuation. As long as the modulation of the planetary signal has a different frequency from the stellar fluctuation, the SVD method efficiently decomposes the stellar leak and planetary signal. The PSSD method provides three advantages in planet detection compared to previous methods that find a modulation of the planetary signal during baseline rotation in both the wavelength and the time domains. One is robustness against a large systematic OPD error. Because the stellar leak has a different wavelength dependence from the planetary signal, only the correlation of the signal efficiently decomposes the stellar leak and the planetary signal. The second is that PSSD does not correlate with a long-term fluctuation of the stellar leak because a 2D image is formed by summing over 1D images that require only two-phase chop states. The third advantage is robustness against a limited number of baselines.

We performed numerical simulations to investigate the feasibility of PSSD under various conditions. We put three terrestrial planets with semi-major axes of 0.73, 1, and 1.5 AU, corresponding to those of Venus, Earth, and Mars, respectively, around a Sun-like star at 10 pc. The simulation included both statistical and systematic noises. PSSD successfully detected the three planets and reconstructed their spectra for an OPD rms error of 0.75 nm, which is the same as the baseline requirement of LIFE.

We also confirmed that PSSD has robustness against a large systematic OPD error. We increased the amplitude of the systematic OPD error by a factor of 5, 10, and 15. PSSD successfully detected the three planets almost without being affected by the stellar leak, even under the largest systematic OPD error of 11.3 nm. This is because the stellar leak is inversely proportional to approximately the fourth power of wavelength, which is largely different from the modulation of the planetary signal in the wavelength domain.

In contrast, the reconstructed spectra were more affected by the long-term systematic noise than planet detection. This is because PSSD uses the modulation of the planetary signal during baseline rotation to reconstruct the spectra of the three planets. The long-term systematic OPD error more easily correlates with the modulation of the planetary signal. The spectra could not be accurately extracted under an OPD rms error of 7.5 nm, corresponding to ten times the baseline requirement of LIFE. The signal-to-noise ratio significantly decreases in the wavelength range shorter than 7.5 µm.

Here, focusing on the fact that the planetary signals are much smaller than the stellar leak at shorter wavelengths than 6 µm, we can measure the fluctuation of the stellar leak and its wavelength dependence for the data at short wavelengths. Because PSSD can successfully obtain the planetary signals even under a large systematic noise, PSSD utilizes the planet position information and tells us what type of sources contribute to the signal at short wavelengths. If we confirm from the reconstructed 2D image that the stellar leak is dominant over the planetary signals at short wavelengths, the wavelength dependence of the stellar leak can be modeled from the data in the short wavelength range. After extrapolating the estimated stellar leak model to the longer wavelength range, the stellar leak could be subtracted from the data over the entire wavelength range. The spectra of the three planets were successfully reconstructed by applying the subtracted data to the SVD process. The S/N were almost the same as those for the ideal condition.

Finally, we compared PSSD with a previous method that reconstructs a 2D image by simultaneously fitting the modulation of the planet in the time and the wavelength domains after the baseline is rotated. Because the long-term noise is correlated with the planetary signal in the time domain, systematic patterns were formed in the reconstructed image and covered the planetary signals under systematic OPD errors larger than 3.8 nm. The S/N for a planet detection significantly decreased for large OPD errors compared to PSSD. This is because the noise floor increased in the inner region due to the systematic OPD error. In addition, limited azimuth coverage in the U-V plane impacted planet detection because the modulation in the time domain was lost. In contrast, PSSD can reconstruct a 2D image from fewer baselines. Even in the case where the azimuth coverage of the baseline is limited to 8%, the three planets could be discovered by PSSD.

Thus, PSSD is more robust against a large OPD error and a limited number of baselines, which could relax the requirements of LIFE regarding the OPD error and the stability duration. However, this numerical simulation was performed as the first step under an ideal case that makes detection of terrestrial planets easier. As a next step, we will investigate various impacts not considered in this study, such as asymmetric exozodiacal structure (Defrère et al. 2010), and the other systematic errors, on planet detection and characterization.

Acknowledgements

We express our sincere gratitude to Dr. Lacour for many valuable comments and suggestions on this study. Part of this work was supported by JST FOREST Program, Grant Number JPMJFR202W. Part of this work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. F.D. and S.P.Q. acknowledge the financial support of the SNSF. Part of this work was carried out within the project SCIFY, which has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement CoG – 866070). A.B.K. acknowledges funding provided by University of Belgrade-Faculty of Mathematics (the contract 451-03-47/2023-01/200104), through the grants by the Ministry of Science, Technological Development and Innovation of the Republic of Serbia.

References

All Tables

All Figures

	Fig. 1 Rectangular array used for this study. Each filled circle represents a 2 m-diameter telescope. The ratio of the imaging to nulling baselines is 6:1.
In the text

Fig. 2 Procedure of image reconstruction. Panel a shows the planet position for this simulation. The planet is positioned at 1 AU in parallel to the x-axis, where (x, y) is the coordinate system of the object plane in the unit of AU. The planetary system at 10 pc consists of only one Earth-radius planet with an effective temperature of 285 K (without its host star). Panel b shows the directions of the imaging baseline for the two reconstructed 1D images (panels c and d). Panel e shows a 2D image converted from the set of the 1D images. The vertical axis of panels c and d and the color bar of panel e represent the number of photoelectrons. The parameters of the telescope and instrument are compiled in Table 2.

In the text

	Fig. 3 Target planetary system. Left: Configuration of three planets. Right: Signals of planet P1 (black), planet P2 (blue), planet P3 (green), the nulled host star (red), local zodiacal light (brown), and exozodiacal light (purple) with a resolving power of 50 per a unit of time. Target system and instrument parameters are compiled in Tables 1 and 2.
In the text

Fig. 4 Image and spectrum reconstruction under photon-noise limited condition. (a) Reconstructed 2D image. The white arrows denote positions of planet P1, planet P2, planet P3, respectively. The integration time was set to 55 h. The unit of the color bar is the number of photo-electrons. Reconstructed spectra of three planets (b) P1, (c) P2, and (d) P3 for an integration time of 75 days. The grey line and grey vertical bar of each panel show the input model and the standard deviation of each data point derived through 100 numerical simulations, respectively.

In the text

Fig. 5 Systematic stellar noise. Panel a shows comparison of signals of the three planets P1 (black), P2 (blue), and P3 (green) with the nulled stellar leaks left in the subtraction of the two chop states over the entire wavelength range. The red, brown, gray, and light gray lines represent the stellar leaks for systematic OPD rms errors of 0.75, 3.8, 7.5, and 11.3 nm. Panels b, c, and d show fluctuations in the stellar leak due to an OPD rms error of 0.75 nm (red) and the demodulated signal of planet P1 (black) at wavelengths of 4, 8, and 12 µm, respectively. The integration time of each data point is 550 s.

In the text

	Fig. 6 Same as Fig. 4 except for a systematic OPD rms error of 0.75 nm.
In the text

Fig. 7 Image reconstruction under large systematic noises. Panels a, b, and c show reconstructed images for systematic OPD rms errors of 3.8, 7.5, and 11.3 nm, respectively. Panel d shows the noise floors for systematic OPD errors of 3.8 (black), 7.5 (yellow), and 11.3 nm (brown) were compared with the planetary signals (star symbol). The noise floor for each OPD error was calculated for a reconstructed 2D image without the planetary signals. The gray line shows the noise floor for the ideal case (i.e., only the shot noise) as a reference. Because the reconstructed planetary signals are slightly affected by the OPD error, the black, yellow, and brown star symbols represent the planetary signals for OPD rms errors of 3.8, 7.5, and 11.3 nm, respectively.

In the text

Fig. 8 Spectrum reconstruction under large systematic noise. Panels a, c, and e show the reconstructed spectra of planet P1, planet P2, and planet P3 only through the SVD process. Panels b, d, and ſ show the spectra of planet P1, planet P2, and planet P3 reconstructed through subtracting the stellar leak from the demodulated signal before the SVD method. The OPD rms error was set to 7.5 nm for both the two cases.

In the text

	Fig. 9 Same as Fig. 7 except for using the cross-correlation method. The noise floor derived under the LIFE baseline requirement (dashed line) is added as a reference.
In the text

Fig. 10 Images reconstructed through PSSD (left panels) and the previous method (right panels) under the condition that the number of baselines is limited by a, b 100%, c, d 28%, and e, f 8%. Only a long-term systematic error was included in the simulations; shot noise was not considered while investigating how the limited number of the baselines affects the image reconstruction. The OPD rms error was set to 3.8 nm.

In the text