Optical ISAC with Single-Photon Avalanche Diode (SPAD) Arrays

Introduction

Optical integrated sensing and communication (O-ISAC) recently emerged as the optical wireless counterpart of integrated sensing and communication (ISAC) in radio-frequency (RF) systems, reusing one intensity-modulated waveform for data transmission and environment sensing [1, 2]. Direct-current-biased optical orthogonal frequency-division multiplexing (DCO-OFDM) supports both high-rate communication and delay-based sensing [2]. As in RF ISAC, the sensing performance of such a system is characterized by the Cramér-Rao bound (CRB) and the CRB-rate region [3, 4, 5].

Corner-cube reflectors (CCRs) strengthen the echo of cooperative targets [6, 7]. Diffuse targets, in contrast, re-radiate the intercepted light as a Lambertian source, so that two spreading losses compound into a $1/d^{4}$ round-trip loss at distance $d$ [7, 8]. The low echo power motivates single-photon avalanche diode (SPAD) arrays as the sensing receiver [9]. Prior O-ISAC work, however, employs linear receivers [1, 2, 6, 7, 10], and how the nonlinearity of a photon-counting front end affects the sensing-communication trade-off has not been analyzed.

In SPAD-based optical wireless communication (OWC), the dead time compresses and eventually paralyzes the response and makes the counts sub-Poisson [9, 11]. In single-photon light detection and ranging (LiDAR), pile-up [12, 13] and the resulting range bias, similar to the range-walk error of pulsed LiDAR, are compensated by dead-time-aware estimators [14, 15]. For pulsed SPAD LiDAR, an optimal incident photon rate is also known, both from an attenuation rule [16] and from a CRB under dead time [15]. None of this transfers directly to ISAC, where the photon rate follows a data-bearing DCO-OFDM intensity rather than a pulse train. As a result, the sensing optimum is coupled with the rate through the shared transmit power. Moreover, the SPAD count models of OWC [9] assume a photon rate that is quasi-static over a dead time. Ranging, however, demands a sampling interval as short as the dead time, where this assumption fails and the pile-up correction must be used instead [12, 13].

Motivated by the above, we analyze ranging and the rate-constrained transmit power of a monostatic DCO-OFDM O-ISAC system that combines the sensing channels of CCRs and diffuse targets [6, 7] with the passive-quenching SPAD-array model of [9]. We formulate the pile-up-corrected ranging observation and derive the joint delay-amplitude Fisher information and the range CRB. Through simulations, we show that the maximum-likelihood (ML) estimator based on this model is unbiased and, up to moderate saturation, close to the bound, whereas an estimator based on the quasi-static model exhibits a power-dependent range bias. Finally, we show that, once the echo saturates, the ranging accuracy of the SPAD array degrades beyond a sensing-optimal transmit power, unlike for a linear photodiode (PD) receiver, which benefits from more power in both functions.

System Model

We consider the system in Figure 1, with a dual-function base station (BS) mounted on the ceiling and facing downwards, which carries a light-emitting diode (LED) transmitter and a co-located SPAD array with $N_{\mathrm{a}}$ pixels. A communication user at line-of-sight range $d_{\mathrm{c}}$ is equipped with a conventional PD facing upwards, and a passive target at range $d_{\mathrm{s}}$ on the axis of the BS reflects the transmitted waveform back to it [2]. The direct LED-to-array leakage is assumed to be optically suppressed, with its residual absorbed in the background rate, as in the monostatic O-ISAC literature [2], [6].

DCO-OFDM Waveform and Communication Link

A DCO-OFDM block consists of the $N$ intensity samples [9]

