Robust Beam Position Estimation with Photon Counting Detector Arrays in Free-Space Optical Communications

Optical beam center position on an array of detectors is an important parameter that is essential for estimating the angle-of-arrival of the incoming signal beam. In this paper, we have examined the beam position estimation problem for photon-counting detector arrays, and to this end, we have derived and analyzed the Cramér-Rao lower bounds on the mean-square error of the unbiased estimators of the beam position. Furthermore, we have also derived the Cramér-Rao lower bounds of other beam parameters such as peak intensity, and the intensity of background radiation on the array. In this sense, we have considered a robust estimation of the beam position in which none of the parameters are assumed to be known beforehand. Additionally, we have derived the Cramér-Rao lower bounds of beam parameters for observations based on both pilot and data symbols of a pulse position modulation (PPM) scheme. Finally, we have considered a two-step estimation problem in which the peak intensity and background radiation are estimated using a method of moments estimator, and the beam center position is estimated with the help of a maximum likelihood estimator.


I. Introduction
Free-space optical (FSO) communications has typically been employed in deep-space communications due to the low divergence of optical beam that can help transmit data over much longer distances. However, the urgent need to provide connectivity to almost 4 billion people who currently do not have internet access has prompted the integration of FSO and millimeter wave/Terahertz communications in the backhaul of 6G wireless communication systems. Thus, FSO/millimeter/Terahertz wave links will be deployed in an integrated network of satellites, drones, high altitude platform and balloons that will form the backhaul of the next generation of wireless communications. In this regard, FSO is an emerging candidate for future communications due to its ability to support high data rates.
However, the problem of pointing, acquisition and tracking is significant in the FSO (as well as millimeter/Terahertz systems) domain because of the narrow beam widths associated with the optical signal. Acquisition is the process in which the two terminals acquire the initial location of each other before the actual data communication begins. However, after the acquisition is achieved, the system still needs to maintain the alignment between the transmitter and receiver assemblies due to physical factors such as random effects associated with atmospheric turbulence, the mechanical jitter introduced in the transmitter/receiver assemblies due to outside disturbances, or building sways receiver. Finally, the authors in [11] consider time synchronization schemes based on an array of detectors.
Furthermore, we also briefly discuss the literature on pointing and tracking in FSO systems that examine the tracking problem from the perspective of a single detector. In this regard, [12] develops the pointing error statistics for a circularly shaped detector and a Gaussian beam, and the outage capacity is optimized as a function of beam radius. The authors in [13] investigate a slightly different optimization problem concerning pointing: They have considered the maximization of link availability as a function of beam radius (for fixed signal power). Additionally, they also explore the minimization of transmitted power by tuning to the optimal beam radius under the constraint of a fixed link availability. In addition to these papers, the interested reader may be directed to [14]- [17] for a detailed study on the performance of FSO systems when the optical channel suffers degradation due to pointing errors for a single detector receiver.
Readers who might be interested in deep space optical communications with photon-counting detector arrays are referred to [18]- [22].

B. Contributions of This Paper
Even though the authors in [6] have proposed beam position estimation algorithms with an array of detectors, they have not considered the derivation of the Cramér-Rao lower bounds on the variance of unbiased estimators. We believe that an understanding of the Cramér-Rao bounds of beam position and other parameters is an important problem as these bounds will give us important insights into the behavior of estimators under different channel conditions. Additionally, beam parameters other than the beam position-such as the peak intensity, the beam radius and the background radiation level 3are assumed by these authors as known quantities in their derivation of estimators. However, even though the beam radius on the focal plane array may be considered constant, the peak intensity and the background radiation intensity may change significantly with time, and they have to be estimated in real-time in order to improve the performance of beam position estimators. Additionally, it goes without saying that a knowledge of peak intensity and background noise is important for allocating power in different channels of a multiple-input-multiple-output (MIMO), or a multiple-input-singleoutput (MISO), FSO communication system as these two quantities specify the signal-to-noise ratio of the channel [23].
In this study, we have derived and analyzed the Cramér-Rao lower bounds of the beam position estimators for an array of detectors. The estimation problem discussed in this paper is robust since we also estimate the signal intensity as well as background radiation/noise power levels. In this regard, Cramér-Rao lower bounds are derived for the beam position, peak intensity and noise intensity for various scenarios. Moreover, two types of observations are considered in this estimation problem: i) observations based on pilot symbols and ii) observations based on data symbols. Using energy of the signal based on data symbols for our estimation problem leads to a more bandwidth/energy efficient scheme. However, as we will see later in this study, the estimation performance (in terms of meansquare error) corresponding to data symbols suffers more at poorer signal-to-noise ratio as compared to the pilot symbols.
For observations based on data symbols, we additionally consider two different observation intervals for estimation: i) a pulse position modulation symbol period, or, ii) the PPM slot period containing the signal. The estimation based on the slot period assumes that the signal pulse is present in the slot. If this assumption is true, the signal-to-noise ratio in the slot is much higher (as much as K times) than the signal-to-noise ratio in one K-PPM symbol period. This high signal-to-noise ratio leads to a better estimation performance. However, if the receiver makes an error, we use the "noise only" slot as our observation which results in a severe degradation in estimation performance. We discuss these ideas in more detail in Section V.

