pointwise product is always lower. Conversely, if there is a significant
overlap between the results of the two convolutions, then the objective
function $\Omega(\mathbf{Z})$ is significantly higher.

The spatial spectrum can be estimated by computing the DFT of the
autocorrelation function of the corresponding signal. It is worth noting that
for a given sparse configuration, the autocorrelation function of the
selection vector is given by the corresponding redundancy of the
autocorrelation lags. Therefore, unlike the structured sparse array design
which seeks to maximize the contiguous correlation lags, the MaxSINR sparse
design is directly linked to the DFT of the autocorrelation sequence of the
lag redundancy.

We illustrate the proposed approach with the help of the following example.
Consider an 8-element sparse array on the 24 point equally spaced grid
locations that are potentially available for sensor selections. The minimum
spacing among the sensors is $d=\lambda/2$. Consider a source signal located
at $60^{0}$ and four unwanted intereferneces located at $154^{0}$, $55^{0}$,
$117^{0}$ and $50^{0}$ with the INRs ranging from $10-20$ dB. The sparse array
configuration achieving the best performance is found through enumeration and
shown in Fig. 7.4. The associated correlation lag redundancy of this
configuration and the corresponding spatial spectrum is depicted in Figs. 7.5
and 7.6, respectively. The spatial spectrum of the desired source at $60^{0}$
is depicted in the Fig. 7.7. The resultant normalized spectrum of the two
convolved spectra in Figs 7.6,7.7 is shown as solid lines (blue) in the Fig.
7.8. The normalized spectrum for each interfering signal is shown as dotted
lines in the same figure. Fig. 7.9 plots the normalized spectrum for the
sparse array worst case scenario for comparison purposes. Note the maximum of
the convolved desired signal spectrum is at 35th DFT bin. For the best case
scenario, all the convolved interfering signals assume minimum power at the
aforementioned DFT position. This is in contrast to the worst case scenario
depicted in Fig. (7.9) where the convolved interfering signals have
considerable power at the maxima of the desired source. Apart form the maxima
location, it is noted that for the best case design, there is minimum
overlapping between the desired signal and interfering signals at the DFT bin
locations which is clearly in contrast to the worst case design.

In order to further understand the effectiveness of the proposed approach,
Fig. (7.10) plots the SINR performance of all possible sparse configurations
after sorting the array topologies in ascending order of the output of the
proposed objective function. It is clear that the average SINR in the plot is
higher and more desirable towards the left side where the objective function
is minimum vis a vis the right side where the objective function is high. It
is also noted that the best enumerated result of MaxSINR performance does not
correspond to the array configuration with the smallest objective function.
This is because, the optimum sparse configuration also depends on the
beamformer weights to minimize the interfering signals. This is also clear
from the high variance of the curve. Similarly, the worst performing array is
very close to the right side of the curve where the proposed objective
function is high due to the significant overlap of the desired and undesired
signal spectra. The significantly reduced performance towards the right side
of the plot can now be explained intuitively in terms of the inherent
inability of these configurations to mitigate the significant overlap in the
desired and interfering signals which is an artifact of the sparse beamformer
explained in the frequency domain.

Therefore, in light of the above discussions, it becomes prudent to seek
sensor arrangements that minimize the proposed objective function in an
attempt to find more suitable sparse configuration for MaxSINR design. Towards
this end, we propose an iterative algorithm that minimize the proposed
objective function by successive sensor selection, hence deciding on one senor
location at a time. For the initial iteration, a sensor location is chosen
randomly on an $N$ grid points. For the $i$th subsequent iteration, the
proposed objective is evaluated at the remaining $N-i$ locations and then
selecting the sensor location that yields the minimum objective function. The
procedure is iterated $P$ times for selecting the $P$ locations from $N$
possible locations. Due to high variance of the curve in the above figure, it
is best to initialize the algorithm with different sensor location or DFT size
and find the corresponding configurations for each initialization, and
eventually consider the one with the best SINR performance. Steps of the
algorithm are detailed in Table 7.1.

Table 7.1: SBSA Algorithm

0: $N$, $P$, Look direction DOA $\theta$, Interfernce DOAs, SNR and INRs

0: Sparse beamformer $\mathbf{w}_{o}$

Initialize $\mathbf{z}$=$[0\,\,1\,\,0\,...\,0]$ where all entries of
$\mathbf{z}$ are zero except an arbitrarily selected entry

Compute the spatial spectrum of desired source and interfering signals

For (j=1 to $P$-1)

For (i=1 to $P$-1-j)

Select the ith sensor from $P$-1-j remaining locations

Compute the lag redundancy of this sparse array consisting of j+1 sensor

Compute the spatial spectrum of the j+1 sensor sparse array

Convolve the spatial spectrum of the j+1 sensor sparse array with the spectrum
of the desired source and the interfering sources

Compute the overlapping power in the spatial spectra by computing the proposed
metric in Eq. 7.8 EndFor

Select the $i$th sensor from the inner for loop which results in the minimum
overlapping power computed by Eq. 7.8

Update $\mathbf{z}$ by setting the jth location in $\mathbf{z}$ to 1 EndFor

After finding the sparse configuration find $\mathbf{w}_{o}$ by running Eq.
7.5 for reduced size correlation matrix while ignoring the sensor locations
corresponding to $\mathbf{z}$

####  7.3.2 DNN based learning of the SBSA and enumerated designs

Modeling the behaviour of optimization algorithms by training a DNN is
justifiable, as DNNs are universal approximators for arbitrary continuous
functions and have the required capacity to accurately map almost any
input/output formula. For effective learning, it makes sense that the DNN can
generalize to a more broader class that is represented by the finite number of
drawn training examples. From Capon beamforming perspective, it is instructive
to realize that a given arrangement of a desired source direction,
interference DOAs and respectives SNR/INRs constitute one particular example.
The class, in this case, would constitute any arbitrary permutation of the
interference DOAs and respective powers while keeping the desired source DOA
fixed. The DNN task is, therefore, to learn from a dataset, characterized by a
set of different training examples and corresponding optimum sparse array
predictions.

Here, an important question arises that whether we could aim for a more
stronger notion of generalization, i.e., instead of training for a fixed
desired source, is it possible to generalize over all possible desired source
DOAs. To answer this query, we note that for a given desired source and
interference setting, the received signal remains the same even if one of the
interfering signal is assumed as a desired signal and the desired signal is
treated as an interference. In this case, although the received signal is
identical in both scenarios yet the Capon beamformer task is flipped and is
now required to be pointing to an entirely different direction. In so doing,
the corresponding optimum sparse configuration and beamformer weights could be
very different than before. Therefore, instead of relying entirely on the
information of the received data correlation it is imperative to incorporate
the knowledge of the desired source or look direction, which is always assumed
in Capon beamforming formulation. For DNN learning, this information can
either be incorporated by exclusively training the DNN for each desired source
DOA or the desired source DOA can be incorporated as an additional input
feature to DNN. In this chapter we adopt the former approach.

Our proposed approach uses a fully connected neural network with the input
layer of size $2N-1$ and the output layer of size $N$. Although there are $N$
unique correlation lags corresponding to the $N$ sensor locations on the grid,
the dimensionality of the input layer is $2N-1$. This is due to concatenating
the real and imaginary entries of the generally complex valued correlation
lags except the zeroth lag. We have 3 hidden layers with 450, 250 and 80
nodes, respectively. The ReLU activation function is used for all hidden
layers activation. The network is input the correlation values of the received
data assuming a stationary environment. The network output is a binary vector
such that 1 indicate sensor selection and 0 indicates the absence of the
corresponding sensor location. For the scope of this chapter, we assume that
we have an estimate of all the filled co-array correlation lags corresponding
to the correlation matrix of the full aperture array.

The received data is generated in the following manner. For a given desired
source location, the $i$th training realization is simulated by randomly
selecting $L$ interfering signals from a DOA grid spanning the range of
$10^{0}$ to $170^{0}$, with a grid spacing of $1^{0}$. The interfering jammers
are allocated random powers uniformly distributed with INR from 10 dB to 20
dB. For this given scenario, the received correlation function, which includes
the desired source signal, is calculated corresponding to the full sensor
configuration. The corresponding optimum configuration is found through
enumeration and also through the proposed SBSA algorithm. The process is
repeated 30000 times against a given desired source DOA to complete the
training data set. Similarly, a small sized validation data set is generated
for hyperparamter tuning for model selection, minibatch size and learning
rate.

For the training stage, the weights of the neural network are optimized to
minimize the mean squared error between the label and the output of the
network. The ADAM algorithm is used to carry out parameter optimization. This
algorithm used an efficient implementation of mini-batch stochastic gradient
descent called the ADAM algorithm [130], that improves learning as a function
of the most recent gradient values. The learning rate is set to 0.001 and
dropout regularization with the keep probability of 0.9 is used. The initial
weights are initialized using the Xavier initialization.

The robustness of the learned models is demonstrated by generating the testing
data that is different from the training stage, and is generated by assuming
the DOA of the interfering signals off grid. This is simulated by adding the
Gaussian noise to the interference DOAs on the grid. We also present the
results under limited and unlimited data snapshots. Additionally, to ensure
the selection of $P$ antenna locations at the output of DNN, we declare the
$P$ highest values in the output as the selected sensor locations. Therefore,
generalization in this context means that the learned DNN works on different
interference settings which can change according to changing environment
conditions.

###  7.4 Simulations

In this section, we show the effectiveness of the proposed approach for sparse
array design achieving MaxSINR. The results are examined first by training the
DNN to learn the enumerated optimum array configurations. Then, we demonstrate
in the follow-on examples the effectiveness of the DNN when trained by the
labels drawn from SBSA algorithm.

####  7.4.1 Enumerated design

In this example, we pose the problem as selecting $P=6$ antennas from $N=12$
possible equally spaced locations with inter-element spacing of $\lambda/2$.
For all numerical results, we use a network with three hidden layers, one
input layer, and one output layer. Accordingly, the input to the network is of
size 23 and output is size 12.

Figure 7.11 shows the output SINR performance comparisons for different array
configurations. The horizontal axis is the DOA of the desired point source,
and the performance is computed at six different source DOAs varying from
$15^{0}$ to $90^{0}$ in steps of $15^{0}$. The SNR of the desired signal is
$0$ dB, and the interference-to-noise-ratio (INR) for each interference is
chosen randomly between $10-20$ dB for a given realization. The results
presented in Fig. 7.11 are obtained by using unlimited number of data
snapshots (USS) and employing enumerated labels to train the DNN. The network
performance is reported by averaging over $900$ testing scenarios for each
desired source DOA. It is evident that the DNN-EN design (prediction using the
enumerated labels) performs close (0.45 dB trade off) to the performance of
the optimum array found by enumeration ($980$ possible configurations). The
latter gives the highest SINR performance but involves expensive singular
value decomposition (SVD) for each enumeration and is also not scale-able with
the problem size, facing the curse of dimensionality.

The DNN-EN design performance is also compared with the NNC (nearest neighbour
correlation) design which returns the label corresponding to the input nearest
neighbour correlation function (in terms of mean square error). NNC design is
simply a lookup table such that for a given test data it returns the label of
the closest training example by sifting over the entire training set. It is
noted that the DNN-EN design outperforms the NNC design with an average
performance gain of 0.5 dB. The former approach not only offers superior
performance but also is more economically viable from ‘Edge computing’
perspective as compared to the nearest neighbour design which requires to
maintain a large dictionary and run exhaustive search over the entire training
set for each prediction. Similar results are obtained for the NNC design
obtained by minimizing the mean absolute error instead of the mean square
error. For the underlying case, the DNN has around 88 $\%$ accuracy for the
training data, and has around 54$\%$ accuracy on the test data (meaning that
all sensor locations are correctly predicted). This implies that the superior
performance of the DNN is not simply because it memorizes the optimum but due
to the ability of DNN to generalize the learning as applied to the test set
which is not seen before by the DNN. The superior performance of DNN over NNC
design again signifies that the DNN doesn’t memorize a lookup table
implementation to locate the nearest training data point and output the
corresponding memorized optimal sparse configuration. It is also clear that
the proposed design yields significant gains over a compact ULA, sparse ULA
and randomly selected array topology. The utility of the effective sparse
design is also evident from the worst case performance which exhibits
significantly low output SINR of around -5 dB on average.