for $n = 0,\ldots,N-1$, where $s[n]$ is the real-valued multicarrier signal obtained from $N$ Hermitian-symmetric subcarriers with $M$-ary quadrature amplitude modulation ($M$-QAM). The bias is $B = \beta\sigma_s$, with $\sigma_s$ being the standard deviation of $s[n]$ and $\beta$ the normalized bias¹. The operator $(\cdot)^{+} = \max\{\cdot,0\}$ enforces the non-negativity of the intensity. The block is preceded by a cyclic prefix (CP) of $L_{\mathrm{cp}}$ samples, and consecutive blocks are transmitted back to back. The samples are spaced by $T_s$, and the LED emits the intensity $x(t) = \left(s(t) + B\right)^{+}$, with $s(t)$ being the band-limited interpolation of $\{s[n]\}$, so that $x(nT_s) = x[n]$. The average emitted optical power is $P_t = \mathbb{E}[x(t)] \approx B$, where $\mathbb{E}[\cdot]$ denotes the expectation over the data symbols. The PD user observes $y_c[n] = h_cx[n] + w_c[n]$. Here, $h_c = (m+1)A_{\mathrm{pd}}\rho_{\mathrm{pd}}\kappa_{\mathrm{opt}}\cos^m\phi_c\cos\psi_c/(2\pi d_c^2)$ is the Lambertian line-of-sight gain [8], where $m=-\ln2/\ln\cos\Psi_{1/2}$ is the Lambertian order for half-power semi-angle $\Psi_{1/2}$, $\phi_c$ and $\psi_c$ are the irradiance and incidence angles, and $A_{\mathrm{pd}}$, $\rho_{\mathrm{pd}}$, and $\kappa_{\mathrm{opt}}$ denote the PD area, its responsivity, and the optical filter-concentrator gain, respectively. Moreover, $w_c[n]$ is the zero-mean additive white Gaussian noise of the thermal-noise-limited PD front end, with power spectral density $N_0$ [6] and variance $\sigma_c^2 = N_0B_{\mathrm{occ}}$, where $B_{\mathrm{occ}} = (N/2-1)/(NT_s)$ is the occupied bandwidth. The electrical signal-to-noise ratio (SNR) per data subcarrier is thus $\mathrm{SNR}_c = h_c^2\sigma_s^2/\sigma_c^2$.

Sensing Channel and Photon Rate

The CCR and diffuse channel gains at range $d_s$ are given by [6], [7]

respectively, where $A_{\mathrm{spad}}$ is the array aperture, $k_{\mathrm{eff}}$ denotes the effective CCR reflectance, and $A_{\mathrm{tgt}}$ and $\rho_{\mathrm{diff}}$ are the area and the reflectivity of the diffuse target, respectively. In addition, $\phi_s$ and $\psi_s$ are the irradiance and incidence angles of the sensing link, and $\psi_t$ and $\chi$ denote the incidence and re-radiation angles at the diffuse target. Let $\tau=2d_s/c$ denote the round-trip delay, with speed of light $c$, and let $D=c\tau/2$ be the range to be estimated. The echo reaches the array with optical power $\zeta x(t-\tau)$, $\zeta\in\{\zeta_{\mathrm{CCR}},\zeta_{\mathrm{diff}}\}$, so that the incident photon rate is [9]

where $\alpha_s=\zeta C_s$. Here, $C_s=\Upsilon_{\mathrm{PDE}}(1+\phi_{\mathrm{AP}}+\phi_{\mathrm{CT}})/E_{\mathrm{ph}}$ converts optical power into the photon rate, with $\Upsilon_{\mathrm{PDE}}$ being the photon detection efficiency, $\phi_{\mathrm{AP}}$ and $\phi_{\mathrm{CT}}$ the afterpulsing and crosstalk probabilities, and $E_{\mathrm{ph}}$ the photon energy at the operating wavelength. The rate $\lambda_b=(\vartheta_{\mathrm{DCR}}+\Upsilon_{\mathrm{PDE}}P_B/E_{\mathrm{ph}})(1+\phi_{\mathrm{AP}}+\phi_{\mathrm{CT}})$ collects the dark counts of the array, for total rate $\vartheta_{\mathrm{DCR}}$ and background light power $P_B$. As in [9], afterpulsing and crosstalk are counted as additional independent arrivals, and the array is uniformly illuminated, so that the photon arrivals at each of its $N_a$ pixels form an independent Poisson process of rate $\lambda(t)/N_a$.