C. Model Assumptions
The major assumption regarding beam position estimation with detector arrays is that the array area is large enough so that the beam footprint is smaller than the footprint of the array. This assumption is required in order to avoid any outage in the received signal in case the beam wanders within the bounds of the array. Furthermore, a large array area is also required in order to track the movement of the beam and align it-perhaps through a fast steering mirror (FSM) assembly-to the center of the array.
Secondly, all arrays are assumed to be of a square shape and each detector in the array is also assumed to be of a square shape as well.
Finally, the focus of this study is on non-Bayesian estimation techniques for beam position estimation. This is due to the fact that unless we are certain about the parameters of the prior random motion model of the beam on the array, we are likely going to incur a significant loss in performance if there is mismatch in our assumptions and the real world parameters 4 [24].

D. Organization of This Paper
This paper is organized as follows. Section III defines the beam profile and the Poisson model that governs the occurrence of photodetections 5 in the array of detectors. Section IV discusses the derivation of the Cramér-Rao lower bound of the beam parameter estimation problem when pilot symbols are used as an observation. In Section V, we derive the Cramér-Rao bounds for observations based on data symbols. Section VI considers the two-step estimation (method of moments and maximum likelihood estimators) algorithm to estimate beam parameters. The simulation results are explained in Section VII and Section VIII briefly discusses the complexity of the two estimators. The conclusions of this study are summarized in Section IX.

III. System Model
The received optical signal on the receiver aperture gives rise to photoelectrons or photodetections in each detector of the array due to photoelectric effect. The emission of these photoelectrons during the signal pulse interval helps us detect transmitted symbols. The photon count Z m in the mth detector or cell of the array-during some specified observation interval-is modeled as a (Poisson) discrete random variable. Its probability mass function is characterized by the following expression: where λ s (x, y) is the scaled beam intensity 6 profile on the detector array, λ n is the scaled noise intensity profile, A m is the region of the mth detector on the detector array, Z 1 , Z 2 , . . . , Z M are independent Poisson random variables and M is the total number of detectors in the array. As may have been discerned by the reader, the coordinate (x, y) stands for any point inside the region of the detector array. Moreover, λ n is a constant factor that accounts for the background radiation and the thermal effects of the detector array [25]. We assume that the airy pattern of the beam on the focal plane array is well-approximated by a Gaussian function. The received (scaled) signal and noise intensity at the detector array is given by the expression where I 0 /ρ 2 is the peak intensity measured in terms of number of signal photons measured during an observation interval. Furthermore, λ n is also measured in terms of number of noise photons generated during the same observation interval. The quantity ρ is known as the beam radius measured in mm, and (x 0 , y 0 ) is the center of the Gaussian beam on the detector array. The function 1 A (·) represents the indicator function over some (measurable) set A, and A is the region of the detector array. Furthermore, it is a general assumption in the following sections that the center of the array has the coordinates (0, 0). Additionally, the area of A m is denoted by A since all detectors are assumed to have an equal area. The area of the detector array is denoted by |A|. The length of one side of the array is denoted by (A).