Figure 7.11: Performance comparison enumerated design under unlimited
snapshots Figure 7.12: Performance comparison SBSA design under unlimited
snapshots

####  7.4.2 DNN based SBSA design

Although the training phase is entirely offline, it is infeasible to train the
DNN while relying on the enumerated results. This is because the number of
possible sparse solutions can be considerably larger even for the modest size
problem. For instance, choosing 14 elements out of 24 possible locations
results in the order of $10^{6}$ candidate spare configurations for each
training data, i.e., for each environment scenario. In order to circumvent
this problem and generate a large amount of training data labels, we resort to
the proposed SBSA design. Fig. 7.12 shows the performance of SBSA design which
is merely 0.3 dB down than the design obtained through enumeration. Quite
similar to the DNN-EN design, the DNN-SBSA design is around 0.4 dB down from
DNN-SBSA design. However, this places the DNN-SBSA design further 0.8 dB
suboptimal to the enumerated design. This is still a reasonable performance
yielding significant dividends over the commonly used compact ULA, sparse ULA
and random sparse topology.

####  7.4.3 Robust design

In order to gauge the robustness of the DNN based scheme, the performance is
evaluated under a limited number of data snapshots. Also, the desired source
DOA is perturbed with Gaussian noise of 0.250 variance, in generating the test
data to account of possible uncertainty around the desired source DOA. For
simulating the limited snapshot scenario, 1000 snapshots are generated
assuming the incoming signals (source and interfering signals) are independent
BPSK signals in the presence of Gaussian noise. The correlation matrix under
limited snapshots doesn’t follow the Toeplitz structure. Therefore, we average
along the diagonal entries of the correlation matrix to calculate the average
correlation values. Figs. 7.13 and 7.14 show the performance of DNN-EN and
DNN-SBSA designs under the limited data snapshots. It is clear from the
figures that the performance is largely preserved with an SINR discrepancy of
less than 0.01 dB demonstrating the robustness of the proposed scheme. The NNC
design, in this case is suboptimal with more than 0.3 dB additional
performance loss.

####  7.4.4 Performance comparisions with state-of-the-art

The performance of the proposed SBSA, DNN-EN and DNN-SBSA are compared with
existing work on sparse array design which is based on SDR and SCA approaches
[20, 71]. It is clear from Fig. 15 that the SBSA algorithm outperforms the
other designs and is also more than 100X computationally efficient as compared
to the SDR and SCA (Wang etal. [20]) approach. However, it is only fair to
compare the SBSA design with the SCA (Wang etal.) approach because both
incorporate the apriori knowledge of interference parameters. Therefore, in
comparing the data dependent designs, it is found that SDR design (also the
SDR-symmetric [71]) is comparable to the DNN-EN design, whereas the DNN-SBSA
is marginally suboptimal with the average performance degradation of 0.37 dB.
This slight performance trade off is readily leveraged by the real time
implementation of the DNN-SBSA algorithm implementing the Capon beamformer in
time frames of the order of few milli-seconds.

###  7.5 Conclusion

This chapter considered sparse array design for maximizing the beamformer
output SINR for a desired source in an interference active scenario. A sparse
beamformer spectral analysis (SBSA) algorithm was proposed which provided an
insightful perspective for MaxSINR beamformer design. Also a DNN based
approach was developed to configure a data driven sparse beamfomer by learning
the enumerated algorithm as well as SBSA design. The proposed methodology
combined the merits of the data dependent designs and the designs assuming
prior information of interference parameters. It was shown through the design
examples that the proposed schemes are robust in the case of limited data
snapshots and showed promising performance with reduced computational
complexity.

Figure 7.13: Performance comparison enumerated design under 1000 snapshots
Figure 7.14: Performance comparison SBSA design under 1000 snapshots Figure
7.15: Performance comparisons with the state-of-the-art

## Chapter 8CONCLUSIONS AND RECOMMENDATIONS

In this work, sparse array design was proposed to improve active and passive
sensing for radar and communication applications. Novel algorithms for sparse
beamfomer design were presented for both receive and transmit signal
platforms. The proposed approaches were analyzed analytically as well as with
the help of realistically simulated data.

The thesis addressed the problem that the optimization of the array
configuration requires full data correlation matrix which is not readily
available in practice. Two different design approaches were considered; one
assumed prefixed position of subset of sensors so as to provide full array
augmentation, referred to as the hybrid-design approach, whereas the other
used the Toeplitz estimation of the autocorrelation at the missing spatial
lags. We also considered maximizing the beamformer output SINR for the case of
wideband signal models. Two different implementations, namely the TDL and the
DFT implementation schemes, were presented for the wideband sparse array
design. The sparse array configuration optimized in the DFT domain and later
imported to the TDL implementation scheme was proposed, yielding promising
results. In all cases, the MaxSINR optimum sparse array delivered considerable
dividends over other suboptimal sparse arrays and compact uniform array.

A frequency domain sparse array algorithm was proposed which provided an
insightful perspective for MaxSINR beamformer design. Also a DNN based
approach was developed to configure a data driven sparse beamfomer by learning
various sparse design algorithms to enable real time implementation. We also
considered active sensing using sparse array configuration and sought optimum
transmit beamforming for radar applications. It was shown that the active
sparse array design can achieve desirable beampattern characteristics such as
directing a high proportion of the transmitted power towards the perspective
target locations with minimal cross-correlation among target signals.

Future work will use the findings of this thesis to address some of the
practical issues for sparse beamformer design such as evaluating the
performance on real data and incorporating physical constraints such as
antenna radiation pattern and mutual coupling effects [131, 132].