SPAD Array and Pile-Up-Corrected Binomial Model

After each detected photon, a SPAD pixel is quenched and remains blind to incident photons for a fixed interval, the dead time $\tau_d$ [9]. With passive quenching, the standard low-complexity array implementation, every photon arriving during the dead time restarts it, i.e., the dead time is paralyzable, the more severe form of saturation [9], [11]. For a constant photon rate $\lambda$, the mean detected photon count of the array per sampling interval $T_s\geq\tau_d$ and its variance are [9]

respectively, which we refer to as the quasi-static model. The mean peaks at $\lambda^\star=N_a/\tau_d$ and decreases beyond it, i.e., the array paralyses, and the variance is sub-Poisson [11]. The normalized photon rate $\epsilon=\lambda\tau_d/N_a=\lambda/\lambda^\star$ measures the saturation, with $\epsilon=1$ at the peak.

The array counts continuously, and the count of a pixel in a sampling interval depends on the photons of the preceding dead time. At the start of a block, these belong to the echo that precedes it. We therefore choose the CP longer than the round-trip delay plus one dead time, $\tau+\tau_d\leq L_{\mathrm{cp}}T_s$, and use the $N$ sampling intervals that follow it, $t\in[0,NT_s)$, as the ranging observation, so that their counts depend only on the data of the current block.

Ranging calls for a short sampling interval, as the range bin is $cT_s/2$, and we set it equal to the dead time, $T_s=\tau_d$. This is the shortest interval for which a dead time reaches back only into the preceding interval, and the longest for which a pixel detects at most one photon per interval. Let $y[n]$ denote the number of pixels that detect a photon in the $n$-th interval $[nT_s,(n+1)T_s)$, $n=0,\ldots,N-1$, and $\tilde q[n]$ the probability that a given pixel does so, i.e., the firing probability. Since the pixels have independent and identical arrival processes, the array count follows the binomial distribution with $N_a$ trials and success probability $\tilde q[n]$, which we write as

A paralyzable pixel detects a photon arriving at $t$ only if no photon arrived in $(t-\tau_d,t)$ [11]. For Poisson arrivals of rate $\lambda(t)/N_a$, the probability of no arrival in $(t-\tau_d,t)$ is $\exp\left(-\frac{1}{N_a}\int_{t-\tau_d}^{t}\lambda(u)\,\mathrm{d}u\right)$. Since a pixel detects at most one photon per interval, its firing probability equals its mean number of detections in the interval, i.e., the pile-up integral [12], [13]

with $\lambda(t)$ given in (3). For a photon rate that is constant over the interval and the dead time before it, (7) reduces to the quasi-static value $\mu_a(\lambda[n])/N_a$ of (4), where $\lambda[n]=\frac{1}{T_s}\int_{nT_s}^{(n+1)T_s}\lambda(t)\,\mathrm{d}t$ is the photon rate averaged over the $n$-th interval. With DCO-OFDM, at $T_s=\tau_d$, however, the rate varies within a dead time, and the exponent of (7) weights the rate of the preceding interval. The two models thus differ from interval to interval, and an estimator based on (4) incurs a power-dependent range bias [14], [15].

Sensing Performance and ISAC Operating Point

Range CRB

The delay and the reflected amplitude, $\boldsymbol{\theta}=[\tau,\alpha_s]^{\mathsf{T}}$, where $(\cdot)^{\mathsf{T}}$ denotes the transpose, are unknown and estimated jointly, whereas $\lambda_b$ is known, e.g., measured with the LED off. Successive counts are treated as independent, following the per-symbol observations of RF ISAC [4], [5]. This neglects the negative correlation between adjacent counts that the dead time induces.