IV. Cramèr-Rao Lower Bounds for Beam Parameter Estimation Based on Pilot
Symbols In this section, we derive and analyze the Cramér-Rao lower bounds for the beam parameter estimation problem based on pilot symbols. The pilot symbol is transmitted as a known pulse position modulation symbol. For instance, we may transmit only the '0' symbol (signal pulse only in the first slot) of a K-PPM scheme. The observation interval in this case is the first slot of every pilot symbol.
Let θ x 0 y 0 I 0 λ n 7 . The likelihood function is given by where Λ m and the random vector Z Z 1 Z 2 · · · Z M T . Let us define the total incident power on the array Λ s M m=1 Λ m . Then,

A. Estimation Based on Pilot Symbols
In this section, we derive Cramér-Rao lower bounds based on the observations corresponding to pilot symbols. 1) Cramèr-Rao Lower Bound of I 0 As a first step in computing the Cramér-Rao lower bound for any unbiased estimatorÎ 0 , we compute the first partial derivative of (5): and then, Now, taking the expectation with respect to Z m and taking the negative of the resulting quantity, we have that 2) Cramèr-Rao Lower Bound of λ n In this case, we assume that the background radiation is estimated at the receiver while the transmitter is turned off (no signal is present at the receiver). Therefore, in this case, Λ m = λ n A. Using the same line of argument as used in the derivation of (8), it can be easily shown that the Cramèr-Rao lower bound on the variance of any unbiased estimatorλ n for one noise-only slot is given by 3) Cramèr-Rao Lower Bounds of (x 0 , y 0 ) The Cramèr-Rao lower bound forx 0 andŷ 0 is derived in Section A of the appendix. The final expressions are produced here as follows: and where Ψ(x 0 , y 0 , I 0 , ρ) is defined in (65).

4) Cramèr-Rao Lower Bounds for Joint Estimation of
In this section, we state the Cramér-Rao lower bounds for the three parameter estimation problem in which the three beam parameters are x 0 , y 0 and I 0 8 . We denote the 3×3 Fisher Information Matrix by I(x 0 , y 0 , I 0 ). The Cramér-Rao lower bounds are given by In order to lower the complexity of the estimation problem, we can estimate λn independently of x0, y0 and I0. In this case, all we need to do is to estimate the average number of noise photons by occasionally turning the transmitter off.
The determinant of the Fisher information matrix is given by where We know that each detector in the array counts or reports the photodetections that occur inside its region in a given observation interval for the purpose of beam position estimation. However, the detector does not specify the exact location of the photodetection inside its region. In the ideal case when M → ∞ for fixed array area, the true location of each photodetection can be reported by the infinitesimally small detector. When M → ∞, we call this limiting array a "continuous" array. This asymptotic case is of interest since the probability of error/tracking performance of a practical array can be reasonably approximated with the continuous array when the number of detectors is large enough (M ≥ 64) [6], [9]. Therefore, in this section, we look at the Cramér-Rao lower bound of (x 0 , y 0 ) for the M → ∞ case for the low and high signal-to-noise-ratio regimes, and the convergence rates of the Cramér-Rao lower bounds are derived in terms of beam radius ρ.
In the following analysis, let us analyze the Cramér-Rao lower bound ofx 0 only since the same analysis will holdŷ 0 due to the symmetric nature of the Gaussian beam.

1) Estimation of x 0 : High Signal-to-Noise Ratio
For high signal-to-noise ratio, λ n A << Am is the center of the mth small detector, and ∆ M is its infinitesimal area. Then, the numerator of (10) simplifies as where, in the last approximation of (19) we have used the fact that . The positive term in the denominator (see (65)) can be simplified in a similar fashion. The square root of the term with minus sign can be simplified as where in the last approximation of (20), we have used the fact that we have used the fact that where X and Y are independent Gaussian random variables with the same variance ρ 2 , but with different means: We note that the Cramér-Rao lower bound is minimized by minimizing ρ (a more focused beam) for fixed signal power. The Cramér-Rao lower bounds goes to zero at the rate O(ρ 2 ) as ρ → 0, where O represent the "big O" notation. Moreover, the Cramér-Rao lower bounds goes to zero in terms of I 0 at the rate O(I −1 0 ).