## References

  * [1] G. R. Lockwood, J. R. Talman, and S. S. Brunke, “Real-time 3-D ultrasound imaging using sparse synthetic aperture beamforming,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 45, no. 4, pp. 980–988, July 1998.
  * [2] O. Mehanna, N. D. Sidiropoulos, and G. B. Giannakis, “Joint multicast beamforming and antenna selection,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 10, pp. 2660–2674, May 2013.
  * [3] Y. He and K. P. Chong, “Sensor scheduling for target tracking in sensor networks,” in _2004 43rd IEEE Conference on Decision and Control (CDC) (IEEE Cat. No.04CH37601)_ , vol. 1, Dec 2004, pp. 743–748 Vol.1.
  * [4] W. V. Cappellen, S. J. Wijnholds, and J. D. Bregman, “Sparse antenna array configurations in large aperture synthesis radio telescopes,” in _2006 European Radar Conference_ , Sept 2006, pp. 76–79.
  * [5] S. Joshi and S. Boyd, “Sensor selection via convex optimization,” _IEEE Transactions on Signal Processing_ , vol. 57, no. 2, pp. 451–462, Feb 2009.
  * [6] H. Godrich, A. P. Petropulu, and H. V. Poor, “Sensor selection in distributed multiple-radar architectures for localization: A knapsack problem formulation,” _IEEE Transactions on Signal Processing_ , vol. 60, no. 1, pp. 247–260, Jan 2012.
  * [7] A. Moffet, “Minimum-redundancy linear arrays,” _IEEE Transactions on Antennas and Propagation_ , vol. 16, no. 2, pp. 172–175, March 1968.
  * [8] P. Pal and P. P. Vaidyanathan, “Nested arrays: A novel approach to array processing with enhanced degrees of freedom,” _IEEE Transactions on Signal Processing_ , vol. 58, no. 8, pp. 4167–4181, Aug. 2010.
  * [9] ——, “Coprime sampling and the music algorithm,” in _2011 Digital Signal Processing and Signal Processing Education Meeting (DSP/SPE)_ , Jan. 2011, pp. 289–294.
  * [10] S. Qin, Y. D. Zhang, and M. G. Amin, “Generalized coprime array configurations for direction-of-arrival estimation,” _IEEE Transactions on Signal Processing_ , vol. 63, no. 6, pp. 1377–1390, March 2015.
  * [11] X. Wang, E. Aboutanios, M. Trinkle, and M. G. Amin, “Reconfigurable adaptive array beamforming by antenna selection,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 9, pp. 2385–2396, May 2014.
  * [12] X. Wang, M. G. Amin, X. Wang, and X. Cao, “Sparse array quiescent beamformer design combining adaptive and deterministic constraints,” _IEEE Transactions on Antennas and Propagation_ , vol. PP, no. 99, pp. 1–1, 2017.
  * [13] N. D. Sidiropoulos, T. N. Davidson, and Z.-Q. Luo, “Transmit beamforming for physical-layer multicasting,” _IEEE Transactions on Signal Processing_ , vol. 54, no. 6, pp. 2239–2251, June 2006.
  * [14] V. Roy, S. P. Chepuri, and G. Leus, “Sparsity-enforcing sensor selection for DOA estimation,” in _2013 5th IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)_ , Dec. 2013, pp. 340–343.
  * [15] W. Roberts, L. Xu, J. Li, and P. Stoica, “Sparse antenna array design for MIMO active sensing applications,” _IEEE Transactions on Antennas and Propagation_ , vol. 59, no. 3, pp. 846–858, March 2011.
  * [16] G. G. Raleigh and V. K. Jones, “Multivariate modulation and coding for wireless communication,” _IEEE Journal on Selected Areas in Communications_ , vol. 17, no. 5, pp. 851–866, May 1999.
  * [17] H. Bolcskei, D. Gesbert, and A. J. Paulraj, “On the capacity of OFDM-based spatial multiplexing systems,” _IEEE Transactions on Communications_ , vol. 50, no. 2, pp. 225–234, Feb. 2002.
  * [18] P. Chavali and A. Nehorai, “Cognitive radar for target tracking in multipath scenarios,” in _2010 International Waveform Diversity and Design Conference_ , Aug. 2010, pp. 000 110–000 114.
  * [19] J. L. Krolik, J. Farrell, and A. Steinhardt, “Exploiting multipath propagation for GMTI in urban environments,” in _2006 IEEE Conference on Radar_ , April 2006, pp. 4 pp.–.
  * [20] X. Wang, M. Amin, and X. Cao, “Analysis and design of optimum sparse array configurations for adaptive beamforming,” _IEEE Transactions on Signal Processing_ , vol. PP, no. 99, pp. 1–1, 2017.
  * [21] R. B. Ertel, P. Cardieri, K. W. Sowerby, T. S. Rappaport, and J. H. Reed, “Overview of spatial channel models for antenna array communication systems,” _IEEE Personal Communications_ , vol. 5, no. 1, pp. 10–22, Feb. 1998.
  * [22] G. Raleigh, S. N. Diggavi, A. F. Naguib, and A. Paulraj, “Characterization of fast fading vector channels for multi-antenna communication systems,” in _Proceedings of 1994 28th Asilomar Conference on Signals, Systems and Computers_ , vol. 2, Oct. 1994, pp. 853–857 vol.2.
  * [23] A. Abdi and M. Kaveh, “A space-time correlation model for multielement antenna systems in mobile fading channels,” _IEEE Journal on Selected Areas in Communications_ , vol. 20, no. 3, pp. 550–560, April 2002.
  * [24] A. Abdi, J. A. Barger, and M. Kaveh, “A parametric model for the distribution of the angle of arrival and the associated correlation function and power spectrum at the mobile station,” _IEEE Transactions on Vehicular Technology_ , vol. 51, no. 3, pp. 425–434, May 2002.
  * [25] S. Shahbazpanahi, A. B. Gershman, Z.-Q. Luo, and K. M. Wong, “Robust adaptive beamforming for general-rank signal models,” _IEEE Transactions on Signal Processing_ , vol. 51, no. 9, pp. 2257–2269, Sept. 2003.
  * [26] R. M. Gray, “Toeplitz and circulant matrices: A review,” _Foundations and Trends® in Communications and Information Theory_ , vol. 2, no. 3, pp. 155–239, 2006.
  * [27] W. C. Y. Lee, _Mobile Communications Engineering: Theory and Applications_. New York, NY, USA: McGraw-Hill, Inc., 1997.
  * [28] X. Wang, M. G. Amin, and X. Cao, “Optimum adaptive beamformer design with controlled quiescent pattern by antenna selection,” in _2017 IEEE Radar Conference (RadarConf)_ , May 2017, pp. 0749–0754.
  * [29] S. Boyd and L. Vandenberghe, _Convex Optimization_. New York, NY, USA: Cambridge University Press, 2004.
  * [30] S. P. Chepuri and G. Leus, “Sparsity-promoting sensor selection for non-linear measurement models,” _IEEE Transactions on Signal Processing_ , vol. 63, no. 3, pp. 684–698, Feb 2015.
  * [31] E. BouDaher, Y. Jia, F. Ahmad, and M. G. Amin, “Multi-frequency co-prime arrays for high-resolution direction-of-arrival estimation,” _IEEE Transactions on Signal Processing_ , vol. 63, no. 14, pp. 3797–3808, July 2015.
  * [32] H. Unz, “Linear arrays with arbitrarily distributed elements,” _IRE Transactions on Antennas and Propagation_ , vol. 8, no. 2, pp. 222–223, March 1960\.
  * [33] R. Harrington, “Sidelobe reduction by nonuniform element spacing,” _IRE Transactions on Antennas and Propagation_ , vol. 9, no. 2, pp. 187–192, March 1961.
  * [34] A. Maffett, “Array factors with nonuniform spacing parameter,” _IRE Transactions on Antennas and Propagation_ , vol. 10, no. 2, pp. 131–136, March 1962.
  * [35] Y. Lo and S. Lee, “A study of space-tapered arrays,” _IEEE Transactions on Antennas and Propagation_ , vol. 14, no. 1, pp. 22–30, 1966.
  * [36] P. Jarske, T. Saramaki, S. K. Mitra, and Y. Neuvo, “On properties and design of nonuniformly spaced linear arrays (antennas),” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 36, no. 3, pp. 372–380, March 1988.
  * [37] R. M. Leahy and B. D. Jeffs, “On the design of maximally sparse beamforming arrays,” _IEEE Transactions on Antennas and Propagation_ , vol. 39, no. 8, pp. 1178–1187, Aug 1991.
  * [38] D. Caratelli and M. C. Viganó, “Analytical synthesis technique for linear uniform-amplitude sparse arrays,” _Radio Science_ , vol. 46, no. 4, 2011\.
  * [39] R. L. Haupt, “Thinned arrays using genetic algorithms,” _IEEE Transactions on Antennas and Propagation_ , vol. 42, no. 7, pp. 993–999, July 1994\.
  * [40] A. Trucco and V. Murino, “Stochastic optimization of linear sparse arrays,” _IEEE Journal of Oceanic Engineering_ , vol. 24, no. 3, pp. 291–299, July 1999.
  * [41] B. Fuchs, “Application of convex relaxation to array synthesis problems,” _IEEE Transactions on Antennas and Propagation_ , vol. 62, no. 2, pp. 634–640, Feb. 2014.
  * [42] S. Eng Nai, W. Ser, Z. Liang Yu, and H. Chen, “Beampattern synthesis for linear and planar arrays with antenna selection by convex optimization,” _Antennas and Propagation, IEEE Transactions on_ , vol. 58, pp. 3923 – 3930, 01 2011.
  * [43] G. Prisco and M. D’Urso, “Maximally sparse arrays via sequential convex optimizations,” _IEEE Antennas and Wireless Propagation Letters_ , vol. 11, pp. 192–195, 2012.
  * [44] M. B. Hawes and W. Liu, “Sparse array design for wideband beamforming with reduced complexity in tapped delay-lines,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 22, no. 8, pp. 1236–1247, Aug 2014\.
  * [45] A. Goldsmith, _Wireless Communications_. New York, NY, USA: Cambridge University Press, 2005.
  * [46] H. L. V. Trees, _Detection, Estimation, and Modulation Theory: Radar-Sonar Signal Processing and Gaussian Signals in Noise_. Melbourne, FL, USA: Krieger Publishing Co., Inc., 1992.
  * [47] J. Li, P. Stoica, and Z. Wang, “On robust Capon beamforming and diagonal loading,” _IEEE Transactions on Signal Processing_ , vol. 51, no. 7, pp. 1702–1715, July 2003.
  * [48] X. Wang and M. Amin, “Design of optimum sparse array for robust MVDR beamforming against DOA mismatch,” in _2017 IEEE 7th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)_ , Dec 2017, pp. 1–5.
  * [49] Y. I. Abramovich, D. A. Gray, A. Y. Gorokhov, and N. K. Spencer, “Positive-definite Toeplitz completion in DOA estimation for nonuniform linear antenna arrays. I. Fully augmentable arrays,” _IEEE Transactions on Signal Processing_ , vol. 46, no. 9, pp. 2458–2471, Sep 1998.
  * [50] Y. Abramovich, N. Spencer, and A. Gorokhov, “Positive-definite Toeplitz completion in DOA estimation for nonuniform linear antenna arrays. II. Partially augmentable arrays,” _Trans. Sig. Proc._ , vol. 47, no. 6, pp. 1502–1521, Jun. 1999.
  * [51] A. Bertrand and M. Moonen, “Efficient sensor subset selection and link failure response for linear MMSE signal estimation in wireless sensor networks,” in _2010 18th European Signal Processing Conference_ , Aug 2010, pp. 1092–1096.
  * [52] J. Szurley, A. Bertrand, M. Moonen, P. Ruckebusch, and I. Moerman, “Energy aware greedy subset selection for speech enhancement in wireless acoustic sensor networks,” in _2012 Proceedings of the 20th European Signal Processing Conference (EUSIPCO)_ , Aug 2012, pp. 789–793.
  * [53] J. Zhang, S. P. Chepuri, R. C. Hendriks, and R. Heusdens, “Microphone subset selection for MVDR beamformer based noise reduction,” _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 26, no. 3, pp. 550–563, March 2018.
  * [54] P. Stoica and R. L. Moses, _Introduction to spectral analysis_. Upper Saddle River, N.J. : Prentice Hall, 1997.
  * [55] D. L. Donoho, “For most large underdetermined systems of linear equations the minimal $l_{1}$-norm solution is also the sparsest solution,” _Communications on Pure and Applied Mathematics_ , vol. 59, no. 6, pp. 797–829, 2006.
  * [56] A. M. Bruckstein, D. L. Donoho, and M. Elad, “From sparse solutions of systems of equations to sparse modeling of signals and images,” _SIAM Rev._ , vol. 51, no. 1, pp. 34–81, Feb. 2009. [Online]. Available: http://dx.doi.org/10.1137/060657704
  * [57] A. Y. Yang, S. S. Sastry, A. Ganesh, and Y. Ma, “Fast $l_{1}$-minimization algorithms and an application in robust face recognition: A review,” in _2010 IEEE International Conference on Image Processing_ , Sep. 2010, pp. 1849–1852.
  * [58] E. J. Candès, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted $l_{1}$ minimization,” _Journal of Fourier Analysis and Applications_ , vol. 14, no. 5, pp. 877–905, Dec. 2008.
  * [59] M. Bengtsson and B. Ottersten, “Optimal downlink beamforming using semidefinite optimization,” 1999.
  * [60] S. A. Hamza, M. G. Amin, and G. Fabrizio, “Optimum sparse array beamforming for general rank signal models,” in _2018 IEEE Radar Conference (RadarConf18)_ , April 2018, pp. 1343–1347.
  * [61] Z. q. Luo, W. k. Ma, A. M. c. So, Y. Ye, and S. Zhang, “Semidefinite relaxation of quadratic optimization problems,” _IEEE Signal Processing Magazine_ , vol. 27, no. 3, pp. 20–34, May 2010.
  * [62] B. Recht, M. Fazel, and P. A. Parrilo, “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,” _SIAM Rev._ , vol. 52, no. 3, pp. 471–501, Aug. 2010.
  * [63] K. Mohan and M. Fazel, “Iterative reweighted algorithms for matrix rank minimization,” _J. Mach. Learn. Res._ , vol. 13, no. 1, pp. 3441–3473, Nov. 2012.
  * [64] R. Rajamäki and V. Koivunen, “Symmetric sparse linear array for active imaging,” in _2018 IEEE 10th Sensor Array and Multichannel Signal Processing Workshop (SAM)_ , July 2018, pp. 46–50.
  * [65] P. Pal and P. P. Vaidyanathan, “Nested arrays in two dimensions, Part 1: Geometrical considerations,” _IEEE Transactions on Signal Processing_ , vol. 60, no. 9, pp. 4694–4705, Sept 2012.
  * [66] K.-C. Huarng and C.-C. Teh, “Adaptive beamforming with conjugate symmetric weights,” _IEEE Transactions on Antennas and Propagation_ , vol. 39, no. 7, pp. 926–932, July 1991.
  * [67] A. Ahmed, Y. D. Zhang, and J. Zhang, “Coprime array design with minimum lag redundancy,” in _ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , May 2019, pp. 4125–4129.
  * [68] S. A. Hamza and M. G. Amin, “Two dimensional sparse array design for wideband signals ,” in _Big Data: Learning, Analytics, and Applications_ , vol. 10989, International Society for Optics and Photonics. SPIE, 2019, pp. 135 – 142.
  * [69] S. A. Hamza and M. G. Amin, “Sparse array dft beamformers for wideband sources,” in _2019 IEEE Radar Conference (RadarConf)_ , 2019, pp. 1–5.
  * [70] ——, “Optimum sparse array receive beamforming for wideband signal model,” in _2018 52nd Asilomar Conference on Signals, Systems, and Computers_ , Oct 2018, pp. 89–93.
  * [71] ——, “Hybrid sparse array beamforming design for general rank signal models,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 24, pp. 6215–6226, Dec 2019.
  * [72] ——, “Hybrid sparse array design for under-determined models,” in _ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , May 2019, pp. 4180–4184.
  * [73] ——, “Sparse array design utilizing matrix completion,” in _2019 Asilomar Conference on Signals, Systems, and Computers_ , Nov 2019.
  * [74] M. S. Ibrahim, A. Konar, M. Hong, and N. D. Sidiropoulos, “Mirror-prox sca algorithm for multicast beamforming and antenna selection,” in _2018 IEEE 19th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ , June 2018, pp. 1–5.
  * [75] C. Zhou, Y. Gu, Z. Shi, and Y. D. Zhang, “Off-grid direction-of-arrival estimation using coprime array interpolation,” _IEEE Signal Processing Letters_ , vol. 25, no. 11, pp. 1710–1714, Nov 2018\.
  * [76] C. Zhou, Y. Gu, X. Fan, Z. Shi, G. Mao, and Y. D. Zhang, “Direction-of-arrival estimation for coprime array via virtual array interpolation,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 22, pp. 5956–5971, Nov 2018.
  * [77] C. Liu, P. P. Vaidyanathan, and P. Pal, “Coprime coarray interpolation for DOA estimation via nuclear norm minimization,” in _2016 IEEE International Symposium on Circuits and Systems (ISCAS)_ , May 2016, pp. 2639–2642.
  * [78] S. M. Hosseini and M. A. Sebt, “Array interpolation using covariance matrix completion of minimum-size virtual array,” _IEEE Signal Processing Letters_ , vol. 24, no. 7, pp. 1063–1067, July 2017.
  * [79] H. Qiao and P. Pal, “Gridless line spectrum estimation and low-rank Toeplitz matrix compression using structured samplers: A regularization-free approach,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 9, pp. 2221–2236, May 2017.
  * [80] L. Yang and G. B. Giannakis, “Ultra-wideband communications: an idea whose time has come,” _IEEE Signal Processing Magazine_ , vol. 21, no. 6, pp. 26–54, Nov 2004.
  * [81] M. S. Brandstein and D. B. Ward, “Microphone arrays - signal processing techniques and applications,” in _Microphone Arrays_ , 2001.
  * [82] C. Paulson, J. Chang, C. Romero, J. Watson, F. Pearce, and N. Levin, “Ultra-wideband radar methods and techniques of medical sensing and imaging,” in _Proceedings of SPIE - The International Society for Optical Engineering_ , B. Cullum and J. Carter, Eds., vol. 6007, 2005.
  * [83] E. Pancera, “Medical applications of the ultra wideband technology,” in _2010 Loughborough Antennas Propagation Conference_ , Nov 2010, pp. 52–56.
  * [84] E. D. Di Claudio and R. Parisi, _Robust Wideband Beamforming_. John Wiley $\&$ Sons, Ltd, 2005, ch. 7, pp. 353–415.
  * [85] W. Liu and S. Weiss, _Wideband Beamforming: Concepts and Techniques_. John Wiley & Sons, 2010.
  * [86] K. Buckley, “Spatial/spectral filtering with linearly constrained minimum variance beamformers,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 35, no. 3, pp. 249–266, March 1987.
  * [87] B. D. V. Veen and K. M. Buckley, “Beamforming: a versatile approach to spatial filtering,” _IEEE ASSP Magazine_ , vol. 5, no. 2, pp. 4–24, April 1988.
  * [88] O. L. Frost, “An algorithm for linearly constrained adaptive array processing,” _Proceedings of the IEEE_ , vol. 60, no. 8, pp. 926–935, Aug 1972.
  * [89] M. Er and A. Cantoni, “Derivative constraints for broad-band element space antenna array processors,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 31, no. 6, pp. 1378–1393, Dec 1983.
  * [90] K. Buckley and L. Griffiths, “An adaptive generalized sidelobe canceller with derivative constraints,” _IEEE Transactions on Antennas and Propagation_ , vol. 34, no. 3, pp. 311–319, March 1986.
  * [91] L. C. Godara, “Application of the fast fourier transform to broadband beamforming,” _The Journal of the Acoustical Society of America_ , vol. 98, no. 1, pp. 230–240, 1995. [Online]. Available: https://doi.org/10.1121/1.413765
  * [92] J. Sheinvald and M. Wax, “Direction finding with fewer receivers via time-varying preprocessing,” _IEEE Transactions on Signal Processing_ , vol. 47, no. 1, pp. 2–9, Jan 1999.
  * [93] Moon-Sik Lee, V. Katkovnik, and Yong-Hoon Kim, “System modeling and signal processing for a switch antenna array radar,” _IEEE Transactions on Signal Processing_ , vol. 52, no. 6, pp. 1513–1523, June 2004.
  * [94] Li Yang, Liang Liwan, Pan Weifeng, Chen Yaqin, and Feng Zhenghe, “Signal processing method for switch antenna array of the fmcw radar,” in _Proceedings of the 2001 IEEE Radar Conference (Cat. No.01CH37200)_ , May 2001, pp. 289–293.
  * [95] Y. Asano, S. Ohshima, T. Harada, M. Ogawa, and K. Nishikawa, “Proposal of millimeter-wave holographic radar with antenna switching,” in _2001 IEEE MTT-S International Microwave Sympsoium Digest (Cat. No.01CH37157)_ , vol. 2, May 2001, pp. 1111–1114 vol.2.
  * [96] M. Guo, Y. D. Zhang, and T. Chen, “DOA estimation using compressed sparse array,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 15, pp. 4133–4146, Aug 2018.
  * [97] X. Wang, M. G. Amin, and X. Wang, “Optimum sparse array design for multiple beamformers with common receiver,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing_ , April 2018.
  * [98] M. H. Bae, I. H. Sohn, and S. B. Park, “Grating lobe reduction in ultrasonic synthetic focusing,” _Electronics Letters_ , vol. 27, no. 14, pp. 1225–1227, July 1991.
  * [99] J. Lu and J. Greenleaf, “Study of two-dimensional array transducers for limited diffraction beams,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 41, no. 5, pp. 724–739, 9 1994\.
  * [100] J. H. Doles and F. D. Benedict, “Broad-band array design using the asymptotic theory of unequally spaced arrays,” _IEEE Transactions on Antennas and Propagation_ , vol. 36, no. 1, pp. 27–33, Jan 1988.
  * [101] V. Murino, A. Trucco, and A. Tesei, “Beam pattern formulation and analysis for wide-band beamforming systems using sparse arrays,” _Signal Processing_ , vol. 56, no. 2, pp. 177 – 183, 1997.
  * [102] V. Murino, A. Trucco, and C. S. Regazzoni, “Synthesis of unequally spaced arrays by simulated annealing,” _IEEE Transactions on Signal Processing_ , vol. 44, no. 1, pp. 119–122, Jan 1996.
  * [103] D. B. Ward, R. A. Kennedy, and R. C. Williamson, “Theory and design of broadband sensor arrays with frequency invariant far-field beam patterns,” _The Journal of the Acoustical Society of America_ , vol. 97, no. 2, pp. 1023–1034, 1995.
  * [104] A. Trucco, “Synthesizing wide-band sparse arrays by simulated annealing,” in _MTS/IEEE Oceans 2001. An Ocean Odyssey. Conference Proceedings (IEEE Cat. No.01CH37295)_ , vol. 2, Nov 2001, pp. 989–994 vol.2.
  * [105] F. Anderson, W. Christensen, L. Fullerton, and B. Kortegaard, “Ultra-wideband beamforming in sparse arrays,” _IEE Proceedings H - Microwaves, Antennas and Propagation_ , vol. 138, no. 4, pp. 342–346, Aug 1991.
  * [106] I. S. Reed, J. D. Mallett, and L. E. Brennan, “Rapid convergence rate in adaptive arrays,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. AES-10, no. 6, pp. 853–863, Nov. 1974.
  * [107] S. A. Hamza and M. G. Amin, “Sparse array beamforming design for wideband signal models,” _IEEE Transactions on Aerospace and Electronic Systems_ , pp. 1–1, 2020.
  * [108] M. G. Amin, P. P. Vaidyanathan, Y. D. Zhang, and P. Pal, “Editorial for coprime special issue,” _Digital Signal Processing_ , vol. 61, no. Supplement C, pp. 1 – 2, 2017, special Issue on Coprime Sampling and Arrays.
  * [109] S. A. Hamza and M. G. Amin, “Sparse array design for maximizing the signal-to-interference-plus-noise-ratio by matrix completion,” _Digital Signal Processing_ , p. 102678, 2020. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1051200420300233
  * [110] E. J. Candès and B. Recht, “Exact matrix completion via convex optimization,” _Foundations of Computational Mathematics_ , vol. 9, no. 6, p. 717, Apr 2009.
  * [111] P. Stoica and A. Nehorai, “Music, maximum likelihood, and Cramer-Rao bound,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 37, no. 5, pp. 720–741, May 1989.
  * [112] S. A. Hamza and M. Amin, “Planar sparse array design for transmit beamforming (Conference Presentation),” in _Radar Sensor Technology XXIV_ , K. I. Ranney and A. M. Raynal, Eds., vol. 11408, International Society for Optics and Photonics. SPIE, 2020. [Online]. Available: https://doi.org/10.1117/12.2557923
  * [113] S. A. Hamza and M. G. Amin, “Sparse array design for transmit beamforming,” in _2020 IEEE International Radar Conference (RADAR)_ , 2020, pp. 560–565.
  * [114] A. Hassanien, M. G. Amin, Y. D. Zhang, and F. Ahmad, “Signaling strategies for dual-function radar communications: an overview,” _IEEE Aerospace and Electronic Systems Magazine_ , vol. 31, no. 10, pp. 36–45, 2016\.
  * [115] A. Hassanien, M. G. Amin, E. Aboutanios, and B. Himed, “Dual-function radar communication systems: A solution to the spectrum congestion problem,” _IEEE Signal Processing Magazine_ , vol. 36, no. 5, pp. 115–126, 2019.
  * [116] S. A. Hamza and M. G. Amin, “Sparse array receiver beamformer design for multi-functional antenna,” in _2020 28th European Signal Processing Conference (EUSIPCO)_ , 2021, pp. 1836–1840.
  * [117] P. Stoica, J. Li, and Y. Xie, “On probing signal design for MIMO radar,” _IEEE Transactions on Signal Processing_ , vol. 55, no. 8, pp. 4151–4161, Aug 2007.
  * [118] S. A. Hamza and M. G. Amin, “Learning sparse array capon beamformer design using deep learning approach,” in _2020 IEEE Radar Conference (RadarConf20)_ , 2020, pp. 1–5.
  * [119] S. A. Hamza and M. Amin, “A method for optimum sparse array beamforming design (Conference Presentation),” in _Big Data II: Learning, Analytics, and Applications_ , F. Ahmad, Ed., vol. 11395, International Society for Optics and Photonics. SPIE, 2020. [Online]. Available: https://doi.org/10.1117/12.2561526
  * [120] K. Li and J. Malik, “Learning to optimize,” in _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_ , 2017.
  * [121] M. Andrychowicz, M. Denil, S. G. Colmenarejo, M. W. Hoffman, D. Pfau, T. Schaul, B. Shillingford, and N. de Freitas, “Learning to learn by gradient descent by gradient descent,” in _Proceedings of the 30th International Conference on Neural Information Processing Systems_ , ser. NIPS’16. Red Hook, NY, USA: Curran Associates Inc., 2016, p. 3988–3996.
  * [122] T. J. O’Shea, T. C. Clancy, and R. W. McGwier, “Recurrent neural radio anomaly detection,” _CoRR_ , vol. abs/1611.00301, 2016.
  * [123] H. Sun, X. Chen, Q. Shi, M. Hong, X. Fu, and N. D. Sidiropoulos, “Learning to optimize: Training deep neural networks for interference management,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 20, pp. 5438–5453, 2018.
  * [124] Y. Shen, Y. Shi, J. Zhang, and K. B. Letaief, “LORM: learning to optimize for resource management in wireless networks with few training samples,” _IEEE Trans. Wireless Communications_ , vol. 19, no. 1, pp. 665–679, 2020\.
  * [125] P. Sprechmann, R. Litman, T. B. Yakar, A. M. Bronstein, and G. Sapiro, “Supervised sparse analysis and synthesis operators,” in _Advances in Neural Information Processing Systems 26: 27th Annual Conference on Neural Information Processing Systems 2013, Lake Tahoe, Nevada, United States_ , 2013, pp. 908–916.
  * [126] T. J. O’Shea, T. Erpek, and T. C. Clancy, “Deep learning based MIMO communications,” _CoRR_ , vol. abs/1707.07980, 2017.
  * [127] J. R. Hershey, J. L. Roux, and F. Weninger, “Deep unfolding: Model-based inspiration of novel deep architectures,” _CoRR_ , vol. abs/1409.2574, 2014\.
  * [128] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” _SIAM Journal on Imaging Sciences_ , vol. 2, no. 1, pp. 183–202, 2009.
  * [129] K. Gregor and Y. LeCun, “Learning fast approximations of sparse coding,” in _Proceedings of the 27th International Conference on International Conference on Machine Learning_ , ser. ICML’10. Madison, WI, USA: Omnipress, 2010, p. 399–406.
  * [130] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , Y. Bengio and Y. LeCun, Eds., 2015.
  * [131] E. BouDaher and A. Hoorfar, “Electromagnetic optimization using mixed-parameter and multiobjective covariance matrix adaptation evolution strategy,” _IEEE Transactions on Antennas and Propagation_ , vol. 63, no. 4, pp. 1712–1724, 2015.
  * [132] E. BouDaher, F. Ahmad, M. G. Amin, and A. Hoorfar, “Mutual coupling effect and compensation in non-uniform arrays for direction-of-arrival estimation,” _Digital Signal Processing_ , vol. 61, pp. 3 – 14, 2017, special Issue on Coprime Sampling and Arrays.