Since the dead time correlates adjacent counts, (10) is the CRB of the independent-count model rather than a strict bound for the actual counts. Moreover, as $\bar{\mathbf{J}}$ averages the FIM over the data before inversion, (10) lower-bounds the CRB averaged over the data [5].

In the quasi-static special case, $\tilde q[n]=\mu_a(\lambda[n])/N_a$ and, at $T_s=\tau_d$, $N_a\tilde q[n](1-\tilde q[n])=\sigma_a^2(\lambda[n])$. Hence, (9) becomes $[\mathbf{J}]_{ij}=\sum_n\Xi(\lambda[n])\frac{\partial\lambda[n]}{\partial\theta_i}\frac{\partial\lambda[n]}{\partial\theta_j}$, with $\Xi=\mu_a^{\prime2}/\sigma_a^2$ and $\mu_a^\prime=\mathrm{d}\mu_a/\mathrm{d}\lambda=T_s e^{-\lambda}(1-e^{-\lambda})$. Since $\lambda[n]$ averages $\lambda_b+\alpha_s x(t-\tau)$ over the $n$-th interval, $\partial\lambda[n]/\partial\tau=-[\lambda((n+1)T_s)-\lambda(nT_s)]/T_s$, as for (8b). In both models, the delay information thus enters through the change of a rate across the edges of the interval. In the quasi-static model, this is the incident rate $\lambda(t)$, scaled by the slope $\mu_a^\prime$ of the mean response, and in (9), it is the detected-photon rate $r(t)$, cf. (8b).

Since $\mu_a^\prime$ vanishes at $\lambda^\star$ and $\Xi$ decays as $e^{-\epsilon}$ beyond, a saturated sample carries vanishing delay information. For the PD benchmark of the same bare aperture, which observes $\rho_{\mathrm{pd}}\xi x(nT_s-\tau)$ in the thermal noise of variance $\sigma_c^2$, this Gaussian observation yields the classical linear delay bound [17], which is monotone in $P_t$.

Rate-Constrained Operating Point

Communication and sensing share the average optical power $P_t\approx B$. Since $\sigma_s=P_t/\beta$, the electrical SNR is quadratic in the optical power, $\operatorname{SNR}_c=\gamma_cP_t^2$ with $\gamma_c=h_c^2/(\beta^2\sigma_c^2)$. On each data subcarrier, adaptive $M$-QAM at a target bit error rate $\mathrm{BER}_t$ operates at the SNR gap $\Gamma= -\ln(5\mathrm{BER}_t)/1.5$ from capacity [18], so that the achievable rate is

Simulation parameters
$\tau_d = T_s$	10 ns	$N_a$	$\{256, 1024, 4096, 8192\}$
$\Upsilon_{\mathrm{PDE}}$	0.35	$\vartheta_{\mathrm{DCR}}, P_B$	0.5 MHz, 10 nW
$\phi_{\mathrm{AP}}, \phi_{\mathrm{CT}}$	0.75%, 2.5%	$A_{\mathrm{spad}} = A_{\mathrm{pd}}$	1 cm$^2$
$k_{\mathrm{eff}}$	0.92	$\rho_{\mathrm{diff}}, A_{\mathrm{tgt}}$	0.7, 100 cm$^2$
$N, L_{\mathrm{cp}}, L$	64, 8, 30	$\beta, \Psi_{1/2}\ (\mathrm{m})$	3, $60^\circ$ (1)
$d_c, d_s$	3 m, 2 m	$N_0, \rho_{\mathrm{pd}}$	$10^{-22}$ A/$\sqrt{\mathrm{Hz}}$, 0.25 A/W
Wavelength	450 nm	$\kappa_{\mathrm{opt}}, \phi_c = \psi_c$	2.5, $48^\circ$
$\mathrm{BER}_t$	$10^{-3}$	$P_{\max}$	2 W
$\Gamma$ at $\mathrm{BER}_t$	5.5 dB	$\lambda_b$	$8.2 \times10^9\ \mathrm{s}^{-1}$