2) Estimation of x 0 : Low Signal-to-Noise Ratio
In this case, let us assume that λ n A >> Am In this case, the square root of the term with the minus sign in the denominator is which is zero due to the symmetric nature of the Gaussian beam. Therefore, by further simplification, which goes to as M → ∞. In this case, the Cramér-Rao lower bound goes to zero at a rate O(ρ 4 ) as ρ → 0. This is a faster rate of convergence than O(ρ 2 ) for the high signal-to-noise ratio case. Additionally, Cramér-Rao lower bound converges to zero at the rate O(I −2 0 ) in the low signal-to-noise ratio regime.

V. Cramèr-Rao Lower Bounds for Beam Parameter Estimation Based on Data Symbols
Since pilot symbols incur a loss in energy and bandwidth, there is a motivation to use data symbols for the estimation of beam parameters even though that may result in some loss in estimation performance. In this section, we derive the Cramér-Rao lower bounds of beam parameters based on PPM data symbols. As discussed in Section II-B, we can either use either one PPM symbol period or a signal slot period as our observation interval. We first look at the Cramér-Rao lower bounds related to the symbol period based observation in the next section.

A. Observations Based on Symbol Period
In this case, the noise power goes up K times where K is the number of slots in PPM. Thus, the new λ n Kλ n , and Thus, in this case,

B. Observations Based on Slot Period: A Decision-Aided Scheme
The motivation behind choosing the slot period is to maximize the signal-to-noise ratio in the sufficient statistic. If the slot containing the signal is chosen, the resulting signal-to-noise ratio is K times bigger than the signal-to-noise ratio available in a symbol period. However, for the slot period case, we depend on the correct decision of the receiver to choose the "right" slot that contains the signal. If the receiver makes a mistake, we end up choosing a "noise-only" slot and the resulting noise photons do not give us any information about the beam parameters. Thus, if the receiver starts making too many errors, the estimation performance will take a significant hit. Thus, in the slot period case, the correct symbol decision is the key to a good estimation performance, and we term the estimation based on slot period alternatively as decision-aided estimation of beam parameters. Fig. 5 shows the block diagram of the decision-aided beam position estimation scheme in which the output of the equal gain combiner is fed into the beam position estimation block. If the equal gain combiner declares some symbol j as the transmitted (K-PPM) symbol for 0 ≤ j < K, the beam position estimation block chooses the jth slot as its observation interval.
For observation based on one slot, Z m ∼ P oisson (Λ m ) , with probability P c and Z m ∼ P oisson(λ n A) with probability (1 − P c ). Thus, Z m is a doubly stochastic Poisson process or a Cox process whose intensity varies randomly according to the Bernoulli distribution as follows: where δ(·) is the Dirac delta function. Therefore, the likelihood function becomes The quantity P c is the probability of a correct decision of the equal gain combiner. It can be shown that for a maximum a posteriori probability detector that operates on a K-PPM symbol, we have that In (29), Z s ∼ P oisson(Λ s ) and Z n ∼ P oisson(λ n |A|). The random variable Z Z s −Z n is a (discrete) Skellam random variable whose distribution is where I z (·) is the modified Bessel function of the first kind (not to be confused with peak intensity I 0 ). Thus, Fig. 4 shows the probability of correct decision P c for different values of beam radius ρ. A large beam radius results in some loss of energy since some of the beam energy falls off the edge of the array. This leads to a lower probability of correct decision for larger beam radii. Finally, since P e = 1 − P c , we have that (31)

1) Monte Carlo Expectation
It is not straightforward to compute the probabilistic expectations E ∂ 2 ln p(Z|θ) and E ∂ 2 ln p(Z|θ) ∂x 0 ∂y 0 for the likelihood function in (28). Thus, we resort to the Monte Carlo simulations to compute these expectation. The simulations are carried out as follows: 1) Sample 1 with probability P c and 0 with probability 1 − P c .

3) Substitute the Z m 's obtained from
Step 2 into each of the second order partial derivatives:

VI. Two-Step Estimation of Beam Parameters
In this section, we will look at a two-step estimation algorithm that is used for estimating the beam parameters. The two-step estimation algorithm is as follows: 1) In the first step, the peak intensity I 0 and background radiation λ n are estimated using a method of moments estimator. 2) The estimatesÎ 0 andλ n obtained from Step 1 are substituted into the loglikelihood function ln p(Z|θ) and the estimate of (x 0 , y 0 ) is obtained by maximizing the loglikelihood function (maximum likelihood estimation). Alternatively, all the four parameters (x 0 , y 0 , I 0 , λ n ) can be estimated via the maximum likelihood estimator. However, since no closed-form expressions for the maximum likelihood estimator are available, we have to resort to numerical optimization techniques (such as a genetic algorithm) in order to find the peak of the loglikelihood function. This incurs a much higher computational complexity if all the four parameters are estimated with the maximum likelihood estimator. The two-step estimation algorithm reduces the complexity since two of the four parameters (I 0 and λ n ) can be estimated via the computationally efficient method of moments estimator without any knowledge of (x 0 , y 0 ), and the numerical search for the maximum of loglikelihood function is limited to just two dimensions in order to find (x 0 ,ŷ 0 ).

A. Method of Moments Estimator of I 0 and λ n 1) Pilot Symbol Case
The method of moments estimator of I 0 for the pilot symbol case iŝ where Z

(s)
i,m is a Poisson random variable with mean Λ m . The method of moments estimator of λ n iŝ where Z (n) i,m is a Poisson random variable with mean λ n A. It can be easily shown that E[λ n ] = λ n , and Thus, bothÎ 0 andλ n are unbiased estimators of I 0 and λ n , respectively. It is straightforward to verify that The mean-square error betweenÎ 0 and I 0 is

2) Decision-Aided Estimation (Observations Based on Slot Period)
In this case, the generation of photon counts are governed by a doubly stochastic Poisson process. Thus, where E[Z m |c] = Λ m = Λ (s) m + λ n A and E[Z m |e] = λ n A. Therefore, Therefore, Moreover, Finally, since Var Î 0 − I 0 = Var Î 0 , we have that B. Maximum Likelihood Estimation of (x 0 , y 0 ) For the pilot symbol scheme, the maximum likelihood estimator of beam position (x 0 , y 0 ) on the array is given by [6]: where Φ(·) is the cumulative distribution function of a standard normal random variable, andÎ 0 and λ n are the method of moments estimates of I 0 and λ n , respectively. The quantity (x m 2 , y m 2 ) is the location of the north east corner of the mth detector, and (x m 1 , y m 1 ) is the position of south west corner.
For estimation based on data symbols, the maximum likelihood estimate is obtained by maximizing (26) (symbol period) or by maximizing (28) (slot period).
Regarding the maximization of the loglikelihood functions (3), (26) and (28), we have to resort to a genetic algorithm to look for the global maximum.
The mean-square error of the maximum likelihood estimator is computed via Monte Carlo simulations. The average of the squared errors is computed by repeating the experiment many times and then computing the sample average of the squared errors.