# Implicit biases in transit models using stellar pseudo-density

Gregory J. Gilbert  Department of Physics and Astronomy, University of
California, Los Angeles  Mason G. MacDougall  Department of Physics and
Astronomy, University of California, Los Angeles  Erik A. Petigura  Department
of Physics and Astronomy, University of California, Los Angeles

###### Abstract

The transit technique is responsible for the majority of exoplanet discoveries
to date. Characterizing these planets involves careful modeling of their
transit profiles. A common technique involves expressing the transit duration
using a density-like parameter, $\tilde{\rho}$, often called the “circular
density.” Most notably, the Kepler project – the largest analysis of transit
lightcurves to date – adopted a linear prior on $\tilde{\rho}$. Here, we show
that such a prior biases measurements of impact parameter, $b$, due to the
non-linear relationship between $\tilde{\rho}$ and transit duration. This bias
slightly favors low values ($b\lesssim 0.3$) and strongly disfavors high
values ($b\gtrsim 0.7$) unless transit signal-to-noise ratio is sufficient to
provide an independent constraint on $b$, a criterion that is not satisfied
for the majority of Kepler planets. Planet-to-star radius ratio, $r$, is also
biased due to $r{-}b$ covariance. Consequently, the median Kepler DR25 target
suffers a $1.6\%$ systematic underestimate of $r$. We present a techniques for
correcting these biases and for avoiding them in the first place.