The rate is strictly increasing in $P_t$, so that a target rate $R_0$ is equivalent to the power floor $P_t^c(R_0) = \sqrt{\frac{\Gamma(2^{R_0(N+L_{\mathrm{cp}})/(NB_{\mathrm{occ}})}-1)}{\gamma_c}}$. Since the ranging precision is not monotone in $P_t$, the operating point is the power that minimizes $\mathrm{CRB}(D)$ over $0 < P_t \leq P_{\max}$, with LED power limit $P_{\max}$, subject to $R_c(P_t) \geq R_0$, which is feasible for $R_0 \leq R_c(P_{\max})$. Let $P_t^s$ denote the minimizer without the rate constraint, the sensing-optimal power. If the CRB does not decrease on $[P_t^s, P_{\max}]$, the operating point over the feasible powers $[P_t^c(R_0), P_{\max}]$ is

i.e., the system is sensing-limited when $R_0 \leq R_c(P_t^s)$ and rate-limited otherwise. For a linear receiver or a non-saturating target, the sensing quality is monotone in the power, so that (13) collapses to the maximum power.

Numerical Results

The simulation parameters are listed in Table I and follow [2, 6, 7, 9]. The user’s PD is 2 m below the BS and 2.24 m off its axis, i.e., $\phi_c = \psi_c = 48^\circ$, and for $\beta= 3$, the clipping probability $Q(\beta) \approx1.3 \times10^{-3}$ is negligible. Figure 2(b) and Figure 3 are reported versus the received optical power $P_{\mathrm{rx}} = \zeta P_t$, or equivalently the mean normalized photon rate $\epsilon= (\lambda_b + C_sP_{\mathrm{rx}})\tau_d/N_a$. All markers are obtained by Monte Carlo (MC) simulation. For the SPAD, the photon arrivals at each pixel, including those during the CP of each block, are drawn from a Poisson process of rate $\lambda(t)/N_a$, a pixel detects a photon only if no photon arrived in the preceding $\tau_d$ seconds, and $y[n]$ is the number of pixels that detect a photon in the $n$-th interval. Each marker in Figure 2(b)–Figure 4 is the root-mean-square error (RMSE) of the range estimate, averaged over 200 frames with new data. For each frame, the ML estimator maximizes the log-likelihood of the binomial counts (6) of its $L$ blocks,

where $y_l[n]$ is the count in the $n$-th interval of the $l$-th block and $\tilde{q}_l[n]$ is its firing probability for the parameters $\theta$, given by (7) or, for the estimator based on the quasi-static model, by

$\mu_a(\lambda[n])/N_a$ from (4). For the PD, the ML estimator uses its samples $\rho_{\mathrm{pd}}\mathcal{X}(nT_s-\tau)$ in Gaussian noise, and its RMSE is averaged over 300 frames. The analytical curves in Figure 2(b)–Figure 4 are averaged over $10^4$ blocks of 16-QAM data.

Figure 2 compares the two count models in (4) and (7) at $T_s = \tau_d$ and shows their impact on ranging. Figure 2(a) shows the mean count of samples $n = 20$ to 44 of one data block at $P_{\mathrm{rx}} = 0.41\ \mu\mathrm{W}$, i.e., in saturation ($\epsilon\approx3.4$). For this block, the pile-up model gives $N_a\tilde{q}[n]$ from (7), the quasi-static model gives $\mu_a(\lambda[n])$ from (4), and the MC markers average $y[n]$ over 4000 photon realizations. The pile-up model follows the MC, whereas the quasi-static model misses the sample-to-sample variations. The reason is that, at $T_s = \tau_d$, the count in an interval depends on the photons of the preceding interval, whereas the quasi-static model assumes the same rate in both. After a dark interval, few pixels are blind, so that a bright interval yields more counts than (4) predicts. For instance, after the intensity dip at $n = 27$, 723 pixels fire at $n = 28$, whereas (4) never exceeds its maximum $N_a/e \approx377$.

