The coherent differential imaging on speckle area nulling (CDI-SAN) method for high-contrast imaging under speckle variation

Differential imaging is a postprocessing method to obtain high contrast, often used for exoplanet searches. The coherent differential imaging on speckle area nulling (CDI-SAN) method was developed to detect a faint exoplanet lying beneath residual speckles of a host star. It utilizes image acquisitions faster than the stellar speckle variation synchronized with five shapes of a deformable mirror repeatedly. By using the only the integrated values of each of the five images and square differences for a long interval of observations, the light of the exoplanet could be separated from the stellar light. The achievable contrast would reach to almost the photon-noise limit of the residual speckle intensities under appropriate conditions. The CDI-SAN can be applied to both ground-based and space telescopes.


INTRODUCTION
Various high-contrast imaging methods have been developed to detect the light of exoplanets for characterization. Planet intensities observed from reflected light have been estimated to be about 10 −8 to ∼ 10 −10 of the intensity of its host star. To achieve such high-contrast, diffracted starlight in the telescope should be removed by a coronagraph, and residual starlight should be suppressed by wave-front control. Observed images would be often reduced by postprocessing techniques called differential imaging, which would use two or more images at different wavelengths, polarizations, or orientation angles to further suppress the residual starlight speckles.
Real-time focal-plane wave-front sensing and control methods utilized speckle intensities with some coherent-light modulations and nulled the speckles to produce high-contrast regions mainly in laboratory experiments leading to space telescopes (e.g., Guyon 2004;Trauger et al. 2004;Bordé & Traub 2006;Baudoz et al. 2006;Give'on et al. 2011;Oya et al. 2015). In this process, speckle intensities were derived separately from the incoherent component (planet light). Coherent differential imaging (CDI) is a method that uses the same principle in postprocessing (e.g., Jovanovic et al. 2018). In a ground-based telescope using wave-front control by a wave-front sensor, the coherent-light technique was applied to reduce quasi-static diffracted light and speckles (e.g., Martinache et al. 2014;Bottom et al. 2017).
Recently, a fast operation in the CDI method mentioning the potential of long exposure using a self-coherent camera was proposed (Gerard et al. 2018). In the present paper, a fast synchronized detection process and long-exposure technique based on the electric field conjugation (EFC; Give'on et al. 2011) and the speckle area nulling (SAN; Oya et al. 2015) algorithms is presented; this is a kind of CDI, and we call it the CDI on speckle area nulling (CDI-SAN) method. The principle behind CDI-SAN will be described in section 2. Numerical simulations and discussions will be shown in section 3. Conclusions will be described in section 4.

PRINCIPLE BEHIND CDI-SAN
The EF at the focal plane can be modulated by controlling a deformable mirror (DM) at the pupil plane. Let the intensity without modulation at a detector pixel in the final focal plane be then those with some modulations can be written as where I s and E s are the intensity and an EF of the speckle with zero modulation, respectively, i.e., I s = |E s | 2 , I p is the planet intensity including other light incoherent to the starlight, and ∆E 1 and ∆E 2 are the nonzero modulated EFs. The pair-wise modulations used in the EFC and the SAN, ±∆E 1 and ±∆E 2 , were generated by the ±sin and ±cos shapes of the wavefront at the pupil plane, respectively, to make ∆E 1 and ∆E 2 perpendicular in the complex plane (Figure 1 (a)). Let us define a set as a sequence of the five intensities (I 0 , I + 1 , I − 1 , I + 2 and I − 2 ; see Figure 1 (b)). All intensities are treated as contrasts, normalized by the peak of a point spread function (PSF) of a telescope, in the present paper. After solving the simultaneous equation of the Equations (1) and (2) on E s , the DM is controlled to delete E s by producing −E s at many pixels to make a dark area-thus called dark-field control-and a final measurement is made without the modulations after reducing I s to the lowest possible level. The details of the control are not the interest of the present paper, although the control is recommended (not required) in parallel with the CDI-SAN observation. The solution of E s was not resolved into real and imaginary parts as in the EFC but into ∆E 1 and ∆E 2 , which were used in the SAN and the present method (Figure 1 (a)) as where • is the inner product when the EF is considered as a two-dimensional vector in Cartesian coordinates.
The five intensity measurements should be made within a timeframe where the variation of E s (hereafter dE s ) is negligibly small, or the simultaneous equations will be disturbed. In the case of a ground-based telescope, atmospheric turbulence uncompensated by an adaptive optics (AO) system would be the cause of dE s at the tens of millisecond scale. In a future space telescope, deformation of the telescope optical structure might produce dE s under a long exposure of minutes and hours aiming at detections of faint exoplanets.
To overcome this problem, the CDI-SAN repeatedly utilizes measurements of the set that are faster than dE s and integrates the data for a long interval well beyond the time constant of dE s in postprocessing (Figure 1 (b)). Let us calculate the average of zeromodulation intensities as and those with some modulations as where means the average of all observed sets. Assuming that ∆E 1 , ∆E 2 and I p are constant, the relations were derived. Then, we obtain the averaged speckle intensity as where the square of the absolute value of E s in Equation (3) and ∆E 1 • ∆E 2 = 0 were used. Finally, from Equation (4), we estimate the unknown planet intensity by removing the averaged speckle intensity I s from the averaged zero-modulation intensity I 0 as  which is the essence of CDI-SAN, an exact solution of the planet intensity just derived using only the average values over a long interval, i.e., the five averaged intensities and the two averaged square differences. Using Equations (7) and (8) need not be limited to the star and planet light in varying wave-front environments. This refers to methods to estimate I s (averaged coherent-light intensity) and I p (incoherent intensity) from the integrated values when additional coherent EFs of ∆E 1 and ∆E 2 are available. Equations (7) and (8) are the exact solutions and can provide a contrast of infinity, i.e., I p1 =I p even for 10 −100 , if E s does not change within a set (perfect correlation), where different E s can be acceptable for different sets. Rather less correlation-randomness-with changing E s is the problem. So, in the worst case, the performance of the equations themselves under the no-correlation condition (random change of E s ) will be good to study first. Intensity changes owing to fast dE s during a set produce errors in the simultaneous equations, which can be considered to act as a kind of noise (hereafter speckle variation noise, SVN).
Other than that, the photon shot noise (hereafter PSN) and readout noise (hereafter RON or R) of a detector will affect to the intensity measurements significantly, when photon flux becomes low. Then, the Equation (8) could be developed into another essential equation for CDI-SAN with bias corrections for the mean square difference values in the numerators as where V + 1 , V − 1 , V + 2 , and V − 2 are the noise variances of I + 1 , I − 1 , I + 2 , and I − 2 , respectively. The noise variances would be estimated by additional mean square differ-ences as (10) where I Y Xa and I Y Xb are two half-exposures of I Y X (see Figure 1 (b)), which satisfy where the denominators, 4 and 2, are required because the half-exposures are the contrast normalized by the PSF's peak. It is expected that the averaging process for I p1 or I p2 for a long interval will produce a highcontrast image by suppressing these noises. The amplitude of ∆E 1 and ∆E 2 should be adjusted carefully because the solutions diverge if I + 1 + I − 1 − 2 I 0 or I + 2 + I − 2 − 2 I 0 becomes close to zero deu to noise. Because the noise behaviors in the Equations (8) and (9) are not simple, we need numerical simulations to know their characteristics.
After revealing how the equations themselves performe, we need to investigate some issues in implementation. We know how to generate ∆E 1 and ∆E 2 for a wide area in the focal plane using a DM operation with the shape formed by the sum of many sin and cosine waves with various wavenumbers (Give'on et al. 2011;Oya et al. 2015); however, the target area should be limited to suppress the maximum DM stroke when ∆E 1 and ∆E 2 become large amplitudes for bright speckle conditions. Here, the orthogonality of ∆E 1 and ∆E 2 should also be considered. The contrast estimation with an appropriate spatial distribution of I s is interesting. E s correlations after t 0 as in the atmospheric wave-front variations, are welcome and would be advantageous relative to the no-correlation conditions.

SIMULATIONS AND DISCUSSIONS
In this section the CDI-SAN method will be examined using numerical simulations. The equations are solvable for each pixel so that the performance of the solutions of the equations themselves are independent of the spatial distributions of I s , ∆E 1 , and ∆E 2 . Therefore, it is effective to estimate the standard deviation (SD) and rms (root-mean-square) from the measured pixel intensities calculated in parallel at many pixels within a defined area assuming uniform I s , ∆E 1 , and ∆E 2 to know the noise level of a pixel after an observation. It would not be good to estimate them by combining them with many nonuniformities, which would make the performance of the equations themselves hard to understand.
A measurement sequence is shown in Figure 1 (b). The set consisted of the five intensities, and they are divided into two half-exposures suffixed with a and b; the set then included 10 half-exposures. Let us introduce a time constant, t 0 , and have E s change to a new value every t 0 . The number of sets obtained during t 0 is M , and the observation was made up to N × t 0 . As mentioned in Section 2, a fixed E s does not generate any errors in the equations and can be removed exactly using the present method; then, the strong-correlation condition would be an advantage and the SVN would become small. But in the simulation here, E s after t 0 was assumed to be 100% independent, which would be the worst condition for the equations.
The E s was generated as a focal-plane EF by a fast Fourier transform (FFT) of a pupil function with a wavefront error, the diameter of which was 128 in a 256 × 256 array, where the EF without the wave-front error (the PSF of the pupil) was subtracted as a coronagraph. The pupil function was unity within the diameter of the circular aperture and zero outside. The wave-front phase error (hereafter φ, omitting the pupil-plane coordinates, generated in the full of array) was assumed to have a spectrum of a flat-power up to a radius of 80 followed by a power law of -4 outside, and random arguments, where the complex conjugate numbers were forced at axisymmetric positions to make the phase error a real number. Its absolute amplitude was then adjusted to have an SD of 10 √ C waves in the pupil, where C was the desired mean contrast of a rectangular extracted region (explained later) and C = 1 × 10 −5 (hereafter case 1, for a ground-based telescope). The phase error was applied to the pupil function as e iφ . A wave-front amplitude error was not considered. The subtraction of the PSF after the FFT was equivalent to a subtraction of unity in the pupil before the FFT. Under a small wave-front-error approximation, the residuals after the unity subtraction would become iφ in the pupil area, and then I s would have the same intensity distribution as the spectrum of φ, flat within the radius of 80 pixels and not depending on the distance from the center. Therefore, we can ex-tract anywhere from the flat-spectrum area to study the present method, except for the correlated pixels at the axisymmetric positions. A focal-plane area of 100 × 50 pixels (−50 ≤ x ≤ 49, 11 ≤ y ≤ 60) was extracted. At this point, the average of I s within the rectangular area was about C, then it was adjusted to C by multiplying a correction ratio to E s . In addition to case 1, C = 1 × 10 −9 for a space telescope will be discussed as case 2.
A new independent E s was generated at t 0 intervals and nonlinearly interpolated for the 10M half-exposures between the intervals using weights of (1 − w)/W and w/W for the previous wave-front error and the new one, respectively, where w = k/10M (k = 1 to 10M ) and W = (w 2 + (1 − w) 2 ) 1/2 . The area was, as shown in Figure 2 (a), divided into five regions every 20 pixels in the x direction, (i), (ii), (iii), (iv), and (v), considering a regional contrast coefficient (hereafter B) of 100, 10, 1, 0.1, and 0.01, respectively, to cover wide range contrasts at a time. Now √ B was multiplied to E s , and the contrast was forced to stepwise BC. Figure 2 (a) is an example I s of one exposure without modulation and noise.
A photon flux (hereafter P ) per t 0 at C must be introduced, and P =100 was adopted, then the flux BP =10,000, 1000, 100, 10, and 1 at each region. The photon flux depends on a target star's magnitude, a telescope diameter, and observation conditions, and P =100 would be roughly obtained, e.g., at GJ 411 (M2V, m J =4.2) using a 30m telescope, a wavelength of 1.2µm, and t 0 =36 ms in case 1 or at a solar-like star of 10 pc distance using a 9 m telescope, a wavelength of 500 nm, and t 0 =3600 s in case 2, under the common conditions of a total photon detection efficiency of 0.2, an acceptance of a quarter energy of the PSF at a pixel, and a band width of 2%.
For the modulation EF, constant values of ∆E 1 = A √ I n and ∆E 2 = ∆E 1 i were adopted for each region, where A was a modulation amplitude coefficient, i was an imaginary unit, and I n was an expected noise intensity estimated by I n = (B 2 P 2 + BP + 10M R 2 ) 1/2 C/P considering the SVN, the PSN, and the RON. Here A=2 was adopted as it showed the best contrast among A = 1, 2, 3, 4, 5, and 10 in most conditions. The contrast becomes worse with A less than 1, which will be discussed later.
Half-exposures were made including the PSN and the RON, and the data were integrated. An image of I p2 for the case 1 was shown in Figure 2 (b) under the conditions of M =2, R=0.1, and N =100,000, which corresponded to a 3600 s observation when t 0 =36 ms. We could recognize bias intensities in bright regions, and then the image appeared to be improved by subtracting the bias intensities calculated using the median values at individual regions as shown in Figure 2 (c) and at four horizontal profiles plotted in Figure 2 (d). The artificial planets observed in the figures were included in the halfexposures. The brightest planet in each region had an intensity of BC and the best appeared above the noise levels were down to BC/1000, BC/100, and BC/10, at (i) and (ii), (iii), and (iv) and (v), respectively.
Next, the I p2 was recalculated by adding either the SVN, the PSN, or the RON individually without the planets to separate their contributions to the noise level, and an SD was estimated. An SD of 50 pixels in y at each x (hereafter SD x ) is shown in Figure 3 (a), where I 0 indicates the original contrast in reference to the average of I 0 for the corresponding pixels. The effect of the RON was successfully lower than that of the PSN in all regions. It was found that the effect of the SVN was well suppressed to the same level as the PSN at (i) and lower in other regions, which meant that the present method with M =2 could overcome the dE s problem under the conditions evaluated here. The same figures for various N can be seen in Appendix A.
An SD at each region was then estimated and normalized by I 0 , expressed as SD/ I 0 hereafter, which indicated a contrast improvement (hereafter CI), where I 0 is an average of I 0 at the corresponding region. An rms/ I 0 was also estimated as another CI to obtain the bias intensity, where rms was from each region. The SD of I 0 is equal to I 0 and BC statistically because the speckle intensity followed a chi-square distribution with two degrees of freedom. Figure 3 (b) showed the SD/ I 0 with only the SVN, where the values are the average for all regions this time because they did not depend on the region. The larger M and the larger N have better contrast at least within the ranges shown here. The slope at N ≥ 1000 was almost 1/ √ N , and better slopes were obtained at N ≤ 100 for M ≥ 2. The achievable contrast with M ≥ 2 is obviously better than that with M =1 at N ≥ 100. M ≥ 2 is recommended if other conditions permit. The rms/ I 0 is also plotted in Figure 3 (a), where its saturations are found at N ≥ 100 indicating bias intensities produced by the SVN. That was why bias intensities were observed in bright regions in Figure 2 (b). The bias intensity was positive and decreased with a large M . If the bias intensity can be removed by a smooth surface at the focal plane as seen in Figure 2 (c) and (d), we do not have to use a large M to suppress the bias. A small M is effective to slow the data acquisition speed and to reduce parts of the SVN and the RON as discussed later. Therefore, M =2 was used in the evaluations hereafter. Figure 3 (c) showed SD/ I 0 with only the PSN for the five regions with a different photon flux of BP under M =2. Here the dE s during a set stayed zero, i.e., the wave front changed only at the beginning of each set, to get no SVN even under the wave-front change condition. As reference, the SDs of ( a − b )/ I 0 are also shown, where a and b are averages of all the half-exposures suffixed by a and b, respectively. They lie at 2/ √ BP N , which are considered PSN limits of ideal differential images with halves. The SD/ I 0 of I p2 showed a good performance with a small increase of only 1.2 of the PSN limit at BP ≥ 100 between 10 ≤ N ≤ 100, 000, where almost no M dependence was observed. At low-light levels, the increase is still only about 1.3 and 1.8 times at BP =10 and BP =1, respectively, and gradually rises with M . The rms is almost the same as SD which meant that the bias was well removed by the correction terms in Equation (9).  Figure 3 (d) shows the SD/ I 0 with only the RON for the five regions under M =2, where R=1 was used to see whether it is worse than R=0.1. The lines at 5 ≤ N ≤ 100, 000 were almost coincident with the SD of ( a − b )/ I 0 calculated as in the PSN, which lies on ideal RON limits from the function 2R √ 10M /BP √ N . Even with R=0.1 and 10 and M =1-5, the results were close to the function (0.8-1.2 times), e.g., the condition BP/R=1/0.1 was comparable to 10/1. The rms was almost the same as the SD, and again, the bias was well removed by the correction terms in the Equation (9). The successful bias corrections for the RON and the PSN were the reason why Figure 2 (b) showed no bias in the low-light regions. The CIs of I p1 without the bias corrections, Equation (8), are mentioned in Appendix B because they were not the point of the present paper. Now the CIs of the CDI-SAN method against the three noises, which would be considered independent, have been resolved. Also, the CIs in the high-contrast case 2 (P =100 at C = 1 × 10 −9 ) turned out to be the same as those in case 1 described above (Figure 3 (b), (c), and (d)) because they depend on P instead of on C, although N should just be 100 (see the SD x in Figure 9 (a) and (b) of Appendix A). An important result was that CDI-SAN could approach the natural limit of the PSN by suppressing the SVN and the RON with appropriate M and R under the available conditions of (B)P and N .
The orthogonality of the ∆E 1 and ∆E 2 is important to achieve high contrast. An error of the relative argument between them, hereafter θ, would produce an increase in the SD. Assuming that arg(∆E 1 ) = 0 • and arg(∆E 2 )=arg(∆E 1 )+90 • + θ, the SD/ I 0 of I p2 was calculated and shown in Figure 4. It was found that the θ increased the SDs for all the noises. The log(CI) was limited to about -3.6 as seen at the region (i) for θ = 4 • . The orthogonality of the ∆E 1 and ∆E 2 should be considered according to the raw contrast and the target contrast. A similar result was confirmed with arg(∆E 1 )=45 • and arg(∆E 2 )=135 • which will be used later.
When θ is known given by a model, it is possible to use a solution without the orthogonality assumption-a free-argument solution-which is shown in Appendix C. It would be useful if the difference between the model and real optics is small.
The modulation amplitude coefficient, A, would affect the SD. The A-dependent SD/ I 0 of I p2 relative to A=2 is shown in Figure 5 for M =2 and N =100,000. It strongly affects the SD when it is less than 1, which depends on about 1/A 2 , while almost no effect is confirmed between 1 and 10. We recommend using A = 1 to ∼ 2 because a very bright modulation would cause a problem in real optics. Re ( E1) Im ( E1) Re ( E2) Im ( E2) Figure 6. Generation of the modulation EF in arbitrary units. The upper panels are the initial desired EFs, the middle panels are the second desired EFs after the pupil function and the apodization are operated on the pupil plane, and the bottom panels are the EFs generated by the operation e iφ in the pupil and FFT. The above describes the characteristics of the solution of the equation itself. Next, generating the modulation EFs with the DM in real optics will be discussed. A candidate of a modulation wave front (hereafter MWF), i.e., phase map, at a pupil plane can be calculated by the inverse FFT of the desired EF, ∆E 1 or ∆E 2 , at a target area in the focal plane (e.g., Give'on et al. 2011). Let us determine the amplitude and argument of the initial desired modulation EF, and then calculate the real and imaginary parts (see Figure 6 top panels). The EF is the negative conjugate at each axisymmetric position. Let us adopt arg(∆E 1 )=45 • and arg(∆E 2 )=135 • to retain good or-thogonality even at a large amplitude of the modulation EF. After an inverse FFT operation, only the imaginary part is generated. A pupil function and an apodization function are operated on (see Figure 7, top panels), where the latter of which is used to get a flatter amplitude of the modulation EF, and the candidate value of the MWF is obtained in the imaginary part (see Figure  7, bottom panels). The FFT operation of this function generates the second desired modulation EF (see Figure  6, middle panels). Finally, the candidate value is used as the phase of the MWF, and FFT is used to obtain the (final) modulation EF (see Figure 6, bottom panels). The focal-plane modulation EF can be generated by an FFT of the pupil function with the MWF and subtracting the unaberrated pupil function. When the MWF is very small, the (final) modulation EF is almost equal to the second desired modulation EF.
The generated modulation EFs whose desired intensities are 10 −6 , 10 −5 , 10 −4 , and 10 −3 are shown in Figure 8 and characteristic values are listed in Table 1. These are the candidates with relative arguments within 90 ± 3 • .5, taking into account the discussion above, and whose amplitude degradation was better than 0.74. Here an FFT array of 512 × 512 was provided and the diameter was 256 pixels; then, 2 pixels/λ/D was obtained. The modulation area had a radius of 20 pixels (10 λ/D) and |x| ≥ 2 pixels (1 λ/D) for the two darker contrast cases, and others can be seen in Table 1. The modulation EF was evaluated 2 pixels inside the modulated area in all cases, although pixels outside the evaluation area were  Figure 8. Modulation EFs generated by modulation wave-fronts. The top panels are the amplitude and the bottom panels are the relative argument between ∆E1 and ∆E2. Desired contrasts are 10 −6 , 10 −5 , 10 −4 , and 10 −3 from left to right. The modulation area has a radius of 20 pixels (10 λ/D) and |x| ≥ 2 pixels (1 λ/D) for the two darker contrast cases. The color map shows the evaluation area, which is 1 λ/D smaller than the initial desired modulation area. Detailed values are listed in Table  1. not absolutely unusable. Even for the large, nearly 1 radian wave-front phase, which is a difficult condition to approximate in the imaginary part, the relative argument between ∆E 1 and ∆E 2 stayed around 90 • by setting arg(∆E 1 )=45 • and arg(∆E 2 )=135 • , although they changed by more than a dozen degrees at maximum from the initial desired modulation EF. Setting arg(∆E 1 )=0 • and arg(∆E 2 )=90 • did not successfully keep the relative argument at 90 • in the desired intensities above. Because the error of the relative argument and the decrease in the modulation amplitude hinder the improvement of the contrast, careful selection of the modulation area and its amplitude is required. The modulation EF should be arranged while considering observational conditions such as raw contrast, desired area size, target contrast improvement, and total integration time. The MWF for the arranged modulation EF will be added to the AO-corrected wave front in real observations. In the simulation, a perfect coherence of the speckle EF was assumed; however, the incoherent effect of band-width, even at 2%, on the measurements should be studied in a future paper. The stability of the modulation amplitude, incoherent scattered light, and speckle (wave-front) variation characteristics should be studied in real optical conditions. On the other hand, compositions of multiple-band data might help increase the achievable contrast. Combining CDI-SAN with other processing methods might be also interesting. A simulation using a phase screen model is in progress, which would be of advantage for the correlation problem.

CONCLUSION
The CDI-SAN method was developed to detect faint exoplanets lying beneath varying residual speckles. CDI-SAN can be applied to both ground-based and space telescopes. It utilizes image acquisitions faster than stellar speckle variations synchronized repeatedly with five shapes of a deformable mirror. By using only the integrated values of each of the five kinds of images and several square differences for a long interval of ob-servations, the light of the exoplanet can be separated from the stellar speckles. The achievable contrast would reach to almost the photon-noise limit of the residual speckle intensities under appropriate conditions, which would be helpful for exoplanet research missions today and in the future.
This work was partly supported by JSPS KAKENHI grant Nos. 19H01932 and 20H05893. The author would like to thank Dr. N. Murakami for helpful discussions and Profs. S. Tsuneta, Y. Hayano, T. Usuda, and M. Tamura for constant encouragements. APPENDIX A. THE SD X OF I P 2 FOR VARIOUS N Figure 9 showed the SD x of I p2 , as in the Figure 3 (a), with one of the SVN, the PSN, or the RON individually for N =10, 100, 1000, and 10,000 in (a), (b), (c), and (d), respectively, where M =2 and R = 0.1 were used. The CIs are understood to be the SD x levels relative to I 0 at the corresponding regions. The left ordinate indicated the results for the case 1. Figure 9 (a) and (b) also show the contrast in the case 2 in the right ordinates. The CIs in this case 2 are the same as those in the case 1 because they depended on P instead of C, where only N should be 100 if t 0 =3600s and a limited total integration time of 100 hr.

B. CHARACTERISTICS OF I P 1
The characteristics of I p1 , Equation (8), are presented here. The solution I p1 show the bias intensities of the PSN and the RON at low-flux regions, which are seen in saturations of the RMS/ I 0 in contrast to I p2 of Equation (9). They remain because the second and the third terms of the equation of I p1 are larger than those without the noises, and then the uncorrected biases of I p1 by the two noises are negative. Figure 10 (a) shows the rms/ I 0 of I p1 for the PSN, which saturated at a much larger level than the SD/ I 0 levels at all regions. Figure 10 (b) showed the rms/ I 0 saturations of I p1 for the RON, which were found specifically at the low-photon-flux regions.
However, the SD/ I 0 for the PSN shown in Figure 10 (a) was almost the same as I p2 in Figure 3 (b), about 1.2 times ( a − b )/ I 0 at BP ≥ 10 between 10 ≤ N ≤ 100, 000. At a very low photon flux of BP =1, the increase is still only about 1.6 times at M =2 and gradually rises with M , which is slightly better than I p2 . The SD/ I 0 for the RON (R=1) in Figure 10 (b) was almost the same as I p2 in Figure 3 (c) and close to the RON limit of ( a − b )/ I 0 at 50 ≤ N ≤ 100, 000. The SD of I p1 against the RON can be improved by 1/ √ 2 by not dividing the five intensities into two half-exposures, if available.
The biases should be subtracted in some way to know the intensity of the planets. If the bias intensities of the PSN and the RON can be removed by a smooth surface at the focal plane, as in the SVN case shown in Figure 2 (c) and (d) in a real optical condition, we do not have to use I p2 to remove the biases but can use I p1 , which has slightly better SDs in some conditions and has the advantage of slow data acquisition without splitting into two half-exposures. If surface subtraction is not available, it would be good to use I p2 to escape from the bias problem. It was found that the biases in I p1 from the RON could be reduced with a large A, and the large A ≤ 10 did not affect the SDs of the three noises by more than 20% in the simulation. Large A introducing bright modulation intensities, however, might have a risk of additional measurement errors in practice.
On the other hand, on the SVN, the RMS/ I 0 and the SD/ I 0 of I p1 were almost the same as those of I p2 shown in Figure 3 (b). The rms of I p1 was slightly small, about 0.9 times, I p2 at 1000 ≤ N ≤ 100, 000. For either I p1 or I p2 bias intensities exist, as shown in Figure 2 (b), which should be removed by the smooth surface as shown in Figures 2 (c) and (d) or by a large M . The contrast achievable by I p1 and I p2 should be studied in real optical conditions.

C. FREE-ARGUMENT SOLUTION
The solutions, equations (7) and (8), were derived by assuming the orthogonality of the ∆E 1 and ∆E 2 . The error of the relative argument between them, θ, seen in Figure 11, would produce an increase in the SD. When θ is known,    Figure 11. Electric field at a pixel with the relative argument error θ.
given by a model, it would be possible to use a solution without the orthogonality assumption. By using ∆E 2 with θ instead of ∆E 2 , the solution can be written as