††facilities: Kepler††software: astropy (Astropy Collaboration et al., 2018),
exoplanet (Foreman-Mackey et al., 2021), numpy (Harris et al., 2020), PyMC
(Salvatier et al., 2016), scipy (Virtanen et al., 2020), starry (Luger et al.,
2019)

##  1 Introduction

In the two decades since the discovery of the first transiting hot Jupiter
(Charbonneau et al., 2000; Henry et al., 2000), the transit technique has
grown to be the most prolific exoplanet detection method to date, accounting
for 77% of the current census. Contemporary work continues to rely heavily on
the transit technique. To wit, several transit-focused NASA and ESA missions
are either already on-sky (TESS, Ricker et al. 2015) or slated for launch in
the near future (PLATO, Rauer et al. 2014), and next-generation radial
velocity spectrographs have been designed for follow-up characterization of
known transiting planets (e.g. KPF, Gibson et al. 2016; MAROON-X, Seifahrt et
al. 2018). The transit technique will remain indispensable for exoplanet
astronomy for decades to come.

Accurate modeling of the transit lightcurve is a critical step for
characterizing transiting planets. At the most basic level, transit modeling
involves computing the time-dependent flux $F(t)$ of a star obscured by a
transiting planet relative to the unobscured flux $F_{0}$. If one assumes a
spherical planet and star, this computation depends strictly on the planet-to-
star size ratio $r$, the (time-dependent) center-to-center sky-projected
planet-to-star separation $z$ (measured in units of $R_{\star}$), and the
radial dependence of the stellar limb-darkening profile $\\{u\\}$. Early
analyses computed $F(z;r,\\{u\\})$ via numerical integration, but today the
most widely used method is the Mandel & Agol (2002) model, which expresses the
transit lightcurve via an analytic solution to $F(z;r,\\{u\\})$ for several
limb darkening profiles which can be described by a small set of limb-
darkening parameters.

In order to model time-series photometry, one must convert $F(z;r,\\{u\\})$
into $F(t;r,\\{u\\})$. While $z$ is the only parameter that varies with time,
one may choose how to specify the function that maps $t\rightarrow z$. If one
assumes strict periodicity of transits and a constant projected velocity
during transit., then in the limit $r\rightarrow 0$, $z(t)$ may be specified
completely by an orbital period, $P$, a transit mid-point, $t_{0}$, an impact
parameter, $b$, and 1st-to-4th contact transit duration, $T_{14}$.111Several
alternative transit durations may be substituted for $T_{14}$: (1) the 2nd-
to-3rd contact duration, $T_{23}$, (2) the center-to-center contact duration,
$T_{\rm cc}$, also called the 1.5-to-3.5 contact duration, or (3) the full-
width-half-max duration, $T_{\rm FWHM}$, which may be defined in relation to
the transit depth. While each has its merits (see Kipping 2010a for
discussion), we adopt $T_{14}$ throughout this work because it is the transit
duration which is most readily defined for all grazing and non-grazing transit
geometries. This parameterization — $F(t;P,t_{0},r,b,T_{14})$ — is convenient
and is closely linked to the transit geometry.

An alternative approach is to specify $T_{14}$ from a combination of scaled
separation $a/R_{\star}$, orbital eccentricity $e$, argument of periastron
$\omega$, and projected inclination $\cos i$, following Winn (2010) as

$T_{14}\simeq\frac{P}{\pi}\sin^{-1}\Big{(}\frac{R_{\star}}{a}\frac{\sqrt{(1+r)^{2}-b^{2}}}{\sin
i}\Big{)}\Big{(}\frac{\sqrt{1-e^{2}}}{1+e\sin\omega}\Big{)}$ (1)

$b=\frac{a\cos i}{R_{\star}}\Big{(}\frac{1-e^{2}}{1+e\sin\omega}\Big{)}$ (2)

Now the lightcurve is specified by the function
$F(t;P,t_{0},r,a/R_{\star},b,e,\omega)$, which is similar to the
parameterization used by the EXOFAST suite (Eastman et al., 2013; Eastman,
2017).222In practice, EXOFAST uses $\log(a/R_{\star})$ and expresses $b$ as
$\cos i$; $\\{e,\omega\\}$ is usually specified as
$\\{\sqrt{e}\sin\omega,\sqrt{e}\cos\omega\\}$ in order to establish uniform
priors on $e$ and $\omega$ and to a avoid a boundary issue at $e=0$. A related
approach is to replace $a/R_{\star}$ with stellar density by employing
Kepler’s third law. Thus, the light curve may also be parameterized by
$F(t;P,t_{0},r,\rho_{\star},\cos i,e,\omega)$.

These two eccentricity-explicit parameterizations have the advantage that the
lightcurve has been specified completely by properties of the star, planet,
and planetary orbit; the disadvantage is that five parameters have been
replaced by seven, and thus significant degeneracies between
$\\{a/R_{\star},e,\omega\\}$ or $\\{\rho_{\star},e,\omega\\}$ are inevitable.
These degeneracies lead to inefficiencies with light curve fitting and
posterior sampling.

A common shortcut is to fit the lightcurve assuming that $e=0$ even though the
orbit may, in fact, be eccentric. This assumption reduces the number of free
parameters back to five, but $\rho_{\star}$ can no longer be thought of as a
stellar density. Rather, it is a stand-in for duration which merely has units
of density, defined by Seager & Mallén-Ornelas (2003) as

$\tilde{\rho}\equiv\Bigg{(}\frac{4\pi^{2}}{P^{2}G}\Bigg{)}\Bigg{(}\frac{(1+r)^{2}-b^{2}\big{(}1-\sin^{2}[\pi
T_{14}/P]\big{)}}{\sin^{2}[\pi T_{14}/P]}\Bigg{)}^{3/2}$ (3)

where $G$ is Newton’s gravitational constant. This quantity $\tilde{\rho}$ is
sometimes referred to as the “mean stellar density,” the “circular density,”
or the “observed density,” but we prefer to call it the “pseudo-density”
because (1) the other names are confusing, and (2) $\tilde{\rho}$ matches the
true stellar density only when numerous assumptions are met (see Kipping,
2014).

Because the prior expectation for $\tilde{\rho}$ is a complicated function of
$\rho_{\star}$, $b$, $e$, and $\omega$, naïvely placing a flat prior on
$\tilde{\rho}$ and adopting it as a fitting parameter induces undesired biases
on $T_{14}$ and $b$.

To date, $\tilde{\rho}$ has enjoyed widespread use in the exoplanet
literature. For example, the Kepler project (Borucki et al. 2010; the largest
analysis of transit lightcurves to date) fit their lightcurves with the
$F(t;P,t_{0},r,\tilde{\rho},b)$ parameterization (Rowe et al., 2014, 2015;
Mullally et al., 2015; Coughlin et al., 2016; Thompson et al., 2018). We
discuss the effects of that choice in §5. More broadly, this paper
investigates the implicit biases on impact parameter and other light curve
parameters that result from the use of $\tilde{\rho}$.

Throughout this work, we assume that all transit signals under investigation
have been thoroughly vetted such that the detected signal is known to be a
real transit at high confidence. The methods employed in this work are thus
appropriate for parameter estimation but not for transit detection or vetting.

This paper is organized as follows. In §2 we empirically demonstrate the
origin of the $\tilde{\rho}$ bias by fitting a transit lightcurve model to
simulated photometry using the Kepler project parameterization; we then
demonstrate that our preferred parameterization does not suffer from this
bias. In §3 we present a numerical experiment which isolates the effects of
various model assumptions on posterior inferences. In §4 we analytically
derive the Jacobian of the coordinate transformation
$T_{14}\rightarrow\tilde{\rho}$ which explains the origin of the empirical
bias. In §5 we show that the $\tilde{\rho}$ bias has affected most posterior
inferences of $b$ and $r$ derived from Kepler data. In §6 we summarize our
conclusions and discuss other biases which arise from using related
parameterizations such as $a/R_{\star}$.

##  2 Understanding parameter biases with fits to synthetic photometry

To illustrate the $\tilde{\rho}$ bias, we simulated photometric observations
of a warm mini-Neptune ($P=15$ days, $r_{p}=3.3\ R_{\oplus}$) on a circular
orbit around a Sun-like star, transiting at impact parameter $b=0.5$. We
simulated data with a 30 minute observing cadence (matching Kepler’s long
cadence observing mode) within $\pm T$ from the transit center. All
photometric data were oversampled by a factor of 7 and integrated using
Simpson’s rule to account for the effects of finite integration time (Kipping,
2010b). The white noise level was tuned to produce S/N = 16, which is slightly
lower than the median Kepler value and results in a posterior model with
$\sigma_{r}/r\approx 0.10$ and $\sigma_{T}/T\approx 0.05$, where
$\sigma_{r}/r$ corresponds to the fractional posterior measurement, and
similar for $T$. We chose these values in order to produce a transit which is
similar to those found by Kepler. Ground-truth simulation parameters are
listed in Table 1, and simulated photometry is shown in Figure 1.

Table 1: Ground-truth simulation parameters; simulated photometry is shown in Figure 1. Parameter | Value
---|---
$M_{\star}[M_{\odot}]$ | $1.0$
$R_{\star}[R_{\odot}]$ | $1.0$
$u_{1},u_{2}$ | $0.40,0.25$
$P$ [days] | $15.0$
$r$ | $0.03$
$b$ | $0.5$
$T_{14}$ [hrs] | $3.29$
$S/N$ | $16$
Figure 1: Simulated photometry for a mini-Neptune on a circular 15 day orbit
around a Sun-like star, transiting at $b=0.5$. The orange line indicates the
ground truth transit model. Grey points show simulated observations at a one
minute observing cadence; black circles are binned to 30 minutes. The white
noise level was set to S/N = 16, close to the Kepler median. See Table 1 for
ground-truth simulation parameters.