Figure 2(b) shows how this mismatch affects ranging, from the unsaturated regime ($\epsilon\approx0.1$) to deep saturation ($\epsilon\approx40$). From the same counts, $\theta$ is estimated by the ML estimator (13) based on (7) and by the one based on (4), and the markers give the resulting range RMSE. The estimator based on (7) follows $\sqrt{\mathrm{CRB}}$ from (10) except at the highest powers. The RMSE of the estimator based on (4), in contrast, is 0.08 to

1.2 m, up to 400 times larger, and does not follow the trend of the bound, i.e., it grows with the received power in the unsaturated regime, where the bound decreases. The bias of this estimator, i.e., the mean of its range error, is more than 98% of this RMSE, which reveals the power-dependent range bias predicted in Section II.

Figure 3 shows $\sqrt{\mathrm{CRB}}$ from (10) for four array sizes, the quasi-static bound, i.e., (10) evaluated with the quasi-static model (4), $\tilde{q}[n] = \mu_a(\lambda[n])/N_a$, for $N_a = 1024$, and the RMSE of the ML estimator (13) based on (7). The estimator is unbiased, i.e., its mean range error is zero within about two standard errors of the MC average, and up to $\epsilon\approx10$ its RMSE exceeds $\sqrt{\mathrm{CRB}}$ by 9% on average. The minimum of the bound improves roughly as $N_a^{-1/2}$, from 5.0 to 0.8 mm for $N_a = 256$ to 8192, and lies at $P_{\mathrm{rx}} \approx0.11$, 0.45, 1.6, and 3.6 $\mu\mathrm{W}$, respectively, i.e., near $\epsilon\approx3.5$ for all array sizes, beyond which the bound increases. Deep in saturation, the bound (10) becomes loose. The quasi-static bound is up to 2.6 times larger at its blind spot $\epsilon\approx1$, where $\mu'_a = 0$, and within about 25% of (10) from the optimum onward. This bound holds only if the counts follow (4); for the actual counts, the estimator based on (4) is biased, as discussed for Figure 2(b).

Figure 4 shows the resulting root CRB-rate region [5], traced by sweeping $P_t$ up to $P_{\max}$, together with the RMSE of the corresponding ML estimators. The rate grows monotonically with the power, whereas the CCR range CRB attains its minimum at $P_t^s \approx0.22$ W, i.e., at $R_c \approx17$ Mbps and $\sqrt{\mathrm{CRB}} \approx2.2$ mm. Beyond this point, a higher power increases the rate at the cost of ranging accuracy, and at $P_{\max}$, where $R_c \approx203$ Mbps, the bound is about 3 times larger than at $P_t^s$. The CRB thus does not decrease beyond $P_t^s$, as (12) requires. The diffuse target remains unsaturated, so that its CRB decreases monotonically and $P_t^* = P_{\max}$. The PD bound is monotone and is attained within 10%. The SPAD is about 17 times more precise at $P_t^s$, and the two bounds cross only near 160 Mbps, beyond which the PD is more precise. Hence, the trade-off is a receiver-physics effect.

Conclusion

We analyzed ranging with a SPAD array in DCO-OFDM O-ISAC. Sampling the array once per dead time, we modeled the array count in each sampling interval as a binomial random variable that accounts for the intensity variation within a dead time, derived the range CRB, and validated it by MC simulation of the photon counting. At this sampling rate, the SPAD count model of OWC, which assumes a constant photon rate over a dead time, causes a power-dependent range bias. Moreover, saturation makes the ranging accuracy of the SPAD array best at a sensing-optimal transmit power, whereas a linear receiver benefits from more power. Under a rate constraint, the optimal transmit power is then the larger of this power and the minimum power required by the target rate.