VII. Simulation Results and Discussion
In this section, we interpret the simulations results that we have obtained in this study. In these simulations, we have considered the low photon rate regime. In this regard, we have considered-on average-10 signal photons for the entire array during the observation interval. These low photon rate channels are of interest in deep space communications where the received signal energy is so low that we are only able to detect a few signal photons per slot of a PPM symbol [1], [6]. Additionally, the low rate of photons has also to do with the "blocking" phenomenon of avalanche photodetectors that are operated in Geiger mode as photon counters. The blocking occurs because the detection of the first signal photon causes an avalanche of electrons, and this avalanche has to be quenched by an avalanche recovery circuit and the bias has to be restored before the next photon can be detected. Thus, the detector "sleeps" or gets "blocked" for a few microseconds before it is ready to detect the next incoming photon.
Moreover, in the ensuing discussion, λ n is measured in terms of average number of noise photons that occur during an observation interval.
For all the experiments, the area of the detector array |A| = 4 mm 2 . This detector array extends from -1 mm to 1 mm in each of the two dimensions, and the center of the array coincides with the origin. Additionally, we want to emphasize that the area of the array |A| is fixed irrespective of the number of detector M in the array. Thus, a larger M implies a smaller area per detector. In terms of notation, we want to point out that the expression CRLB(x 0 ,ŷ 0 ) denotes the sum of individual Cramér-Rao lower bounds-CRLB(x 0 ) and CRLB(ŷ 0 )-wherex 0 andŷ 0 are any unbiased estimators of x 0 and y 0 , respectively. This is true since these two parameters can be treated independently of each other due to circularly symmetric nature of Gaussian beam. Fig. 6 indicates the Cramér-Rao lower bound plots for (x 0 ,ŷ 0 ,Î 0 ) as a three parameter estimation problem as defined in (12), (13) and (14). These curves are plotted as a function of noise parameter λ n . Fig. 7 depicts the Cramér-Rao lower bound curves as a function of number of pilot symbols used as observations in the estimation of parameters. Fig. 8 shows the Cramér-Rao lower bound plots as a function of beam radius ρ at the point (x 0 , y 0 ) = (0.4, 0.4). We note that the Cramér-Rao lower bound ofÎ 0 increases monotonically with ρ. However, for the CRLB of (x 0 ,ŷ 0 ), we see that there is an optimum value of ρ (lets call it ρ * M for the M -cell array) at which the Cramér-Rao lower bound is minimized. Additionally, ρ * N < ρ * for N > M . Intuitively, these observations are straightforward to explain. For a fixed signal-to-noise ratio, if the beam footprint is small, but at least covers one detector completely, then such a small beam footprint will minimize the mean-square error. This is true since all the power is focused into a small region on the array where the number of noise photons (on average) is relatively small, and this fact will help the estimator in order to estimate the beam position more accurately as opposed to a more "spread out" beam. However, if the beam radius is much smaller than the dimensions of a single detector, then the beam will only give rise to photons in the detector in which it is located, and the neighboring detectors will not register any signal photons. Since we round off the locations of the photons-that occur inside a given detector-to the center of that detector, any movement of the "super thin" beam inside the given detector cannot be tracked. Therefore, the Cramér-Rao lower bound rises if ρ diminishes beyond a certain (optimum) value. Fig. 9 and Fig. 10 show the effect of beam radius on the Cramér-Rao lower bounds of (x 0 ,ŷ 0 ) andÎ 0 as a function of (x 0 , y 0 ) for a 4 × 4 detector array when pilot symbols are used for parameter For (x 0 ,ŷ 0 ), we note that for small ρ, the Cramér-Rao lower bound is highly sensitive to the location of the beam center (x 0 , y 0 ). For example, if we consider the case for ρ = 0.12 mm, we note that the diameter of the beam 2ρ is much smaller than the breadth of a single detector in this case (the breadth of a single detector for a 4×4 array is (A m ) = 0.5 mm). Thus. corresponds to the edges of detectors. We additionally note that the Cramér-Rao lower bound attains its peak value at the centers of detectors and its minimum values on the edges. This pattern is explained by our earlier understanding (as argued during elaboration of Fig. 8) that when the beam is very thin and its center resides on the center of a particular detector, then all the energy of the beam resides in that particular detector, and a small movement of the beam cannot be detected by the array. However, if the beam center of such a thin beam lies on the edge of a detector, its slightest movement can be tracked by detecting the change in the energy difference of the two detectors that share the edge. For the case of M = 64, (A m ) = 0.25 mm, and since 2ρ ≈ (A m ), we observe that the fluctuation of the Cramér-Rao lower bound as a function of (x 0 , y 0 ) is almost negligible. When the beam radius is larger such that 2ρ > (A m ), then the energy of the beam is not confined to a single detector regardless of where the beam center resides on the array. In this case, any movement of the beam will be registered because of change in the detected energy reported by the detectors. However, whether the beam diameter is large or small, once the beam center gets too close to the edge of the array, part of the beam energy will fall off the edge of the array and the detector array will experience a net loss in received signal energy. This leads to a higher Cramér-Rao lower bound at the edges of the detector array. Additionally, even though the fluctuation in Cramér-Rao lower bound is minimized for a higher beam radius, the minimum value of the Cramér-Rao lower bound suffers significantly (the minimum value of Cramér-Rao lower bound goes up) as compared to the smaller beam radius scenario. This is because with a larger beam radius, the signal energy is spread out to a larger number of detectors, which results in a smaller signal-to-noise ratio per detector. With a smaller beam radius, the beam energy is confined to a smaller number of detectors, and the average signal-to-noise ratio in these detectors remains significantly higher compared to the large beam radius case. This observation corroborates our previous assertion that the Cramér-Rao lower bound decays as O(ρ 2 ) and O(ρ 4 ) for the high and low signal-to-noise ratio, respectively, as ρ → 0 (see (21), and (24)).
Regarding the Cramér-Rao lower bound ofÎ 0 , we observe that the Cramér-Rao lower bound is not very sensitive to the beam radius and the performance does not change significantly over the chosen range of beam radii. Additionally, the Cramér-Rao lower bound performance between the 4×4 detector and the 8 × 8 case is not significant. Fig. 11 illustrates the Cramér-Rao lower bound comparisons for the observations based on pilot symbols versus data symbols. We note that, as expected, the Cramér-Rao lower bound computed based on pilot symbols outperforms the Cramér-Rao lower bound based on data symbols (both for slot period and symbol period). We also note that the Cramér-Rao lower bound based on slot period outperforms the Cramér-Rao lower bound related to symbol period for high signal-to-noise ratio. However, the situation reverses for low signal-to-noise ratio in which case the symbol period observation does a better job than the slot period in terms of minimizing the Cramér-Rao lower bound. Fig. 12 show the performance comparison of the maximum likelihood estimator for the pilot and data symbol based observations, and we see a similar trend as the Cramér-Rao lower bound curves in Fig. 11. Fig. 13 and Fig. 14 depict the performance of the method of moment estimator of λ n and I 0 , respectively. The performance of the method of moments estimator is compared with the Cramér-Rao lower bound as well. We note in Fig. 14 (as was the case with the Cramér-Rao lower bounds of (x 0 ,ŷ 0 )) that the Cramér-Rao lower bounds based on slot period is smaller than the Cramér-Rao lower bounds based on symbol period when λ n is less than a particular value, and that the observations based on symbol period outperform the observations related to slot period when λ n exceed that particular (critical) point. Thus, for a system that employs data symbols for estimation of beam parameters, we propose a switching scheme in which the estimation algorithm chooses either the symbol period or the slot period based on whether the value of λ n is higher or lower than the critical value.