The transit model was specified using a standard pseudo-density
parameterization: $\\{P,t_{0},r,b,\tilde{\rho}\\}$. In order to minimize
confounding factors, we held $P$ and $t_{0}$ fixed at their injected values;
we also held the mean out-of-transit flux, $F_{0}$, and photometric white
noise level, $\sigma^{2}_{\rm phot}$, fixed to their true values, which is
equivalent to assuming the raw photometry has been accurately pre-whitened.
For the remaining transit parameters, we adopted broad weakly informative
priors with permissive bounds (see Table 2 for details), for a total of three
free parameters per model: $\\{r,b,\tilde{\rho}\\}$ (the “$\tilde{\rho}$
basis”), or $\\{\log r,b,\log T_{14}\\}$ (the “$\log T$ basis”). We chose the
later basis because $T_{14}$ is typically well constrained by the data and
furthermore may be assigned priors in a sensible fashion; sampling in $\log r$
and $\log T_{14}$ is equivalent to placing log-uniform priors on $r$ and
$T_{14}$ which facilitates exploration of posterior values over different
orders of magnitude. We modeled a circular transit in all cases and held
stellar mass, radius, and limb darkening to their true values during the fit;
there is no loss of generality in this approach, because as long as we ignore
minuscule ingress/egress asymmetry that exists for eccentric transits (Barnes,
2007), there is no difference between a circular and eccentric transit. In
order to avoid complications which arise when modeling grazing transits, we
restricted impact parameters to $b<1-r$.333A common approach (which we did not
adopt) is to draw samples uniformly from the $r-b$ plane using triangular
sampling (Espinoza, 2018). However, naive application of this method induces a
marginal prior on $r$, so caution must be taken to ensure that priors are
established as intended. To confirm that this restriction is permissible, we
explored the parameter space near the limb of the star following the
methodology of Gilbert (2022) and verified that the simulated transit is
inconsistent with a grazing geometry.

Table 2: Priors on model parameters for simulated lightcurve. Parameter | Value | Parameter | Value
---|---|---|---
$r$ | $\mathcal{U}(0.01,0.1)$ | $\log r$ | $\mathcal{U}(-2,-1)$
$b$ | $\mathcal{U}(0,1-r)$ | $b$ | $\mathcal{U}(0,1-r)$
$\tilde{\rho}/\rho_{\odot}$ | $\mathcal{U}(0.1,10)$ | $\log[T_{14}/\rm hr]$ | $\mathcal{U}(1,10)$

We drew samples from the posterior using Hamiltonian Monte Carlo (HMC; Neal,
2011) and the No U-Turn Sampler (NUTS; Hoffman & Gelman, 2011). Each model
iteration consisted of two chains run for 5000 tuning steps and 20,000 draws,
producing an effective number of samples greater than 11,000 for all
parameters for each of the the two parameterizations.

Posterior corner plots for the quantities of interest are shown in Figure 2.
The most notable difference is in the 1D marginalized distribution of impact
parameter. When sampling using the $\tilde{\rho}$ basis, the posterior is
biased toward low $b$; as a point of reference, 74% of the probability mass is
below $b=0.5$, the injected value. When sampling using the $\log T$ basis,
however, the distribution of $b$ is nearly uniform over the allowed range,
reflecting the fact that for a low signal-to-noise transit the impact
parameter is largely unconstrained. Our results did not substantially change
when simulating a one minute observing cadence (matching Kepler’s short
cadence mode), indicating that the $\tilde{\rho}$ bias arises from the model
parameterization and is not an artifact of data binning. We also repeated the
analysis using $r=0.1$ and $r=0.01$ and found that the results did not change.

Clearly, the results are inconsistent between models – which contain identical
underlying physics and differ only in their parameter bases – so at least one
of the two models has produced biased inference. In the sections that follow,
we present both a numerical argument (§3) and an analytic argument (§4) which
demonstrate that the $\log T$ basis has produced the desired result.

Figure 2: Posterior corner plots when sampling in the $\tilde{\rho}$ basis
(orange) vs the $\log T_{14}$ basis (purple). The bias on impact parameter,
$b$, is apparent when sampling with $\tilde{\rho}$ but is resolved when
sampling in $T_{14}$.

##  3 Numerical sampling experiment

We will now demonstrate that the bias on $b$ seen in the previous section
arises solely from the model parameterization and not from vagaries of the
MCMC sampling algorithm or peculiarities of the noise realization in the
photometry.

To do so, we performed a numerical experiment which approximated the
lightcurve modeling procedure from §2 by drawing samples directly from the
prior distributions and then applying an a posteriori importance weighting
designed to mimic the constraints imposed by the photometry. When determining
these importance weights, we employed a Gaussian likelihood function and
approximated the (covariant) parameter constraints from §2 as independent
univariate Gaussians. The key advantage of this method is that we no longer
needed to directly fit the photometry, thereby eliminating potential
confounding factors introduced by the photometry and the sampler.

###  3.1 Experimental setup

We adopted the same fiducial star-planet system as §2, placing a
$3.3R_{\oplus}$ mini-Neptune on a circular 15 day orbit around a solar twin.
We fixed the ephemeris $\\{P,t_{0}\\}$ throughout and placed uniform interval
priors on all other parameters $\\{r,b,\tilde{\rho},\log T_{14}\\}$ as before
(see Table 2), with the small modification that we now allow $b$ to range over
all detectable values, i.e. $b\sim\mathcal{U}(0,1+r)$; this modification is
acceptable because our sampling procedure (see below) avoids the usual issues
which arise when fitting grazing transits (see Gilbert, 2022).

For the first iteration of the experiment we adopted the $\tilde{\rho}$ basis
$\\{r,b,\tilde{\rho}\\}$ and drew random samples directly from the prior
distributions. We next calculated transit duration using

$T_{14}=\frac{PR_{\star}}{\pi a}\Big{(}(1+r)^{2}-b^{2}\Big{)}^{1/2}$ (4)

for each sample. Here, we have approximated Equation 1 by using the small
angle approximation $\sin^{-1}\phi\approx\phi$ and
$i\approx\pi/2\rightarrow\sin i\approx 1$. The scaled separation can be
calculated from Kepler’s Third Law as
$a/R_{\star}=[(GP^{2}\rho)/(3\pi)]^{1/3}$.

For subsequent iterations of the experiment, we modified the procedure to use
use three alternative parameter bases: (1) $\\{\log
r,b,\log{\tilde{\rho}}\\}$, (2) $\\{r,b,T_{14}\\}$, and (3) $\\{\log r,b,\log
T_{14}\\}$. We chose these parameterizations in order to explore the effects
of uniform vs log-uniform priors in addition to the effect of substituting
$\tilde{\rho}\rightarrow T_{14}$. We followed the same sampling procedure as
before, except when drawing samples of $T_{14}$ or $\log T_{14}$ we calculated
$\tilde{\rho}$ following Equation 3.

Mimicking the simulated light curve in §2, we assumed that we could constrain
$r$ to $10\%$ accuracy and $T_{14}$ to $5\%$ accuracy, with independent
Gaussian precision from the photometry (i.e. $\sigma_{r}/r=0.1$,
$\sigma_{T}/T=0.05$). We further assumed that the impact parameter would be
entirely unconstrained by the data. These uncertainties are representative of
typical values, but we have removed the covariance and forced them to be
Gaussian (or unconstrained), which eases interpretation.

We imposed our assumed measurement uncertainties on $r$ and $T_{14}$ by
calculating the log-likelihood of each $i^{\rm th}$ sample

$\log\mathcal{L}_{i}=-\frac{1}{2}\Big{(}\frac{T_{i}-T_{\rm
true}}{\sigma_{T}}\Big{)}^{2}-\frac{1}{2}\Big{(}\frac{r_{i}-r_{\rm
true}}{\sigma_{r}}\Big{)}^{2}$ (5)

which assumes a Gaussian likelihood function. We then weighted each sample by

$w_{i}=\frac{\mathcal{L}_{i}}{\sum_{i}\mathcal{L}_{i}}$ (6)

to produce our synthetic posterior distributions.

###  3.2 Bias on impact parameter

The results of our numerical experiment are summarized in Figure 3. As
expected, when parameterizing the model as $\\{r,b,\tilde{\rho}\\}$ with
uniform priors, we obtain biased results that are qualitatively similar to
those produced in §2 (i.e. by fitting the photometry directly). Notably,
sampling in $\tilde{\rho}$ produces a strong prior on $T_{14}$ (purple) which
is not physically motivated. Because $T_{14}$ is constrained to 5%, the data
overwhelm the prior and the $T_{14}$ posterior distribution (orange) is only
slightly biased. The posterior on impact parameter, however, is clearly
different from the prior even though our model included no information about
impact parameter. Because we have (by construction) placed no measurement
constraint on $b$, the posterior distribution should match the prior. In
reality however, the posterior is tilted toward $b=0$, giving the illusion of
a (modestly) constrained posterior.

The $\tilde{\rho}{-}b$ bias is resolved by using any of the alternative
parameterizations which substitute $\log\tilde{\rho}$, $T_{14}$, or $\log
T_{14}$ for $\tilde{\rho}$. Although using the substitution
$\tilde{\rho}\rightarrow\log\tilde{\rho}$ may seem at first glace to be the
simplest choice (requiring little change from existing practices), we argue
that using either of the duration-based parameterizations is preferable for
two reasons. First, the results are insensitive to the exact choice of
(reasonable) prior placed on $T_{14}$, whereas they are highly sensitive to
the prior placed on $\tilde{\rho}$; insensitivity to priors is in general a
desirable feature of robust inference. Second, setting prior interval bounds
on $\tilde{\rho}$ is a non-intuitive task, requiring careful consideration of
the true stellar density and orbital elements. In contrast, principled priors
may be placed on the transit duration quite simply following inspection of the
transit lightcurve. In fact, setting bounds on $T_{14}$ is so straightforward
that it could even be done algorithmically following the output of a box-least
squares transit search (Kovács et al., 2002). The bottom line is that given
the choice between options which produce equivalent results, we prefer the
simpler of the two.

Figure 3: Results of the numerical sampling experiment described in §3. Each
row corresponds to the prior, likelihood, and posterior for a given model
parameterization. For visual clarity, the height of the $T_{14}$ posterior has
been reduced by a factor of 3 on all plots. The difference in the prior
distribution on $b$ for rows 1 & 3 compared to rows 2 & 4 stems from the use
of $r$ vs $\log r$, respectively. Sampling with a uniform prior on
$\tilde{\rho}$ (top row) produces a nonuniform prior on $T_{14}$ and a biased
posterior for $b$. In contrast, sampling in any of the other parameter bases
produces a posterior estimate of $b$ which matches the prior, except in cases
where constraints on $r$ would produce a non-transiting orbit.

In summary, because we have decoupled the posteriors from complicating factors
(e.g. parameter covariances, sampler inefficiencies, etc.), we conclude that
the differences between posterior distributions obtained under the
$\tilde{\rho}$ basis versus the $\log T_{14}$ basis arise solely due the
parameterization. Furthermore, we conclude that the $\tilde{\rho}$ basis (with
a uniform prior) induces a bias on $b$, whereas the other options we have
presented produce unbiased estimates.

##  4 Mathematical origin of the bias

In the previous sections, we illustrated the biases on $b$ that result from
uniform and log-uniform priors on $\tilde{\rho}$ by exploring synthetic
photometry fits and simple numerical experiments. In this section, we
investigate the mathematical origins of this bias.

