{'fraction_non_alphanumeric': 0.029528724934427618, 'fraction_numerical': 0.011591505203485912, 'mean_word_length': 4.628333333333333, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Point Process Algorithm: A New Bayesian Approach for Planet Signal Extraction with the Terrestrial Planet Finder

The capability of the Terrestrial Planet Finder Interferometer (TPF-I) for planetary signal extraction, including both detection and spectral characterization, can be optimized by taking proper account of instrumental characteristics and astrophysical prior information. We have developed the Point Process Algorithm (PPA), a Bayesian technique for extracting planetary signals using the sine-chopped outputs of a dual nulling interferometer. It is so-called because it represents the system being observed as a set of points in a suitably-defined state space, thus providing a natural way of incorporating our prior knowledge of the compact nature of the targets of interest. It can also incorporate the spatial covariance of the exozodi as prior information which could help mitigate against false detections. Data at multiple wavelengths are used simultaneously, taking into account possible spectral variations of the planetary signals. Input parameters include the RMS measurement noise and the a priori probability of the presence of a planet. The output can be represented as an image of the intensity distribution on the sky, optimized for the detection of point sources. Previous approaches by others to the problem of planet detection for TPF-I have relied on the potentially non-robust identification of peaks in a"dirty"image, usually a correlation map. Tests with synthetic data suggest that the PPA provides greater sensitivity to faint sources than does the standard approach (correlation map + CLEAN), and will be a useful tool for optimizing the design of TPF-I.


Introduction
A key component of the proposed Terrestrial Planet Finder mission is a nulling interferometer (TPF-I) for the wavelength range 7-15 µm, currently envisaged as a free-flying four-element dual-Bracewell array in the "X" configuration as shown in Figure 1. The signals from the two elements in a nulling pair are combined with π phase shift, resulting in destructive interference for a source on the optical axis, i.e., the axis corresponding to equal geometrical path lengths to the two nulling elements. By placing this null at the location of the star, the contrast between planets and star is dramatically increased, thereby overcoming one of the chief obstacles to planet detection (Bracewell 1978). The starlight suppression is further improved by using two such nulling pairs and combining their outputs coherently (see, for example, Beichman & Velusamy (1999)). As the interferometer system is rotated about the line of sight, a modulated signal is then produced, and this signal contains information on the intensity distribution of sources on the sky. Its information content can be increased by alternately applying phase shifts of π/2 and 3π/2 and taking the difference of the phase-shifted signals. This has the effect of suppressing instrumental effects such as low-frequency detector noise and thermal background drifts, and is commonly referred to as "sine chopping" ( (Woolf et al. 1998;Velusamy & Beichman 2001)). Sine chopping also suppresses any structure on the sky which is symmetrical about the star, including that from symmetrical structure in the exo-zodiacal (exozodi) cloud, thus increasing the contrast for the detection of planetary signals.
Since the interferometric signals from such a small number of baselines provide very limited information on spatial Fourier components of the intensity distribution of sources on the sky, it is necessary to supplement the observations with a priori information in order to recover that distribution. One widely-used procedure for this purpose is the maximum entropy method (see Narayan & Nityananda (1986) for a review, and Sutton & Wandelt (2006) for recent developments). However, such procedures are not optimal for the present problem because they do not incorporate an important piece of prior knowledge, namely the pointlike nature of the planetary sources.
Of the various inversion procedures that do incorporate such information (see, for example, Angel & Woolf (1997); Velusamy & Marsh (2004)), the one most widely used involves making a correlation map (referred to as a "dirty image") followed by deconvolution with the CLEAN algorithm (Högbom 1974;Draper et al. 2005). A key step in that procedure is an iterative search for peaks in the dirty image, which can be non-robust when noise bumps fall on "side lobes" of the responses to other sources.
In this paper we propose a Bayesian technique for planet detection which avoids the noise-vulnerable peak-finding step. It can process data at many wavelengths simultaneously, and can incorporate prior information such as the spatial covariance of the exozodi. Since the planetary signals are represented as a random set of points in a suitably-defined state space, we refer to the resulting algorithm as a "Point Process Algorithm" (PPA) following the terminology of Richardson & Marsh (1992).
The PPA, whose mathematical foundation is presented in detail by Richardson & Marsh (1987), is a general technique of which the TPF problem represents a specific application (Velusamy et al. 2005). It may be regarded as a state estimation technique in which the system to be described (for example, an image) consists of the superposition of an indefinite number of objects, each of which is characterized by a set of parameter values. Each object can thus be represented as a single point in a multidimensional state space of which each axis corresponds to one of the parameter values. In our application, the axes of the state space are simply the x and y positions and flux density of an individual point source.
Other possible examples are: (1) the representation of an image of a cluster of galaxies in terms of a 6-dimensional state space in which the axes represent the x and y location, the flux density, major and minor axes and orientation of an individual galaxy, and (2) the representation of a planetary system using a state space whose axes consist of the orbital parameters characterizing an individual planet, using a set of measurements of Doppler shifts or astrometric measurements of the parent star (or both). In each of these examples, the goal of the PPA is to estimate the most probable set of points in the state space given the set of measurements.

Measurement Model
The starting point of our approach is a measurement model which relates a data vector, d (whose components are the complete set of samples of the sine-chop signal at all wavelengths of observation), to the intensity distribution, I(x, y), on the sky. The latter is modeled as the superposition of a set of point sources of unknown number, fluxes and positions, upon an extended background whose intensity at position (x j , y j ) is denoted by ζ j . The distribution of point sources is represented as a set of occupation numbers in a 3-dimensional state space whose axes are flux and x, y position. The state space is divided into a regularly-sampled grid of cells, such that the cell with coordinates (f k , x j , y j ) represents the kth possible flux value, f k , at the jth spatial position, (x j , y j ). If a point source of that particular flux density is present at that particular position, then the occupation number of that cell, Γ n (where indices j and k have been mapped onto a single index, n), will be equal to 1; otherwise, it will equal zero. This representation provides a convenient framework within which to incorporate our prior knowledge of spatial dilution.
We assume that the source distribution is being observed by a rotating dual-Bracewell interferometer consisting of two nulling pairs separated by an "imaging baseline," and that the measurements consist of a time series of sine-chop signal values as the system rotates with respect to the sky. If the ith such measurement is made at an orientation corresponding to instantaneous baseline vectors b i and B i for the nulling and imaging baselines respectively, then the instantaneous response to a point source of unit flux at off-axis sky position (x j , y j ) (denoted by unit vector s j ) is given by: where R i (s j ) is the "primary beam" response, i.e., the response of an individual detector, and the baselines b i and B i are in units of wavelengths.
The measurement model must include a spectral model which relates the planetary fluxes at the different wavelengths. There are several possibilities, one of which is to treat the planets as black bodies at the local radiation temperature. If the orbital inclination and position angle of the orbital tilt axis are known, the planet temperature, T j ≡ T (x j , y j ), can then be determined at each (x j , y j ) from knowledge of the stellar luminosity and radial distance. Our measurement model can then be written as: where λ i represents the wavelength at which the ith measurement is made, B (λ, T ) represents the Planck function, and the f k represent possible values of planetary flux at some suitably-defined reference wavelength, λ 0 , and ν i is the measurement noise. The index, n, is understood to be a function of j and k.
It will be convenient to lump both of the last two terms on the r.h.s. of (2) into a single "noise" vector, µ, whose components are given by: If we then define a matrix A whose components are then we can rewrite (2) in matrix notation as: where Γ is the state vector whose components are Γ n .

Statistical Assumptions
The measurement noise in the TPF signals is distributed as a Poisson process, dominated by the photon statistics of stellar leakage and the local zodi background. Since the expected number of photons per integration period is large (typically ∼ 10 5 ), the statistical distribution accurately attains the limit corresponding to a Gaussian random process (GRP). The measurement noise, ν i , on the sine-chop signal can then be regarded as a zero-mean GRP for which: where E is the expectation operator and δ is the Kronecker delta.
For present purposes we regard the exozodi as an error source rather than as a quantity to be estimated. We model it only to the extent necessary to mitigate its effects on planet detection. If it is distributed symmetrically about the star, then its only effect is to increase σ ν by a factor corresponding to the square root of the increase in the overall photon count, since its coherent contribution cancels out in the sine chop signal. If the exozodi is not symmetrical, one approach for dealing with it is to subtract its estimated contribution from the measurements ahead of time, based on available data from previous (possibly groundbased) observations at lower resolution. The error resulting from this subtraction may be regarded as an additional zero-mean error contribution to our measurement model. It will, in general, be spatially correlated as a result of the finite resolution of the observations and the intrinsic spatial correlation properties of the exozodi itself. We therefore approximate the residual exozodi, ζ j (λ i ), as a GRP with covariance C ζ defined by: From (6) and (7), we can then show that the covariance of µ is given by: The state vector, Γ, is also regarded as a stochastic process; its a priori probability distribution is assumed to be described by: where The quantity P 1 thus represents the a priori probability of occupancy of a cell in the state space of flux and position.

Estimation Procedure
The central operation is to estimate the state vector, Γ, given the observed data. The details of the estimation process depend on the type of optimality criterion selected. For the present algorithm we choose a mean square error criterion, for which the optimal estimate is the a posteriori average value of Γ, given by: where z n is a 3-d vector representing the coordinates, (x j , y j , f k ), of the nth cell in state space, and the conditional probability, P (Γ|d), is given by Bayes' rule: in which P (Γ) is given by (9), P (d) is a normalization constant, and We refer to ρ (z n |d) as a density, since it represents the average local density of occupied cells in the state space of position and flux. Its estimation is a generic problem in statistical mechanics; previous applications have included acoustical imaging (Richardson & Marsh 1987) and target tracking (Richardson & Marsh 1992). The solution procedure involves the solution of a hierarchy of integro-differential equations which fortunately can be truncated, to a good approximation, at the first member. We then obtain: where t is a dimensionless progress variable representing the degree of conditioning on the data, and φ 1 is the conditioning factor given by: where b represents a vector formed from the diagonal elements of A T C −1 µ A. A discretized version of (14) is formed by regarding the process as a series of N c weak conditionings, each of which corresponds to a measurement noise covariance of N c C µ . The progress variable, t, then corresponds to the number of conditionings, and the dt in (14) is replaced by 1. Equation (14) is then to be integrated from t = 0 to a terminal value, t f = N c , representing full conditioning on the data.
The initial condition on ρ for the numerical solution of (14) is ρ (t=0) = ρ 0 , where ρ 0 is the a priori density equal to the constant value P 1 ; its sum over all state space is equal to N 0 , the a priori expectation number of planets present. If we regard the true number of planets, N, as an unknown, then this quantity could, in principle, be estimated by adjusting P 1 for equality between N 0 and the a posteriori number of planets, corresponding to the sum of ρ(z n |d) over all state space (Richardson & Marsh 1987). In practice, however, this would be numerically cumbersome. We have therefore adopted an approximate alternate procedure in which the number of planets is estimated with the aid of χ 2 ν , given by: where M is the total number of data samples. In this procedure, we fix P 1 at a very small value, thereby producing an a priori bias towards a small number of planets being present. We then proceed with the numerical integration (during the course of which, χ 2 ν decreases monotonically), terminating it when χ 2 ν = 1. We thereby arrive at a solution in which the data are fit by the minimum possible number of planets.
Our corresponding estimate of the source intensity distribution is then: Estimates of the planet fluxes themselves may be obtained from the integrated value around each peak in this image; the uncertainties correspond to the standard a posteriori variances of maximum likelihood estimates.

Tests with Synthetic Data
We have tested the PPA using a set of 15 test cases designed to replicate realistic observing scenarios for TPF-I in the wavelength range 7-15 µm. Each case involved 0-5 planets at radial distances of 0.4-5.25 AU from a solar-type star at a distance of 10 or 15 pc, superposed on a 1 Zodi dust distribution. The 15 cases involved a total of 37 planets, 24 of which were of 1 Earth flux, and 2 of which were < 1 Earth flux. The orbital inclination was 60 • except for two face-on cases. All planets were at Earth temperature (260 K). Radial distances and signal-to-noise ratios were distributed as shown in Figure 2.
The measurement configuration was the "X array" consisting of four 4-m detectors at the vertices of a rectangle of width 12 m (corresponding to the nulling baseline length) and a length of either 36 m or 72 m, as shown in Figure 1. We adopted the 12 m nulling baseline in lieu of the nominal design value of 20 m to improve the sensitivity at the short wavelengths. For longer baselines (20 m in particular), the stellar leakage becomes larger because of the narrower null, with particularly severe effects at the shorter wavelengths.
The measurements consisted of the simulated sine chop signals through a full rotation of the array, at wavelengths of 7. 44, 8.50, 9.92, 11.90 & 14.90 µm, with bandwidths ranging from 0.9 µm at 7.44 µm to 3.7 µm at 14.9 µm, and a total integration time of 1 day. Superposed on these signals was Poisson noise due to the various instrumental and astrophysical effects discussed by Beichman & Velusamy (1999), the most important ones being stellar leakage (due to incomplete nulling of the stellar photosphere) and thermal emission from the localand exo-zodiacal dust clouds, the latter assumed geometrically symmetrical.
The data were inverted using the PPA, and also by the standard technique of "correlation map + CLEAN" for comparison purposes. The spatial covariance of the exozodi was ignored (i.e., C ζ was set at zero) since symmetrical sources cancel in the sine chop signal. Full account was, however, taken of the Poisson noise contribution of this component. Figure 3 shows the results for five representative cases involving both array configurations. The results of the PPA and CLEAN algorithms are shown for the 36 m and 72 m configurations; note that in the case of CLEAN, only the 6 strongest source components are shown. Also shown (Column 4) is the result of using the PPA with a combination of data from both configurations, but with the integration time split equally between the two, i.e., the integration time was 0.5 day per configuration. Figure 4 shows the results for the case involving 5 assumed planets (Case 4 of the previous figure). Also shown for comparison (panel (e) of Figure 4) is the result of using noiseless data; it demonstrates that the fidelity of the PPA image is limited, in part, by the gaps in spatial frequency coverage of the measurements.
For all 15 test cases, these inversions yielded images from which we extracted estimates of point source locations, fluxes, and corresponding uncertainties. We identified the six brightest peaks in each image, then applied a SNR threshold and counted the detections. Comparison with the true (assumed) planet positions then led to the number of true detections and false alarms as a function of the detection threshold in sigmas. The fluxes of detected sources were recovered to an accuracy consistent with their theoretical uncertain-ties. For the 5-planet case above, with the 36 m array, these uncertainties ranged from approximately 2% for the brightest source (8 Earth fluxes) to 20% for the weakest (1 Earth flux). Figure 5 shows a plot of true positives vs. false negatives, known as the operating characteristic curve of the detection system, for both PPA and CLEAN. Also included in the figure is the theoretical operating characteristic curve for an isolated source observed with the 36-m array, assuming pure Gaussian statistics. It is apparent that the PPA detected significantly more planets than CLEAN, particularly in the case of the fainter objects at low S/N thresholds. While it is also apparent that the PPA results fall below the theoretical curve, it should be borne in mind that the latter did not take into account the interactions between sources due to the incomplete spatial frequency coverage. Figure 6 shows the PPA performance as a function of detection threshold in sigmas for the 36-m array.

Discussion
The PPA is theoretically a near-optimal approach, and the present results support this. They demonstrate that use of the PPA increases the effective sensitivity for planet detection over standard techniques exemplified by CLEAN.
The results also indicate that the 72 m array would be less sensitive than the 36 m array for this ensemble of planets, presumably due to the smaller relative range of spatial frequencies involved. For this reason, it would be more effective, for a given amount of integration time, to use the 36 m array exclusively rather than split the time between the 36 m and 72 m configurations, even though the inclusion of the 72 m array would increase the overall spatial frequency coverage. The fact that this added coverage did not offset the reduced integration time of the 36 m array is apparent from the results in Figure 3, and is confirmed by the statistical results from the complete set of 15 cases.
Future work in the application of the PPA to the planet detection problem will emphasize the effects of exozodi, including the effects of spatial inhomogeneities which may masquerade as planets in the interferometer signals, and their mitigation via appropriate incorporation of a priori estimates and spatial covariance. For the latter purpose we will make use of recent observations of debris disks using the Spitzer Space Telescope (Beichman et al. 2005).
We will also explore some of the possible applications of the PPA beyond the TPF problem, since, as we have discussed, it is not limited to the detection of point sources. It may, in fact, be regarded as a generic technique for modeling a set of data (which may or may not be image-like) in terms of user-defined primitives (i.e., objects) which are parametrizable.
The main requirement is that the superposition rule be satisfied, i.e. that the measurements represent a linear combination of the contributions from the individual objects in the system. The objects themselves need not be defined in terms of astrophysical entities; they may instead be a set of parametrizable fractal elements convenient for representing an observed image. In such an application, the goal of the PPA would be similar to that of adaptive scale pixel deconvolution techniques such as the Asp-CLEAN technique of Bhatnagar & Cornwell (2004) and the "pixon" technique (Dixon et al. 1996). Unlike the latter, however, the PPA is not restricted to cases involving PSFs of finite support, and it has some potential advantages over the former. One of these is that it may provide a more efficient procedure for spatially-complex images, since a consequence of the PPA's use of an occupation-number representation is that the computational burden does not increase significantly with the number of source components in the image. The second potential advantage of the PPA is that it is fully capable of superresolution, which in practical terms means that it can separate adjacent point sources whose separation is less than the PSF width. It is not clear that Asp-CLEAN would have this capability given the fact that the standard CLEAN algorithm does not. We will investigate these and other issues in order to better characterize PPA's niche in astronomical problems. (a) Assumed geometry: Two 12 m nulling pairs are separated by 36 m or 72 m. Within each nulling pair, the signals are combined with a relative phase shift of π, resulting in destructive interference at the location of the star. The signals from the two nulling pairs are then combined coherently after introducing relative phase shifts, alternately, of π/2 and 3π/2; the difference between these two signals then constitutes the "sine-chop" signal. In practice, the required pair of sine-chop differencing signals is obtained by combining the beams from all four elements simultaneously in a central beam combiner, after introducing a set of phase shifts, to the individual elements, of either {0, π, π/2, 3π/2} or {0, π, 3π/2, π/2}. (b) Angular response on the sky when combining the two nulled signals with π phase shift; the star (nulled out) is at the center, and the planet position is indicated by the white disk. A modulated signal is produced as the interferometer system rotates and the planet passes through peaks and valleys of the response. (c) Example of modulated signals at 12 µm wavelength simulating an observation of a planet 3 Earth radii in size, 1 AU from a G2 star at 10 pc, using an X-array with B = 36 m and b = 9 m.