VIII. A Brief Complexity Analysis of Estimators
It is easy to see from (32)  number genetic algorithm is discussed in detail [6]. We note from [6] that the complexity of the genetic algorithm is a function of number of chromosomes, N c , and the number of generations 9 , N g . The values of N c and N g have to be chosen according to the nature of the objective function-a "spikier" function requires relatively large N g and N c for convergence to the true maximum. Additionally, the larger the number of dimensions of the objective function, the larger the values of N c and N g become in order to speed up convergence. In our simulations, we set N c = 50, and N g = 400. Since the value of the loglikelihood function is computed for each chromosome, the total complexity of the maximum likelihood estimator is approximately N c × N g × O(M ) real multiplications and real additions 10 . The interested reader is referred to [26] for more details on genetic algorithms. . The average signal photon rate is 10 photons/slot period, and the beam radius is 0.2 mm. The performance of the maximum likelihood estimator is compared to the Cramér-Rao lower bound for the pilot symbol observations. IX. Conclusion In this paper, we have analyzed the Cramér-Rao lower bounds for the robust beam position estimation problem for a deep space optical communication system that uses an array of photon counting detectors at the receiver. In this regard, we have derived the Cramér-Rao lower bounds for observations based on pilot symbols as well as data symbols (symbol period and slot period) of the pulse position modulation scheme. Even though the pilot symbols provide a superior performance in terms of estimation of beam parameters, we pay an extra price in terms of a larger bandwidth and a higher energy expenditure for transmission of dedicated pilot symbols. With data symbols, no extra bandwidth and energy expenditure is required. Additionally, for estimation of beam parameters with data symbols, we propose a switching scheme that switches between the symbol period and slot period based observations depending on the noise level at the receiver.