The transit parameter covariance matrix was previously derived by Carter et
al. (2008), but where their treatment prioritized analytic interpretability
(with a small sacrifice to accuracy), our treatment prioritizes accuracy (with
a small sacrifice to interpretability). Most importantly, the covariance
matrix derived by Carter et al. (2008) are least accurate as $b\rightarrow 1$
and in the presence of non-neglible limb darkening, which are precisely the
conditions under which the $\tilde{\rho}$ bias we are investigating become
most important. Thus, our work complements rather than supplants Carter et al.
(2008).

When modeling light curves, our main goal is to derive the posterior
probability density function, $p(\vec{x})$, i.e. the probability that a set of
planet properties $\vec{x}$ resides in an infinitesimal volume element
spanning $\vec{x}$ to $\vec{x}+d\vec{x}$. However, this probability is not
invariant under changes in parameterization. Specifically, for our problem,
$p(T_{14})/dT_{14}\neq p(\tilde{\rho})/d\tilde{\rho}$. To convert $p(T_{14})$
to $p(\tilde{\rho})$, one must account for the change in infinitesimal volume
element resulting from the $T_{14}\rightarrow\tilde{\rho}$ transformation,
i.e. the Jacobian

$J=\frac{d\tilde{\rho}}{dT_{14}}=-\frac{12\pi^{3}}{P^{3}G}\Bigg{(}(1+r)^{2}-b^{2}\Bigg{)}^{3/2}\Bigg{(}\frac{\pi
T_{14}}{P}\Bigg{)}^{-4}$ (7)

which we derive in the Appendix. The Jacobian of the transformation
$T_{14}\rightarrow\log\tilde{\rho}$ is simply

$J^{\prime}=\frac{d\log\tilde{\rho}}{dT_{14}}=-\frac{3}{T_{14}}$ (8)

which is independent of $b$, explaining why using $\log\tilde{\rho}$ in place
of $\tilde{\rho}$ produces unbiased posteriors.

In Figure 4, we show the analytic Jacobian in Equation 7 alongside the
simulated posterior samples of $b$ obtained in §2 and the numerical results
obtained in §3. It is evident from inspection that the distributions are in
close agreement. We conclude that the non-uniform distribution of $b$ arises
from the combination of parameterization and (incorrect) prior, rather than
from any real constraint imposed by the data.

Figure 4: Posterior samples of $b$ from the simulated transit fit (orange
histogram, §2) and the numerical experiment (purple histogram, §3) are nearly
perfectly matched by the expected bias from the analytically derived Jacobian
(grey shaded region, §4).

##  5 Biased Kepler planet properties

We have shown that adopting a linear $\tilde{\rho}$ prior results in a biased
impact parameter. The Kepler project (Borucki et al., 2010; Rowe et al., 2014,
2015; Mullally et al., 2015; Coughlin et al., 2016; Thompson et al., 2018)
used such a parameterization (Jason Rowe, private communication). Therefore,
we expect biased $b$ in all cases except those where $b$ is strongly
constrained by the light curve itself. Because most Kepler planet candidates
exhibit modest transit signal-to-noise (median S/N = 22.4), the characteristic
“hill” shape we have seen for biased posterior $b$ distributions in the
previous three sections is also present in the posterior distributions of
nearly every Kepler planet candidate from DR25 (Thompson et al., 2018). Figure
5 illustrates the presence of the $b$ bias over a grid of orbital periods and
radii. Only the largest (and therefore highest S/N) planets consistently
exhibit meaningful constraints on $b$.

Figure 5: Posterior distributions of impact parameter for a random selection
of KOIs, organized in logarithmic bins on a $P-r$ grid. Data shown are the
posterior MCMC chains from Kepler Data Release 25 (Thompson et al., 2018),
described in detail in Rowe et al. (2014) and downloaded from the NASA
Exoplanet Archive (Akeson et al., 2013). Each posterior distribution is
plotted with 20% opacity so that dark regions indicate where many
distributions overlap; colors correspond to the median $b$ value for a given
KOI. For visual clarity, a maximum of twelve KOIs are plotted per panel. The
horizontal axis of each panel ranges over $b=(0,1.2)$; the vertical range of
each row is different, but the dashed line indicates the same distribution
height. The median S/N in each 2D bin is printed in the upper right-hand
corner of the grid squares. It is clear from inspection that most of the
objects (excluding the largest, highest S/N objects) show qualitatively
similar posterior distributions of $b$. The similarity is particularly
striking for small (low S/N) objects.

Due to signal-to-noise bias which disfavors the detection of high-$b$ transits
(Kipping & Sandford, 2016), the prior expectation on impact parameter is not
exactly flat, and so the posteriors exhibited in Kepler data will not exactly
match the idealized distribution we derived in §2-4. However, most Kepler
detections have $S/N>10$ and fall in the flat part of the detection
completeness curve (Christiansen et al., 2020). Thus, the appropriate prior
for the vast majority of Kepler planets should be nearly flat in $b$, with a
fall off at the value of $b$ that reduces $S/N$ to $\sim\ 10$.

Detection biases notwithstanding, the $\tilde{\rho}$ bias is easily understood
and corrected. Because the relationship between $\tilde{\rho}$, $b$, $r$, and
$T_{14}$ is known analytically (Seager & Mallén-Ornelas, 2003), one needs only
to apply the appropriate Jacobian weighting in order to transform an
unintended prior on $\tilde{\rho}$ into the desired prior on $b$ or $T_{14}$
(or any other basis parameter derivable from these quantities). Unbiased
parameter estimates can then be recovered from existing (biased) posterior
chains by implementing an importance sampling scheme which accounts for this
coordinate transformation, provided the chains are not too sparsely sampled in
their low probability regions. Specifically, one can sample from a
distribution $p_{1}(\vec{x})$ by reweighting samples from a different
distribution $p_{2}(\vec{x})$. An example of this reweighting scheme as
applied to a selection of DR25 targets is shown in Figure 6. A caveat is there
is increased sampling error since $p_{2}(\vec{x})$ is a different distribution
and the samples are not optimally distributed in $p_{1}(\vec{x})$. In essence
there are smaller number of “effective samples” after reweighting. Care must
therefore be taken to ensure that Jacobian-corrected posteriors are reliable,
and the reweighting scheme we have outlined here should not be applied
blindly.

Figure 6: Posterior samples of $b$ from a representative selection of DR25
targets before and after reweighting by the Jacobian to correct for biases
induced by sampling in $\tilde{\rho}$. All targets have $0.02<r<0.04$ and
$10<P<30$ days. Left panel: raw DR25 posteriors chains show a clear (biased)
preference for low values of $b$. Right panel: after reweighting, the
(unbiased) distribution is nearly flat. To minimize spurious peaks and
sampling noise in low probability regions, the lowest density 1% of samples
have been excluded from our reweighting scheme. The slight increase in
probability density near $b\approx 1$ in the reweighted posteriors reflects
the presence of residual importance sampling noise rather than a real feature
of the data. Because there is significant sampling noise (due to the large
implied posterior mass in regions with few samples), our preferred method for
ameliorating the pseudo-density bias is to refit the photometry.

Because $b$ is covariant with $r$ (interacting via non-zero limb darkening),
any bias on $b$ translates to a bias on $r$. For measurements in the final
Kepler data release, DR25, we find this covariance has produced a $1.6\%$
median systematic underestimate of $r$ (Figure 7), extending as high as
$\sim$6% for some targets. This offset is comparable to the fractional
uncertainty on $R_{\star}$ (Gaia Collaboration et al., 2018; Berger et al.,
2018) and so makes up a sizeable portion of the error budget for Kepler
planetary radii. While a few percent difference in planetary radius for a
single planet may be sub-significant, a systematic bias of a few percent on
all planetary radii will significantly impact our interpretation of population
demographics – for example, the precise characteristics of the radius valley
(Fulton et al., 2017) – thereby altering our understanding of the processes by
which planets form and evolve.

Figure 7: Fractional change in median planet-to-star radius ratio for all
planet candidates after correcting posterior chains from DR25 (Thompson et
al., 2018) using the Jacobian reweighting scheme described in §5. $5\sigma$
outliers have been iteratively clipped in order to eliminate spurious values
that are expected to arise due to insufficient sampling of low probability
regions. There is a spike at $\delta_{r}=0$, indicating that some subset of
targets were accurately measured, but the majority of targets are distributed
around $\delta_{r}/r=1.6\%$.

##  6 Summary and conclusions

In this work, we explored the the biases that result from using the popular
stellar pseudo-density, $\tilde{\rho}$, as a parameter in light curve fits.
Adopting a linear prior on this parameter results in a biased distribution on
impact parameter due to the Jacobian that arises from the non-linear
relationship between $\tilde{\rho}$ and transit duration, $T_{14}$. Biased
inferences on $b$ lead to biased inferences on $r$ due to covariances between
the two parameters. We confirmed that the these biases are present in Kepler
modeling that used $\tilde{\rho}$ as a fitting parameter, and we presented a
method for de-biasing the distributions.

Although the $\tilde{\rho}$ bias may be resolved by using $\log\tilde{\rho}$
in place of $\tilde{\rho}$ (or, equivalently, placing log-uniform priors on
$\tilde{\rho}$), we prefer sampling in duration over $\tilde{\rho}$ for
aesthetic and conceptual reasons. To avoid inducing biases, we recommend
sampling directly in duration $T_{14}$ or replacing $T_{14}$ with the true
stellar density and orbital eccentricity vector, i.e.
$\\{\rho_{\star},\sqrt{e}\sin\omega,\sqrt{e}\cos\omega\\}$.

This work focused on the biases induced from using $\tilde{\rho}$ directly as
a fitting parameter; similar biases may arise when using any related
parameterization, for example $a/R_{\star}$, which is a popular choice (e.g.
Crossfield et al., 2015; David et al., 2016; Stassun et al., 2017). As with
$\tilde{\rho}$, adopting a log-uniform prior rather than a linear prior on
$a/R_{\star}$ avoids the unwanted bias. A log-uniform prior is a common
choice, so most analyses which have used $a/R_{\star}$ as a fitting parameter
are probably unaffected by the bias. However, one should always verify what
priors were adopted when interpreting the results of any transit model.

We thank the anonymous referee for reviewing and providing comments which
improved the quality of this manuscript. We are grateful to Jason Eastman, Dan
Fabrycky, Dan Foreman-Mackey, Jason Rowe, Josh Winn, and Jon Zink for helpful
conversations about this work. G.J.G., M.G.M., and E.A.P. acknowledge support
from NASA Astrophysics Data Analysis Program (ADAP) grant (80NSSC20K0457).
E.A.P. acknowledges support from the Alfred P. Sloan Foundation. M.G.M
acknowledges support from the UCLA Cota-Robles Graduate Fellowship. This study
made use of data products from the Kepler mission hosted on the NASA Exoplanet
Archive 2021. Some of the data were obtained from the Mikulski Archive for
Space Telescopes (MAST) at the Space Telescope Science Institute. These data
can be accessed via 10.26133/NEA5 (catalog 10.26133/NEA5). This study also
made use of computational resources provided by the University of California,
Los Angeles and the California Planet Search.

## References

  * Akeson et al. (2013) Akeson, R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125, 989, doi: 10.1086/672273
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Barnes (2007) Barnes, J. W. 2007, PASP, 119, 986, doi: 10.1086/522039
  * Berger et al. (2018) Berger, T. A., Huber, D., Gaidos, E., & van Saders, J. L. 2018, ApJ, 866, 99, doi: 10.3847/1538-4357/aada83
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977, doi: 10.1126/science.1185402
  * Carter et al. (2008) Carter, J. A., Yee, J. C., Eastman, J., Gaudi, B. S., & Winn, J. N. 2008, ApJ, 689, 499, doi: 10.1086/592321
  * Charbonneau et al. (2000) Charbonneau, D., Brown, T. M., Latham, D. W., & Mayor, M. 2000, ApJ, 529, L45, doi: 10.1086/312457
  * Christiansen et al. (2020) Christiansen, J. L., Clarke, B. D., Burke, C. J., et al. 2020, AJ, 160, 159, doi: 10.3847/1538-3881/abab0b
  * Coughlin et al. (2016) Coughlin, J. L., Mullally, F., Thompson, S. E., et al. 2016, ApJS, 224, 12, doi: 10.3847/0067-0049/224/1/12
  * Crossfield et al. (2015) Crossfield, I. J. M., Petigura, E., Schlieder, J. E., et al. 2015, ApJ, 804, 10, doi: 10.1088/0004-637X/804/1/10
  * David et al. (2016) David, T. J., Hillenbrand, L. A., Petigura, E. A., et al. 2016, Nature, 534, 658, doi: 10.1038/nature18293
  * Eastman (2017) Eastman, J. 2017, EXOFASTv2: Generalized publication-quality exoplanet modeling code.  http://ascl.net/1710.003
  * Eastman et al. (2013) Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83, doi: 10.1086/669497
  * Espinoza (2018) Espinoza, N. 2018, Research Notes of the American Astronomical Society, 2, 209, doi: 10.3847/2515-5172/aaef38
  * Foreman-Mackey et al. (2021) Foreman-Mackey, D., Luger, R., Agol, E., et al. 2021, The Journal of Open Source Software, 6, 3285, doi: 10.21105/joss.03285
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109, doi: 10.3847/1538-3881/aa80eb
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gibson et al. (2016) Gibson, S. R., Howard, A. W., Marcy, G. W., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, 990870, doi: 10.1117/12.2233334
  * Gilbert (2022) Gilbert, G. J. 2022, arXiv e-prints, arXiv:2201.08350.  https://arxiv.org/abs/2201.08350
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357–362, doi: 10.1038/s41586-020-2649-2
  * Henry et al. (2000) Henry, G. W., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2000, ApJ, 529, L41, doi: 10.1086/312458
  * Hoffman & Gelman (2011) Hoffman, M. D., & Gelman, A. 2011, arXiv e-prints, arXiv:1111.4246.  https://arxiv.org/abs/1111.4246
  * Kipping (2010a) Kipping, D. M. 2010a, MNRAS, 407, 301, doi: 10.1111/j.1365-2966.2010.16894.x
  * Kipping (2010b) —. 2010b, MNRAS, 408, 1758, doi: 10.1111/j.1365-2966.2010.17242.x
  * Kipping (2014) —. 2014, MNRAS, 440, 2164, doi: 10.1093/mnras/stu318
  * Kipping & Sandford (2016) Kipping, D. M., & Sandford, E. 2016, MNRAS, 463, 1323, doi: 10.1093/mnras/stw1926
  * Kovács et al. (2002) Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369, doi: 10.1051/0004-6361:20020802
  * Luger et al. (2019) Luger, R., Agol, E., Foreman-Mackey, D., et al. 2019, AJ, 157, 64, doi: 10.3847/1538-3881/aae8e5
  * Mandel & Agol (2002) Mandel, K., & Agol, E. 2002, ApJ, 580, L171, doi: 10.1086/345520
  * Mullally et al. (2015) Mullally, F., Coughlin, J. L., Thompson, S. E., et al. 2015, ApJS, 217, 31, doi: 10.1088/0067-0049/217/2/31
  * NASA Exoplanet Archive (2021) NASA Exoplanet Archive. 2021, Kepler Objects of Interest DR25, Version: 2021-08-07 13:20, NExScI-Caltech/IPAC, doi: 10.26133/NEA5
  * Neal (2011) Neal, R. 2011, MCMC Using Hamiltonian Dynamics, 113–162, doi: 10.1201/b10905
  * Rauer et al. (2014) Rauer, H., Catala, C., Aerts, C., et al. 2014, Experimental Astronomy, 38, 249, doi: 10.1007/s10686-014-9383-4
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003, doi: 10.1117/1.JATIS.1.1.014003
  * Rowe et al. (2014) Rowe, J. F., Bryson, S. T., Marcy, G. W., et al. 2014, ApJ, 784, 45, doi: 10.1088/0004-637X/784/1/45
  * Rowe et al. (2015) Rowe, J. F., Coughlin, J. L., Antoci, V., et al. 2015, ApJS, 217, 16, doi: 10.1088/0067-0049/217/1/16
  * Salvatier et al. (2016) Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Comput. Sci., 2, e55.  http://dblp.uni-trier.de/db/journals/peerj-cs/peerj-cs2.html#SalvatierWF16
  * Seager & Mallén-Ornelas (2003) Seager, S., & Mallén-Ornelas, G. 2003, ApJ, 585, 1038, doi: 10.1086/346105
  * Seifahrt et al. (2018) Seifahrt, A., Stürmer, J., Bean, J. L., & Schwab, C. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, ed. C. J. Evans, L. Simard, & H. Takami, 107026D, doi: 10.1117/12.2312936
  * Stassun et al. (2017) Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136, doi: 10.3847/1538-3881/aa5df3
  * Thompson et al. (2018) Thompson, S. E., Coughlin, J. L., Hoffman, K., et al. 2018, ApJS, 235, 38, doi: 10.3847/1538-4365/aab4f9
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Winn (2010) Winn, J. N. 2010, Exoplanet Transits and Occultations, ed. S. Seager, 55–77

##  Appendix A Derivation of Jacobian for $T_{14}\rightarrow\tilde{\rho}$

In this section, we derive the Jacobian of the coordinate transformation
$T_{14}\rightarrow\tilde{\rho}$. The pseudo-density derived by Seager &
Mallén-Ornelas (2003) is

$\tilde{\rho}\equiv\Bigg{(}\frac{4\pi^{2}}{P^{2}G}\Bigg{)}\Bigg{(}\frac{(1+r)^{2}-b^{2}\big{(}1-\sin^{2}[\pi
T/P]\big{)}}{\sin^{2}[\pi T/P]}\Bigg{)}^{3/2}$ (A1)

where all variables are defined as in previous sections. For notational
clarity, we also define $T\equiv T_{14}$ and make the simplifying assumption
$r\approx\sqrt{\Delta F}$, where $\Delta F$ is the fractional change in flux.
Substituting terms

$\displaystyle x=4\pi^{2}/(P^{2}G)$ (A2) $\displaystyle y=(1+r)^{2}$
$\displaystyle z=\sin^{2}[\pi T/P]$

yields

$\tilde{\rho}=x\Bigg{(}\frac{y-b^{2}(1-z)}{z}\Bigg{)}^{3/2}.$ (A3)

By the chain rule,

$\frac{d\tilde{\rho}}{dT}=\frac{d\tilde{\rho}}{dz}\frac{dz}{dT}.$ (A4)

The first term is

$\frac{d\tilde{\rho}}{dz}=-\frac{3x}{2}\Bigg{(}\frac{y-b^{2}}{z^{2}}\Bigg{)}\Bigg{(}\frac{y-b^{2}(1-z)}{z}\Bigg{)}^{1/2}$
(A5)

and the second term is

$\frac{dz}{dT}=\frac{\pi}{P}\sin\Big{[}\frac{2\pi T}{P}\Big{]}$ (A6)

Combining equations A2, A4, A5, and A6 yields the exact Jacobian

$J=\frac{d\tilde{\rho}}{dT}=-\frac{6\pi^{3}}{P^{3}G}\Bigg{(}\frac{(1+r)^{2}-b^{2}}{\sin^{4}[\pi
T/P]}\Bigg{)}\\\ \Bigg{(}\frac{(1+r)^{2}-b^{2}\big{(}1-\sin^{2}[\pi
T/P]\big{)}}{\sin^{2}[\pi T/P]}\Bigg{)}^{1/2}\sin\Big{[}\frac{2\pi
T}{P}\Big{]}$ (A7)

Making the small angle approximation $\sin\phi\approx\phi$ (assuming $\pi T\ll
P$) and collecting terms yields

$J=-\frac{12\pi^{3}}{P^{3}G}\Bigg{(}(1+r)^{2}-b^{2}\Bigg{)}\Bigg{(}(1+r)^{2}-b^{2}\big{(}1-[\pi
T/P]^{2}\big{)}\Bigg{)}^{1/2}\Bigg{(}\frac{\pi T}{P}\Bigg{)}^{-4}.$ (A8)

Once again taking advantage of $\pi T\ll P$ simplifies the expression further
to

$J=-\frac{12\pi^{3}}{P^{3}G}\Bigg{(}(1+r)^{2}-b^{2}\Bigg{)}^{3/2}\Bigg{(}\frac{\pi
T}{P}\Bigg{)}^{-4}.$ (A9)

##  Appendix B Derivation of Jacobian for $T_{14}\rightarrow\ln\tilde{\rho}$

To derive the Jacobian of the transformation $T\rightarrow\ln\tilde{\rho}$, we
note that

$\frac{d\ln\tilde{\rho}}{dT}=\frac{1}{\tilde{\rho}}\frac{d\tilde{\rho}}{dT}.$
(B1)

Adopting our usual approximations $\sin\phi\approx\phi$, $\pi T\ll P$, we may
rewrite Equation A1 in the simplified form

$\tilde{\rho}\equiv\Bigg{(}\frac{4\pi^{2}}{P^{2}G}\Bigg{)}\Bigg{(}(1+r)^{2}-b^{2}\Bigg{)}^{3/2}\Bigg{(}\frac{\pi
T}{P}\Bigg{)}^{-3}.$ (B2)

Combining Equations A9, B1, and B2 and cancelling terms yields

$\frac{d\ln\tilde{\rho}}{dT}=-\frac{3}{T}.$ (B3)

We see that $d\ln\tilde{\rho}/dT$ is independent of $b$.

##  Appendix C Derivation of Jacobian for $T_{14}\rightarrow a/R_{\star}$

To derive the Jacobian of the transformation $T\rightarrow a/R_{\star}$, we
define $\alpha\equiv a/R_{\star}$ and recognize that from Seager & Mallén-
Ornelas (2003) (their Equations 8 & 9),

$\alpha=\Bigg{(}\frac{4\pi^{2}}{P^{2}G}\Bigg{)}^{-1/3}\tilde{\rho}^{1/3}$ (C1)

.

By the chain rule,

$\frac{d\alpha}{dT}=\frac{d\alpha}{d\tilde{\rho}}\frac{d\tilde{\rho}}{dT}.$
(C2)

The first term is

$\frac{d\alpha}{d\tilde{\rho}}=\frac{1}{3}\Bigg{(}\frac{4\pi^{2}}{P^{2}G}\Bigg{)}^{-1/3}\tilde{\rho}^{-2/3}$
(C3)

and the second term we derived previously. Adopting our usual approximations
$\sin\phi\approx\phi$, $\pi T\ll P$ and combining Equations A9, B2, C2, and C3
yields

$\frac{d\alpha}{dT}=-\frac{\pi}{P}\Bigg{(}(1+r)^{2}-b^{2}\Bigg{)}^{1/2}\Bigg{(}\frac{\pi
T}{P}\Bigg{)}^{-2}.$ (C4)

##  Appendix D Derivation of Jacobian for $T_{14}\rightarrow\ln a/R_{\star}$

To derive the Jacobian of the transformation $T\rightarrow\ln a/R_{\star}$, we
note that

$\frac{d\ln\alpha}{dT}=\frac{1}{\alpha}\frac{d\alpha}{dT}$ (D1)

where as before $\alpha\equiv a/R_{\star}$. Following our usual strategies and
combining Equations C1, C4, and D1, we arrive at

$\frac{d\ln\alpha}{dT}=-\frac{1}{T}$ (D2)

which is independent of $b$.