© The Authors 2025

1 Introduction

Photometry is an important technique in astrophysics used to measure the intensity or brightness of celestial objects, such as solar system bodies, stars, and galaxies (Sterken & Manfroid 1992; Elliot & Lub 1995; Budding & Demircan 2007). In recent decades, the advent of space missions dedicated to high-precision photometric monitoring, such as CoRoT (Auvergne et al. 2009), Kepler (Borucki et al. 2010), its extended K2 mission (Howell et al. 2014), and TESS (Ricker et al. 2015), have transformed the field by enabling unprecedented photometric precision and continuity in time series. These capabilities are critical for applications such as detecting and characterizing exoplanets, probing stellar oscillations, and studying long-term stellar variability. In the near future, the PLATO mission (Rauer et al. 2014) and, farther ahead, the proposed HAYDN project (Miglio et al. 2021)¹ , aim to extend this legacy by delivering photometry with even higher precision, longer baselines, and improved characterization of systematic effects. Such missions require highly optimized photometric extraction methods, as their precision goals approach the fundamental limits imposed by photon statistics and the stability of the instrument.

Photometric measurements are essential for various fields of study, from detecting and characterizing exoplanets via transits to monitoring variable stars and supernovae; these measurements can help unravel stellar evolutions or cosmological distances (Gallaway 2020; Kitchin 2020). As light from celestial sources often gets contaminated by background noise, whether from the sky, adjacent objects, or instrumental interference, accurate photometry requires techniques that can effectively isolate the target’s light from this background (Chromey 2016; Sterken et al. 2011). Traditional methods perform these estimations in a sequential two-step manner (Mighell 1999), first estimating the background and then the source flux (see, e.g., DAOPHOT in Stetson (1987), DOPHOT in Schechter et al. (1993), MATPHOT in Mighell (2005), or PSFEx in Bertin (2011) using a parametric approach, and Conseil et al. (2011) for a Bayesian approach, among others). This two-step approach has been demonstrated by Vicuña et al. (2024) to introduce biases and limit the attainable precision, proposing instead a novel methodology for estimating background and source flux simultaneously for aperture photometry, applicable to isolated sources. For this paper we adopted that methodology, but we now address the issue of the impact in photometric precision of the choice of aperture when performing point spread function (PSF) fitting in crowded fields.

Classical photometric techniques have long relied on the optimization of the signal-to-noise ratio (S/N) as a key determinant of accuracy, particularly when using charged-coupled devices (CCD) for point-like source observations. Traditional methods, such as those outlined by Howell (1989); Newberry (1991); Naylor (1998), have made significant contributions to improving photometric precision by carefully managing noise sources, including sky background, dark current, and read-out noise. Howell’s introduction of the CCD equation (see also Merline & Howell (1995), their Eq. (21)), which accounts for various noise terms beyond simple Poisson statistics, and Newberry’s thorough treatment of noise propagation in sky-subtracted data have provided astronomers with valuable tools for improving photometric precision. These methods often involve using a small aperture to deal with crowded fields, where aperture photometry is unfeasible, while maximizing the S/N. As we shall demonstrate below, this approach is not always a good benchmark indicator, and can severely limit photometry precision, particularly when there is some flexibility in choosing the PSF fitting radius.

Modern photometry methods, based on the Cramér-Rao lower bound (CRLB), have been employed as a key indicator of photometric precision Lee & van Altena (1983); Adorf (1996). More importantly, estimators have been proposed that achieve this bound with extremely high accuracy, where the precision is a direct function of the number of observations and, consequently, the aperture radius. This represents a paradigm shift in photometric techniques as it indicates that to achieve the highest possible precision (Mendez et al. 2013, 2014; Espinosa et al. 2018) it is necessary to utilize the maximum amount of information available, which implies using the entire usable image around the target for photometric analysis – if possible – rather than limiting the data to pre-selected pixels or regions. This approach challenges traditional methods and underscores the importance of accessing and incorporating all available relevant pixel data to refine flux and background estimations, leading to more precise and accurate results in astronomical photometry. For example, Perryman et al. (1989) suggested that the precision of flux measurements is improved by better accounting for both the background and the source’s flux simultaneously, leading to more accurate and efficient photometric pipelines.

In Vicuña et al. (2024), they show that the classical approach of determining first the background in an area surrounding the source and subtracting that estimation from the source to quantify the flux of the target is sub-optimal in comparison with a joint simultaneous estimation of these quantities. This was demonstrated from theoretical grounds and based on extensive numerical simulations under realistic observational scenarios. They also tested their methodology on very high-precision photometry based on TESS satellite data on isolated sources, confirming their theoretical predictions. To our knowledge, this is the first completely self-consistent treatment of this problem, on solid statistical grounds, which gives us theoretical assurances of optimality (in the CRLB sense) as well as unbiasedness on the derived quantities of interest, namely, source flux and background.

Given the results presented in the previous paragraph, it becomes essential to examine how the selection of the aperture radius during PSF fitting influences the final photometric precision. In many practical scenarios, such as crowded fields or varying seeing conditions, the choice of aperture radius cannot be arbitrary and often plays a pivotal role in balancing precision and feasibility. Despite its importance, this aspect is commonly guided by heuristics, which may not always reflect the true limits of photometric accuracy.

In this paper, we propose a novel methodology that reframes the problem of aperture selection from an information-theoretic perspective. Instead of relying on S/N, we compute the CRLB for the joint estimation of source flux and background level, thereby quantifying the theoretical precision limit associated with each aperture choice. The Fisher information, which underpins the CRLB, provides a principled measure of the data’s informativeness and allows a rigorous assessment of the trade-offs involved in aperture selection. The main contributions of this work are summarized below.

• We introduce an information-theoretic approach to photometry based on the Fisher information matrix. This approach represents a fundamental paradigm shift in how we approach photometry, moving beyond S/N to a more robust statistical foundation that enables higher precision and flexibility in observational strategies while, at the same time, offering guarantees of precision and unbiased results in a quantifiable way.

• By basing our approach on this statistical framework, we are able to quantify in a precise manner how the choice of a given PSF fitting radius impacts on the accuracy of flux measurements.

• We demonstrate that the proposed method consistently out-performs traditional S/N-based approaches, especially for larger apertures where bias and underestimation of uncertainty are common.

• Our results provide a data-driven tool to inform PSF aperture selection, offering both improved accuracy and flexibility for a wide range of observational constraints.

The paper is organized as follows. In Sect. 2 we revisit the fundamental results of the astrometric model for characterizing point sources on a CCD. In Sect. 3 we formalize the optimal aperture analysis and present the main results. In Sect. 4 we show numerical results analyzing discrepancies between the performance based on Fisher information and some classical practical estimators. Finally, in Sect. 5 we offer our conclusions and final remarks on our findings.

2 Preliminaries

We present a summary of the results obtained by Vicuña et al. (2024) in the context of the joint source and background flux. For completeness, the subsequent section will describe the model in some detail, with a focus on emphasizing the most relevant results that will be critical for the following section.

2.1 Statistical parametric model for photometry

The primary focus is the inference of the flux and background of a point source. In our simplified but realistic model, the source is characterized by two scalar parameters: the object’s position x_c ∈ ℝ within the array, and its aggregated intensity², denoted as F ∈ ℝ⁺. These two parameters define a probability distribution ${\mu }_{x_c,F}$ over an observation space, denoted as 𝕏. Formally, a point source specified by the pair (x_c, F) generates a nominal intensity profile on a photon-integrating device (PID), such as a CCD, which can be expressed as $$F_{x_c}\left(x\right)=F\cdot \varphi \left(x-x_c,\sigma \right),$$(1)

where ϕ(x − x_c, σ) represents the one-dimensional normalized point spread function (PSF), with σ serving as a generic parameter that defines the width (or spread) of the light distribution on the detector. This parameter is generally influenced by the wavelength and the observing site’s quality (Mendez et al. 2013, 2014).

The profile in Eq. (1) is not observable directly; rather, it is subject to the following three sources of perturbations:

• an additive background noise, which accounts for photon emissions from the open (diffuse) sky f_s, as well as noise introduced by the instrument itself, including read-out noise (RON) and dark current (D) (see, e.g., Janesick et al. 1987; Howell 2006; Janesick 2007; McLean 2008), modeled by B in Eq. (2)³: $$B=f_s\cdot \Delta x+D+RO{N}^2;$$(2)

• an uncertainty between the aggregated intensity (comprising the nominal object brightness and the background) and the actual measurements. Independent random variables model this uncertainty following a Poisson distribution (an analysis similar to the one that follows, but for Gaussian noise typical of photographic plates, was done by Lee & van Altena (1983), including the calculation of the Fisher matrix, see their Eq. (8));

• the spatial quantization process associated with the pixel resolution of the PID, as specified in Eqs. (3) and (4).

By modeling these effects, we have a countable collection of independent and non-identically distributed random variables (observations or counts) {I_i: i ∈ ℤ}, where I_i ~ Poisson(λ_i(F, B)), driven by the expected intensity at each pixel element i⁴, given by $${\lambda }_i\left(F,B\right)\equiv \text{𝔼}\left(I_i\right)=F\cdot {𝑔}_i\left(x_c\right)+B_i,\forall i\in \text{\hspace{0.17em}ℤ,}$$(3) $${𝑔}_i\left(x_c\right)\equiv {\displaystyle {\int }_{x_i-\Delta x/2}^{x_i+\Delta x/2}\varphi \left(x-x_c,\sigma \right)dx,\forall i\in \text{\hspace{0.17em}ℤ,}}$$(4) $$\varphi \left(x,\sigma \right)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(\frac{-x_2}{2{\sigma }^2}\right),$$(5)

where 𝔼 {} is the expectation value of the argument and {x_i: i ∈ ℤ} denotes the standard uniform quantization of the real line array with resolution ∆x > 0, i.e., x_i+1 − x_i = ∆x for all i ∈ ℤ. The function ɡ_i(x_c) in Eq. (4) describes the distribution of (source) flux on the detector or, equivalently, the image profile across pixels, integrated over pixel i.

In practice, the PID has a finite collection of measurement elements (or pixels). Then we consider a finite set 𝓝 = {1,…, n} of indices over the pixels (for simplicity, in what follows we assume that n is odd) of a CCD imaging a point source whose relative position x_c is known, and n = |𝓝|. Another usual convention is to characterize the spread (width) of the PSF by its full width at half maximum (FWHM), which relates to the parameter σ through $$FWHM=2\sqrt{2\ln 2}\sigma .$$(6)

Adopting this parametric model, the likelihood function of I = {I₁,…, I_n}, given our target parameter vector (F, B), is given by $$L\left(\mathbf{\text{I}};F,B,\mathcal{N}\right)={\displaystyle \prod _{i\in \mathcal{N}}{P}_{{\lambda }_i\left(F,B\right)}\left(I_i\right),}$$(7) $$P_{{\lambda }_i(F,B)}\left(I\right)=\frac{e^{-{\lambda }_i(F,B)}\cdot {\lambda }_i^I(F,B)}{I!},$$(8)

where $P_{\lambda }\left(x\right)=\frac{e^{-\lambda }\cdot {\lambda }^x}{x!}$ denotes the probability mass function (PMF) of the Poisson law evaluated at the realization of I = x (see, e.g., Kay (1993)). Equation (7) depends on the selected set of pixels. However, due to the independence of the observations, the likelihood will only change based on the observations chosen, not the order in which they are selected. In the following sections, we consider subsets 𝒥_u ⊂ 𝓝 of pixels of the form ${\mathcal{J}}_u=\left\{n_{x_c}-n_u,\dots ,n_{x_c},\dots ,n_{x_c}+n_u\right\}$, where n_u is the number of full pixels within the aperture radius we have, and $n_{x_c}$ is the number of the pixel associated with the position of x_c (typically n/2 if the array is centered). Finally, x_c is assumed to be known and the estimation task reduces to finding the parameter estimator $\widehat{\mathit{\theta }}=\left(\widehat{F},\widehat{B}\right)$, which is the inference of the underlying parameters (F, B) from I. Formally, the inference task consists in defining a regression rule τ(I) : ℕⁿ → Θ = ℝ⁺ × ℝ⁺, such that $\left(\widehat{F},\widehat{B}\right)\equiv \tau \left(\mathbf{I}\right)$.

2.2 The Cramér-Rao lower bound (CRLB)

The well-known CRLB offers a performance bound on the variance estimation error) of the family of unbiased estimators⁵. Here, we revisit the multidimensional⁶ version of this result.

Theorem 1. Rao (1945); Cramér (1946) Let 𝓝 be a finite set of integers with |𝓝| = n; {I_i}_i∈𝓝 be a collection of independent observations, whose likelihood function L(·; θ, 𝓝) is induced by a parameter vector θ = (θ₁,…, θ_m) ∈ Θ, m ∈ ℕ over a parameter space Θ (typically Θ = ℝ^m), such that the following regularity condition is satisfied: $$\text{E}\left\{\frac{\partial \ln L\left(\mathbf{I};\mathit{\theta },\mathcal{N}\right)}{\partial {\theta }_i}\right\}=0,\forall i\in \left\{1,\dots ,m\right\},\forall \mathit{\theta }\in \mathbf{\Theta }.$$(9)

(i) Then, any unbiased estimator $\widehat{\mathit{\theta }}$ of θ, given by a regression rule τ : ℕⁿ → Θ, $\widehat{\mathit{\theta }}\equiv \tau \left(\mathbf{I}\right)$, has a covariance matrix $K_{\widehat{\mathit{\theta }}}\equiv \text{E}\left\{\left(\widehat{\mathit{\theta }}-\mathit{\theta }\right)\cdot {\left(\widehat{\mathit{\theta }}-\mathit{\theta }\right)}^{†}\right\}$ that satisfies⁷ $$K_{{\mathit{\theta }}^{\prime }}-I_{\mathit{\theta }}^{-1}\ge \mathbf{0.}$$(10)

I_θ ∈ ℳ^m×m denotes the Fisher information matrix, whose components are defined by ∀i, j ∈ {1,…, m}, $${\left[I_{\theta }\right]}_{\left(i,j\right)}\equiv \text{𝔼}\left\{\frac{{\partial }^2\ln L\left(I;\mathit{\theta },\mathcal{N}\right)}{\partial {\theta }_i\partial {\theta }_j}\right\}.$$(11)

(ii) Furthermore, if there exists a function h : ℝⁿ → Θ such that $$\frac{\partial \ln L\left(I;\mathit{\theta },\mathcal{N}\right)}{\partial {\theta }_i}={\left[I_{\theta }\left(\mathbf{\text{h}}\left(\mathbf{\text{I}}\right)+\mathit{\theta }\right)\right]}_i,\forall i\in \left\{1,\dots ,m\right\},$$(12)

then the minimal variance unbiased estimator (MVUE) is given by $\widehat{\mathit{\theta }}=\mathbf{\text{h}}\left(\mathbf{\text{I}}\right)$, and its optimal covariance matrix is $I_{\mathit{\theta }}^{-1}$.

From Theorem 1, any unbiased estimator ${\widehat{\theta }}_i$ of θ_i satisfies that (from Eq. (10)) $$Var\left({\widehat{\theta }}_i\right)={\left[K_{\widehat{\theta }}\right]}_{\left(i,i\right)}\ge {\left[I_{\theta }^{-1}\right]}_{\left(i,i\right)},\forall i\in \left\{2,\dots ,m\right\},$$(13)

and therefore ${\left[I_{\mathit{\theta }}^{-1}\right]}_{\left(i,i\right)}$ is a lower bound for the mean squared error (MSE) of ${\widehat{\theta }}_i$⁸.

2.3 CRLB for the joint estimation of source flux and background

First, we revisit a result in which the CRLB is derived for the joint source flux and background estimation. We have the following result:

Lemma 1. (Vicuña et al. 2024, Lemma 3.1) If the astrometry x_c ∈ ℝ is fixed and known, and we want to estimate the pair (F, B) from I ~ L(I; F, B, 𝓝) in (7), then the Fisher information matrix in Eq. (10) is given by $$I_{\theta }=\left[\begin{array}{cc}{I}_{1,1}\left(\mathcal{N}\right)& I_{1,2}\left(\mathcal{N}\right)\\ I_{2,1}\left(\mathcal{N}\right)& I_{2,2}\left(\mathcal{N}\right)\end{array}\right]\equiv \left[\begin{array}{cc}{\displaystyle \sum _{i\in \mathcal{N}}\frac{{𝑔}_i^2\left(x_c\right)}{{\lambda }_i\left(F,B\right)}}& {\displaystyle \sum _{i\in \mathcal{N}}\frac{{𝑔}_i\left(x_c\right)}{{\lambda }_i\left(F,B\right)}}\\ {\displaystyle \sum _{i\in \mathcal{N}}\frac{{𝑔}_i\left(x_c\right)}{{\lambda }_i\left(F,B\right)}}& {\displaystyle \sum _{i\in \mathcal{N}}\frac{1}{{\lambda }_i\left(F,B\right)}}\end{array}\right].$$(14)

Consequently, the following two-dimensional CRLB for any pair of unbiased estimates $\widehat{F}$ and $\widehat{B}$ can be respectively defined as $$Var\left(\widehat{F}\right)\ge {\sigma }_F^2={\left[I_{\theta }^{-1}\right]}_{\left(1,1\right)}=\frac{I_{2,2}\left(\mathcal{N}\right)}{I_{1,1}\left(\mathcal{N}\right)\cdot I_{2,2}\left(\mathcal{N}\right)-I_{1,2}^2\left(\mathcal{N}\right)},$$(15) $$Var\left(\widehat{B}\right)\ge {\sigma }_B^2={\left[I_{\theta }^{-1}\right]}_{\left(2,2\right)}=\frac{I_{1,1}\left(\mathcal{N}\right)}{I_{1,1}\left(\mathcal{N}\right)\cdot I_{2,2}\left(\mathcal{N}\right)-I_{1,2}^2\left(\mathcal{N}\right)}.$$(16)

We present the equivalent joint Fisher information for the flux and background as $$I_F\left(\mathcal{N}\right)=I_{1,1}\left(\mathcal{N}\right)-\frac{I_{1,2}^2\left(\mathcal{N}\right)}{I_{2,2}\left(\mathcal{N}\right)},$$(17) $$I_B\left(\mathcal{N}\right)=I_{2,2}\left(\mathcal{N}\right)-\frac{I_{1,2}^2\left(\mathcal{N}\right)}{I_{2,2}\left(\mathcal{N}\right)}.$$(18)

The statistical interaction (correlation) between components is captured by the off-diagonal elements of the Fisher information matrix in Eqs. (15) and (16). Importantly, in the special case that I_1,2= 0, the estimates become decoupled and one could consider that the joint estimation task reduces to two isolated 1D estimation problems.

2.4 Model assumptions

In this subsection, we outline the assumptions that will aid in the analysis of the photometric errors induced by a given aperture (PSF fitting radius) choice. These assumptions are based on the fact that the observations have good coverage and resolution when capturing the images.

2.4.1 Good spatial coverage of the PSF on the detector

A basic assumption here is that we have a good coverage of the object of interest, in the sense that for a given position x_c, then this condition implies that the set of all n available pixels is such that $${\displaystyle \sum _{i\in \mathcal{N}}{𝑔}_i\left(x_c\right)}\approx {\displaystyle {\int }_{-\infty }^{\infty }\theta \left(x-x_c,\sigma \right)dx=1.}$$(19)

2.4.2 High-resolution regime

This approximation considers that the imaging device allows for a sufficiently high resolution, such that ∆x ≪ σ. In such a scenario, the following approximation can be considered: $${𝑔}_i\left(x_c\right)\approx \varphi \left(x_i-x_c,\sigma \right)\cdot \Delta x.$$(20)

To account for the actual coverage of the chosen aperture, we recognize that the following function from Mendez et al. (2013, 2014) allows us to do that: $$P\left(u\right)\equiv \text{erf}\left(\frac{u}{\sqrt{2}\sigma }\right)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^{\frac{u}{\sqrt{2}\sigma }}{e}^{-{\upsilon }^2}}d\upsilon .$$(21)

Here u is the aperture radius measured in the same units as σ (e.g., in arcseconds). The quantity P(u) represents the fraction of the signal that is captured within the chosen aperture, which is exactly the interpretation of the more general term ${\Sigma }_{i\in {\mathcal{J}}_u}{g}_i\left(x_c\right)$.

2.4.3 Approximations for $g_i^p\left(x_c\right)$

We note that the expressions associated with Fisher information in Eq. (14) involve discrete sums; however, we can calculate approximations to obtain continuous representations as a function of the aperture radius u. The following proposition provides closed-form expressions for the sum of ɡ_i(x_c) across the pixels.

Proposition 1. With the high-resolution assumption, we have the following approximations for the sum of ɡ_i(x_c); in particular, we have that $${\displaystyle \sum _{i\in {𝒥}_u}{𝑔}_i^p\left(x_c\right)\approx \frac{1}{\sqrt{p}}\cdot {\left(\frac{1}{\sqrt{2\pi }}\right)}^{p-1}\cdot {\left(\frac{\Delta x}{\sigma }\right)}^{p-1}\cdot P\left(u\sqrt{p}\right),}$$(22)

for p ≥ 1. Analogously (also for p ≥ 1), it can be found that $${\displaystyle \sum _{i\in {𝒥}_u}\frac{1}{{𝑔}_i^p\left(x_c\right)}\approx \frac{{\left(\sqrt{2\pi }\right)}^{p+1}}{\sqrt{p}}{\left(\frac{\sigma }{\Delta x}\right)}^{p+1}\text{erfi}\left(\frac{u}{\sqrt{2}\sigma }\sqrt{p}\right),}$$(23)

where erfi (x) is defined through the erf function as $$\begin{array}{l}\text{erfi\hspace{0.17em}}\left(x\right)=-i\text{\hspace{0.17em}erf\hspace{0.17em}}\left(ix\right),\\ andi=\sqrt{-1}correspondstotheima𝑔inaryunit.\end{array}$$(24)

These results are demonstrated in Appendix A. The approximations proposed here are useful as they allow us to derive analytical closed-form expressions for the Fisher information terms and to perform an in-depth analysis, which is presented in the following section.

2.5 Jitter and position estimation in the PSF model

A potential source of photometric performance degradation, in particular in space-based photometry, is satellite jitter, which are small pointing errors occurring during an exposure. These fluctuations modify the effective distribution of photons on the detector, altering the shape and localization of the observed PSF. In this work, we assume that the source centroid x_c is known with high precision throughout the exposure. However, if the effect of jitter is to be included explicitly, the appropriate treatment is to consider a joint estimation problem where both F and x_c are unknown parameters. In such case, the Fisher information matrix includes cross-terms that quantify the coupling between photometry and astrometry, and marginalizing over x_c reduces the information available for estimating the flux. This formalism has been studied in detail by Mendez et al. (2014), who conducted a systematic Fisher information analysis to characterize how uncertainty in the source position xc degrades the precision of photometric estimates, thereby quantifying the impact of any perturbation in source localization, such as those induced by satellite jitter.

3 Aperture (PSF fitting radius) analysis

Photometry plays a crucial role in accurately estimating properties of celestial objects, where the selection of pixels as a function of aperture radius is essential. Classical PSF fitting studies in this field have primarily focused on maximizing the S/N to enhance measurement precision. In this section we propose a new paradigm, focusing instead on the Fisher matrix and in how it can be interpreted as an object containing the necessary information (within a certain quantifiable variance) for the task of aperture selection. First, we revisit the formal definition of S/N, and based on this, we compare the aperture selection following the classic S/N-based scheme with our new CRLB-based approach.

3.1 S/N as a function of the aperture radius

The S/N is a common metric used to assess the quality of various measurements performed on a celestial object. We revisit the analytical expression obtained from Mendez et al. (2013) for the S/N and then derive properties that will be useful in the following sections. For the case of the Gaussian function PSF, it makes sense to perform an integration of the PSF that is symmetrical with respect to the center of the source, centered at x_c. On the other hand, the total noise, N, has contributions from the read-out noise of the detector, the noise from the sky, and the noise from the source itself, all of which are assumed to follow Poisson statistics. While we are concerned here with crowded fields, which require a small aperture, we assume, within that selected aperture, that the influence of nearby objects is negligible. Integrating what was previously mentioned, the S/N can be written as a function of the aperture radius (see Gilliland (1992) and Mendez et al. (2013, Eq. (28)). More specifically, we have the following expression: $$S/N\left(u\right)=\frac{P\left(\sqrt{2}\sigma u\right)F}{\sqrt{P\left(\sqrt{2}\sigma u\right)F+4{\sigma }^{2}u\left(B-D\right)}}.$$(25)

It is interesting to note that this function has a maximum with respect to u. This is because as the aperture radius increases, the noise contribution also increases since the flux eventually saturates and there are more pixels contributing with background and noise. The following proposition summarizes the concavity property of the S/N.

Proposition 2. The function S /N(u) has a maximum value $u_{S/N}^{*}$, moreover, this optimal aperture radius is given as the solution to the following equation: $$\begin{array}{l}F\frac{\partial P\left(\sqrt{2}\sigma u_{S/N}^{*}\right)}{\partial u}P\left(\sqrt{2}\sigma u_{S/N}^{*}\right)\hfill \\ =-\left[2\sqrt{2}\sigma u_{S/N}^{*}\frac{\partial P\left(\sqrt{2}\sigma u_{S/N}^{*}\right)}{\partial u}-P\left(\sqrt{2}\sigma u_{S/N}^{*}\right)\right]2\sqrt{2}\sigma \left(B-D\right).\hfill \end{array}$$(26)

The proof is presented in Appendix B. Proposition 2 tells us that there is a criterion for selecting the aperture radius for photometry based on maximizing the S/N (Howell 1989; Newberry 1991).

3.2 Joint Fisher information as a function of the aperture radius

We formalize the concept of photometric uncertainty for a given aperture, based on Fisher information, recognizing its role as an optimal performance indicator for unbiased estimators concerning their variance in the context of the CRLB.

The Fisher information plays a crucial role as it acts as a fundamental limit for the joint estimation problem as a function of the number of pixels. In particular, we see that there is a direct connection between Fisher information and the aperture radius, which motivates a systematic analysis regarding the impact of the selected aperture on the associated photometric precision. We restrict our analysis to aperture-like sets of indices, i.e., sets of the form ${\mathcal{J}}_u=\left\{n_{x_c}-n_u,\dots ,n_{x_c},\dots ,n_{x_c}+n_u\right\}$, where n_u is the number of full pixels within the aperture radius.

First of all, thanks to the additive structure in Eq. (14), we see that the Fisher information in the joint photometric case is given, for an arbitrary set of indices J_u, by $$I{{\alpha }^{★}}_{}\left({𝒥}_u\right)=\left[\begin{array}{cc}{I}_{1,1}\left({𝒥}_u\right)& I_{1,2}\left({𝒥}_u\right)\\ I_{2,1}\left({𝒥}_u\right)& I_{2,2}\left({𝒥}_u\right)\end{array}\right]\equiv \left[\begin{array}{cc}{\displaystyle \sum }_{i\in {𝒥}_u}\frac{{𝑔}_i^2\left(x_c\right)}{{\lambda }_i\left(F,B\right)}& {\displaystyle \sum }_{i\in {𝒥}_u}\frac{{𝑔}_i\left(x_c\right)}{{\lambda }_i\left(F,B\right)}\\ {\displaystyle \sum }_{i\in {𝒥}_u}\frac{{𝑔}_i\left(x_c\right)}{{\lambda }_i\left(F,B\right)}& {\displaystyle \sum }_{i\in {𝒥}_u}\frac{1}{{\lambda }_i\left(F,B\right)}\end{array}\right].$$(27)

Given a total of |𝓝| = n pixels, the flux is dominant around the center of the aperture mask, whereas in pixels outside the aperture radius (or far from the center), the background becomes dominant. By focusing on the flux-dominant pixels at the center, we can estimate the accuracy of photometric measurements for a given PSF fitting radius. Then we have the following expressions for the Fisher information: $$\begin{array}{c}{I}_F\left({𝒥}_u\right)=I_{1,1}\left({𝒥}_u\right)-\frac{I_{1,2}^2\left({𝒥}_u\right)}{I_{2,2}\left({𝒥}_u\right)}\\ ={\displaystyle \sum _{i\in {𝒥}_u}\frac{{𝑔}_i^2\left(x_c\right)}{{\lambda }_i\left(F,B\right)}-\frac{{\left({\displaystyle {\sum }_{i\in {𝒥}_u}\frac{{𝑔}_i\left(x_c\right)}{{\lambda }_i\left(F,B\right)}}\right)}^2}{{\displaystyle {\sum }_{i\in {𝒥}_u}\frac{1}{{\lambda }_i\left(F,B\right)}}},}\end{array}$$(28) $$\begin{array}{c}{I}_B\left({𝒥}_u\right)=I_{2,2}\left({𝒥}_u\right)-\frac{I_{1,2}^2\left({𝒥}_u\right)}{I_{1,1}\left({𝒥}_u\right)}\\ ={\displaystyle \sum _{i\in {𝒥}_u}\frac{1}{{\lambda }_i\left(F,B\right)}-\frac{{\left({\displaystyle {\sum }_{i\in {𝒥}_u}\frac{{𝑔}_i\left(x_c\right)}{{\lambda }_i\left(F,B\right)}}\right)}^2}{{\displaystyle {\sum }_{i\in {𝒥}_u}\frac{{𝑔}_i^2\left(x_c\right)}{{\lambda }_i\left(F,B\right)}}}.}\end{array}$$(29)

We propose a novel information criterion, noting that the Fisher information is a direct quantitative indicator of the maximum achievable precision for a given set of pixels (or aperture). It is well known that having all pixels available for the estimation process improves performance, achieving the best photometric precision (Vicuña et al. 2024) in the CRLB sense, this is the case of isolated targets. In crowded fields, however, only a subset of pixels is chosen to avoid contamination from nearby sources. As a result, the challenge lies in finding an aperture that captures the maximum amount of useful flux while avoiding interference from nearby sources, ensuring accurate photometric and astrometric estimations within these practical limitations. This motivates the need to explore how to select a subset of pixels such that, with an arbitrarily small error, that particular pixel selection ensures that the Fisher information is close to that obtained when using all pixels. This approach optimizes resource use while maintaining a high level of precision in the photometric estimations.

We propose the following methodology. Formally, given a discrepancy level δ > 0, we seek to find the minimum radius u^* with its corresponding subset ${\mathcal{J}}_{u^{*}}\subseteq N$ such that $$\frac{I_F\left(\mathcal{N}\right)-I_F\left({\mathcal{J}}_{u^{★}}\right)}{I_F\left(\mathcal{N}\right)}\le \delta$$(30)

or, equivalently, $$\frac{\underset{u\in \left\{1,\dots ,\frac{n-1}{2}\right\}}{\inf }{I}_F\left({𝒥}_u\right)}{I_F\left(\mathcal{N}\right)}\ge 1-\delta .$$(31)

While Eq. (30) can be implemented numerically, it is a discrete problem since it depends on the pixels and not on the aperture radius. Therefore, it is convenient to make an approximation to reframe the optimization problem and transform it into a continuous optimization problem. To achieve this, we recall that this analysis is a function of the set J_u, and we consider the assumptions from Sect. 2.4, in particular, the high-resolution regime and the approximations for ɡ_i(x_c). Additionally, we note that given N pixels, the pixels closest to the central pixel contain a larger portion of the flux (F ≫ B). Imposing these conditions for the case of photometry and background estimation is equivalent to making a zero-order approximation of the joint Fisher information. Integrating all of the above, we have the following result.

Theorem 2. Under the high-resolution scenario and the approximations for ɡ_i(x_c) in Sect. 2.4, the Fisher information for the flux and background can be written as $$I_F\left({𝒥}_u\right)=\frac{1}{F}\left(P\left(u\right)-\frac{2u^2}{\pi {\sigma }^2\text{\hspace{0.17em}erfi}\left(\frac{u}{\sqrt{2}\sigma }\right)}\right),$$(32) $$I_B\left({𝒥}_u\right)=\frac{1}{F}\left({\left(\frac{\sqrt{2}\pi \sigma }{\Delta x}\right)}^2\text{erfi}\left(\frac{u}{\sqrt{2}\sigma }\right)-\frac{4u^2}{P\left(u\right){\left(\Delta x\right)}^2}\right).$$(33)

Importantly, we have that $\underset{u\to \infty }{\lim }{I}_F\left({\mathcal{J}}_u\right)$ exists, which implies that $\underset{u\to \infty }{\lim }\frac{\partial I_F\left(u\right)}{\partial u}=0$.

The proof is presented in Appendix C. Then, the optimization formulation in Eq. (31) can be rewritten as $$\begin{array}{l}{u}_{I_F}^{*}=\arg \underset{u\in {\text{R+}}^{}}{\min }\left(P\left(u\right)-\frac{2u^2}{\pi {\sigma }^2\text{\hspace{0.17em}erfi}\left(\frac{u}{\sqrt{2}\sigma }\right)}\right)\hfill \\ \text{such\hspace{0.17em}that}\frac{\left(P\left(u_{I_F}^{*}\right)-\frac{2{\left(u_{I_F}^{*}\right)}^2}{\pi {\sigma }^2\text{\hspace{0.17em}erfi}\left(\frac{u_{I_F}^{*}}{\sqrt{2}\sigma }\right)}\right)}{\left(P\left(u_{\max }\right)-\frac{2u_{\max }^2}{\pi {\sigma }^2\text{\hspace{0.17em}erfi}\left(\frac{u_{\max }^2}{\sqrt{2}\sigma }\right)}\right)}\ge 1-\delta .\hfill \end{array}$$(34)

Remark 1. Theorem 2 tells us that both the joint Fisher information in flux and background have a monotonic behavior as a function of the aperture radius. This is important because it reinforces the results presented in Vicuña et al. (2024); specifically, and as expected, to achieve higher precision in estimation, it is necessary to use all available pixels. Any other method based on processing a subset of the available observations will inevitably lead to a loss of information. This is because, according to the data processing inequality (Cover 1999), the information about the parameter decreases with each processing step.

Remark 2. Equation (30) tells us that we can be arbitrarily close to the value that maximizes Fisher information by choosing an appropriate delta. However, the contribution per pixel gradually decreases, leading to a saturation zone. This feature allows for an aperture design by identifying the point where such saturation occurs. The proposed criterion provides flexibility as it allows for the adjustment of tolerance values based on the cost (or convenience) of incorporating more pixels. By doing so, we can determine the gain (or loss) in precision when additional pixels are included (or excluded) in the aperture depending on crowding. In practice, this scheme allows for the free selection of a tolerance level by using a larger number of pixels, which can be costly or inadequate (e.g., in a crowded field). Therefore, the flexibility provided by this method enables the astronomer to conduct a cost-benefit analysis to achieve the highest attainable precision given the observational constraints while using the fewest pixels necessary within the acceptable error budget of the respective science case.

4 Numerical results

In this section, we present numerical results, focusing on two key experiments. The first experiment examines the variation in Fisher information as a function of the aperture radius, and provides insights into how the choice of this radius impacts the information available for the joint estimation. The second experiment analyzes the performance of classical astrometric estimators as a function of the aperture radius, comparing their accuracy and precision under different aperture configurations in comparison with the expected CRLB for that aperture choice.

4.1 Discrepancy analysis as a function of the aperture radius

Here we present a numerical analysis that lends support to the methodology proposed in Eq. (30) and the theoretical expressions derived from Theorem 2. We use Eq. (30) to compute, for a selection of the PSF fitting radius, the Fisher information comparing its predictions to the aperture that maximizes the S/N in Eq. (25). To this end, we developed various experiments with different flux and background levels. The experimental setup (source, device, and sky conditions) considered here and in the forthcoming sections, which we denote as the baseline case, is summarized in Table 1. The parameter dither denotes an offset between the center of the $n_{x_c}^{th}$ pixel and the actual position of the point source x_c (Mendez et al. 2013, 2014), measured in pixels.

For example, for a flux of 10,002 [photo-e⁻], we obtain an S/N of 78⁹. The first experiment involves analyzing the evolution of the discrepancy in the joint CRLB for the photometric case as the aperture radius increases. Specifically, by discrepancy, we mean that for a given PSF fitting (aperture) radius u, we compute $\frac{I_F\left(\mathcal{N}\right)-I_F\left({\mathcal{J}}_u\right)}{I_F\left(\mathcal{N}\right)}$ for the source flux, and $\frac{I_B\left(\mathcal{N}\right)-I_B\left({\mathcal{J}}_u\right)}{I_B\left(\mathcal{N}\right)}$ for the background. We assess how far the (variance) bound deviates from the value obtained when considering all pixels. Additionally, the S/N is plotted as a function of the aperture radius, indicating the point where its value is maximized. Figure 1 shows the results of this experiment for different configurations of F and B.

In terms of the source flux, as expected, we observe a mono-tonic decrease in the CRLB variance (i.e., increasing precision); as the aperture radius increases, more information from the source is involved, leading to a reduction in the variance associated with the flux estimation. The discrepancy also decreases when the difference between the flux F and the background B increases, which is expected. A higher flux value relative to the background means a cleaner signal, thereby improving photometric precision. On the other hand, a larger FWHM implies a flux signal spread wider across the instrument, which results in a slower rate at which the best precision is reached when all pixels are available. Moreover, we note that the aperture radius that maximizes the S/N corresponds to a discrepancy percentage (in the terms defined previously) between 40% and 85%.

Table 2 presents the aperture for several popular choices of the PSF fitting radius for a FWHM of 0.5, 1.0 (baseline case, see Table 1), and 1.5 arcsec, and its corresponding discrepancy values. These findings suggest that there is room for improvement in the joint estimation of background and flux if a more appropriate radius is selected, for example one that limits the discrepancy to a pre-specified discrepancy (e.g., 10%).

In the case of the background, we observe that there is no sharp decline when the aperture radius increases. Obviously, the flux is dominant around the center of the pixel, whereas in pixels outside the aperture radius, the background becomes linearly dominant. Pixels outside the aperture radius are, therefore, more informative for the background than for the flux. This observation motivates the selection of pixels for determining the aperture radius: the background pixels should complement those selected for the flux. By focusing on the flux-dominant pixels in the center and the background-dominant pixels outside the aperture, we can optimize the aperture radius to enhance the accuracy of photometric measurements. This complementary selection ensures a more precise distinction between the flux and background, thereby improving the reliability of the estimates. This is explored in the next subsection.

4.1.1 Discrepancy analysis as a function of the halo

To further the analysis presented above, we consider the complement set of J_u, which includes the pixels located outside the aperture radius. This set corresponds to a halo, as it consists of pixels outside the central region and symmetrically distributed around it. Figure 2 presents an inverse analysis, showing the discrepancy as the halo encompasses more pixels until it includes all of them. Interestingly, achieving a more pronounced decay in the discrepancy requires a higher flux-to-background ratio, F/B. This relationship highlights the critical role of the F/B ratio in optimizing the precision and effectiveness of photometric analyses, particularly when considering contributions from regions beyond the immediate aperture. The Fisher information typically increases for the background estimation with a higher S/N because a stronger signal enhances the ability to disentangle the source flux from the background. Intuitively, one might expect the opposite—that a lower S/N would emphasize the background—because noise dominates in such conditions. However, the relationship is nuanced. When the S/N is high, the observed flux is more distinct from the background, making it easier to accurately model and estimate the contribution from background. This clarity increases the Fisher information for the background because the data provides more reliable clues about the background’s statistical behavior. In contrast, when the S/N is low, the noise dominates, and the background becomes harder to distinguish from the signal. This uncertainty reduces the Fisher information, as the data carries less discernible information about the background parameters.

Fig. 1 Discrepancy analysis of the joint Fisher information as a function of the PSF aperture fitting radius in photometry (blue line, left ordinate on each plot, in %). The discrepancy in the source flux (left column) is defined as ${\Delta }_{max}{I}_F\frac{I_F\left(\mathcal{N}\right)-I_F\left({\mathcal{J}}_u\right)}{I_F\left(\mathcal{N}\right)}$. The discrepancy in the background (right column) is defined as ${\Delta }_{max}{I}_B\frac{I_B\left(\mathcal{N}\right)-I_B\left({\mathcal{J}}_u\right)}{I_B\left(\mathcal{N}\right)}$. The black line (right ordinate on each plot) indicates the corresponding S/N value. The vertical line shows the aperture at which the S/N is maximized. Results are reported for different representative values of F, B, and a FWHM of 1.0 arcsec.

Fig. 2 Discrepancy analysis of the joint Fisher information (in %) as a function of the halo coverage in photometry. The left panel is for the source flux; the right panel is for the background. Results are reported for different representative values of F/B.

4.1.2 Impact of contaminating stars on the aperture optimization

Recent advances in stellar photometry, particularly from Gaia, have made it possible to precisely know the positions of neighboring stars in the field and to estimate their fluxes with high accuracy. This enables the explicit incorporation of contaminating stars into the photometric model for aperture optimization. In the standard framework, the expected counts in pixel i is modeled as in Eq. (3). This assumes that the background is spatially uniform. When contaminating stars are present, the background becomes spatially variable and can be modeled as $${\tilde{B}}_i=B+{\displaystyle \sum _{k=1}^K{F}_k{𝑔}_{ik}\left(x_{ck}\right),}$$(35)

where B is the constant background level, F_k is the flux of contaminating star k, and ɡ_ik(x_ck) is the PSF of star k integrated over pixel i, shifted according to its position relative to the target. The total expected count is therefore $${\lambda }_i\left(F,{\tilde{B}}_i\right)=F{𝑔}_i\left(x_{ck}\right)+{\tilde{B}}_i.$$(36)

In this scenario, B_i is assumed to be known from Gaia-based positions and pre-measured or modeled fluxes, and only F is to be estimated. The Fisher information for F simplifies to $$I_{1F}\left({𝒥}_u\right)={\displaystyle \sum _{i\in {𝒥}_u}\frac{{𝑔}_i^2\left(x_c\right)}{F{𝑔}_i\left(x_c\right)+{\tilde{B}}_i}.}$$(37)

Unlike the constant-background case, the background B_i is a fixed known function of pixel position.

For aperture optimization, the discrepancy at aperture radius u is defined as $$\text{disc}\left(u\right)=1-\frac{I_{1F}\left({𝒥}_u\right)}{I_{1F}\left(\mathcal{N}\right)},$$(38)

where I_1F(𝒥_u) is computed using only pixels within radius u from the target’s center, and I_1F(𝒩) is the Fisher information using the full set of pixels in the analysis window. This discrepancy quantifies the fractional loss of Fisher information when restricting the aperture to a given radius.

To quantify the impact of contamination from nearby stars, the background B_i was modeled including one or more contaminating sources at varying separations from the target. Each contaminant was assumed to share the same Gaussian PSF as the target, with a fixed FWHM, and its flux contribution was scaled according to its total flux F_k. The separation between the target and the contaminant was systematically varied from zero to several arcseconds. For each separation, the discrepancy was computed as a function of aperture radius.

The results are shown in Figure 3, where the horizontal axis represents both the aperture radius and the separation between the contaminating star (SBS) and the target, while the vertical axis displays the color scale of the normalized discrepancy. The left panel of Figure 3 corresponds to the case of a faint contaminant (F = 2000 [photo-e⁻], F_K = 1000 [photo-e⁻], K = 1), while the right panel illustrates the effect of a brighter one (F = 2000 [photo-e⁻], F_K = 2500 [photo-e⁻], K = 1). In both cases, a clear trade-off emerges: at small separations, when the contaminant lies within the core of the PSF, the discrepancy remains high across all aperture sizes. This indicates a significant loss of Fisher information due to the strong overlap between the target and contaminant. As the separation increases, the contamination effect decreases, and the discrepancy curves approach those of the uncontaminated case. For large separations, the contamination is negligible and the optimal aperture radius converges to the same value as in the constant-background scenario.

It is worth noting that, from the perspective of the model, the effect of a single bright contaminant located close to the target can be interpreted as equivalent to the combined effect of multiple fainter stars at similar distances since what ultimately matters is the cumulative flux contribution of contaminants to the pixel i through the term $B_i=B+{\displaystyle {\sum }_{k=1}^K{F}_k{g}_{ik}\left(x_{ck}\right)}$. Consequently, the analysis of the aperture optimization problem in the presence of multiple nearby stars reduces to understanding how these additional flux terms modify the spatial profile of B_i and, therefore, the Fisher information curve as a function of aperture radius.

For the purposes of this work, we restricted our attention to the case of a single well-isolated source, which allowed us to focus on the fundamental behavior of the methodology without the added complexity of crowding effects. A detailed treatment of how nearby sources impact photometric precision and the choice of the optimal aperture will be the subject of future work.

4.2 Discrepancy analysis for practical high dimensional implicit estimators

We now analyze the behavior of some classical estimators used in photometry for flux estimation as a function of the aperture radius. Given the equivalence between aperture radius and the number of pixels used, it is important to analyze how much precision is lost when fewer observations are employed. In these numerical experiments, we examine how much information is available depending on the aperture selection criterion, comparing the use of the S/N-based approach versus the Fisher information-based approach. First, consider an estimator, or regression rule τ_J : ℕⁿ → Θ, τ_J(·) = [τ_J,1(·), …,τ_J,p(·)], which is defined by the solution of $${\tau }_J\left(\mathbf{\text{I}}\right)\equiv \underset{\mathit{\alpha }\in \mathbf{\Theta }}{\mathrm{argmin}}J\left(\mathit{\alpha },\mathbf{\text{I}}\right),$$(39)

where J(α, I) is a general cost function. From this point forward, we consider two classical estimators commonly used in astrometry, along with their associated cost functions: the maximum likelihood estimator and the stochastic weighted least squares estimator. Both estimators have been formally analyzed in detail, and provide theoretical guarantees of optimality in a broad sense across various astrometric and photometric scenarios (see Lobos et al. 2015; Espinosa et al. 2018; Vicuña et al. 2024).

The stochastic weighted least squares estimator is given by the following cost function: $$J_{SWLS}\left(\mathit{\alpha },\mathbf{\text{I}}\right)={\displaystyle \sum _{i=1}^n\frac{1}{I_i}{\left(I_i-{\alpha }_1\cdot {𝑔}_i\left(x_c\right)-{\alpha }_2\right)}^2.}$$(40)

This is a particular case of weighted least squares where the selection of the weights can be seen as a noisy version of 𝔼{I_i} = λ_i(F, B). On the other hand, the ML estimator has been widely used for parameter estimation due to its asymptotic optimality (Kay 1993). For this reason, it has been proposed and used by the astronomical community in photometry and astrometry (Chun-Lin 1993; Gai et al. 2013; Lindegren 2008; Gai et al. 2017; Espinosa et al. 2018), and in many other research areas facing similar inverse image problems, such as fluorescence microscopy (Abraham et al. 2009; Small & Parthasarathy 2014). The likelihood function of I given the source’s relative position x_c and its photometry α ∈ Θ (see Eq. (7)), can be expressed as $$L\left(\mathbf{\text{I;\hspace{0.17em}}}\alpha ,\mathcal{N}\right)={\displaystyle \prod _{i=1}^n\frac{e^{-\left({\alpha }_1\cdot {𝑔}_1\left(x_c\right)+{\alpha }_2\right)}\cdot {\left({\alpha }_1\cdot {𝑔}_1\left(x_c\right)+{\alpha }_2\right)}^{I_i}}{I_i!}.}$$(41)

Then, the ML estimator is given by $${\tau }_{ML}\left(\mathbf{\text{I}}\right)=\underset{\alpha \in \mathbf{\Theta }}{\mathrm{argmin}}\underset{=-\ln L\left(\mathbf{\text{I;\hspace{0.17em}}}\alpha ,\mathcal{N}\right)}{\underbrace{J_{ML}\left(\mathit{\alpha },\mathbf{\text{I}}\right)}}.$$(42)

Using these estimators, the second experiment involves quantifying the change in performance of these estimators as a function of the aperture radius. The aim is to see how the use of fewer pixels or a smaller aperture radius affects the precision of these estimators. Figure D.1 (see Appendix D) shows the numerical results where we analyze the empirical variance of each estimator and compare it with their variance when all pixels are available. The black vertical line represents the aperture radius associated with the value that maximizes the S/N, while the blue vertical line corresponds to the radius selection based on the Fisher information (i.e., the CRLB), allowing a discrepancy of 20% according to Eq. (30). The green vertical line corresponds to the value of the aperture radius following the approximated methodology in Eq. (34). We compared the SWLS and ML estimators as a function of the aperture radius. More precisely, we analyzed the following merit function (expressed as % of discrepancy with respect to the case of the full pixel¹⁰); for visualization, the variances of the source’s flux are expressed in terms of magnitudes (Mendez et al. 2014): $${\sigma }_m^2\left(\mathbf{\text{I}}\right)=\frac{2.5}{2}\cdot \log \left(F+\sqrt{Var\left({\tau }_{J,F}\left(\mathbf{\text{I}}\right)\right)}\right)-\log \left(F-\sqrt{Var\left({\tau }_{J,F}\left(\mathbf{\text{I}}\right)\right)}\right),$$(43) $${\Delta }_{\max }{\widehat{\sigma }}_m\left(u\right)=\frac{{\sigma }_m^2\left(\mathbf{\text{I}}\right)-{\sigma }_m^2\left(\mathbf{\text{I}}\right)}{{\sigma }_m^2\left(\mathbf{\text{I}}\right)}\cdot 100,$$(44) $${\Delta }_{\max }{\widehat{\sigma }}_B\left(u\right)=\frac{\sqrt{Var\left({\tau }_{J,B}\left(I_u\right)\right)}-\sqrt{Var\left({\tau }_{J,B}\left(I\right)\right)}}{\sqrt{Var\left({\tau }_{J,B}\left(I\right)\right)}}\cdot 100$$(45)

Here τ_J,F(·) and τ_J,B(·) denote the estimated values of the source flux and background, respectively; J represents the cost function of each estimator, namely stochastic weighted least squares or maximum likelihood; and I_u denotes the observations located within the aperture radius u. As anticipated, we can see that the discrepancy significantly decreases as the aperture radius increases. This result is consistent with the findings from the Fisher information case. Thus, in empirical terms, the performance comparison of the estimators aligns with the fundamental limits. We also note that, for these scenarios, the discrepancy can reach up to 180% when using the aperture that maximizes the S/N, indicating that the estimation process could be significantly improved if Fisher information is considered. In such a case, the empirical variance would be much closer to the optimal value when more pixels are used.

Fig. 3 Normalized discrepancy (vertical axis, Eq. (38)) as a function of aperture radius and separation between the contaminating star (SBS) and the target (horizontal axis), for two contamination scenarios. The left panel corresponds to a faint contaminant (F = 2000 [photo-e−], FK = 1000 [photo-e−], K = 1), while the right panel shows the effect of a brighter one (F = 2000 [photo-e−], FK = 2500 [photo-e−], K = 1). Colors indicate the fractional loss of Fisher information relative to the full-pixel case.

5 Conclusions

This research underscores the relevance of using the Fisher information matrix to significantly enhance photometric analysis through an optimal (in the CRLB sense) PSF fitting radius (aperture selection) in crowded fields. By calculating the joint Fisher information matrix for photometric measurements as a function of the aperture radius, we gain valuable insights into how different aperture radii influence the precision of source flux estimates. The joint Fisher information effectively quantifies the amount of information that each pixel or measurement point within an aperture contributes to the final photometric estimate, enabling a more robust data-driven approach to defining the aperture radius.

Our findings show that while increasing the aperture radius improves the joint Fisher information, and thus enhances the accuracy of flux measurements, there is a practical trade-off. Larger apertures incorporate more pixels, increase computational costs, and eventually run into crowding problems. Consequently, the framework developed here allows for an informed balance between aperture size and the precision of flux estimation, enabling the selection of an optimal aperture that minimizes the loss of information from fewer pixels without unnecessarily inflating computational requirements or incurring in crowding issues. The results suggest that selecting a more appropriate aperture radius (e.g, one that limits the discrepancy to pre-specified percentage) can lead to improvements in the flux estimation process. In some cases, using S/N-based criteria results in discrepancies as large as 80% compared to choosing a larger aperture radius. This indicates that a more precise radius selection can significantly reduce estimation errors and improve the reliability of the photometric measurements.

The results presented in this paper are particularly relevant, considering that the joint estimation photometric methodology presented in Vicuña et al. (2024) assures variances that are very close to the CRLB and that can be attained by practical estimators (as demonstrated by our realistic numerical simulations). Hence, the Cramér-Rao limit is a meaningful benchmark to provide a prescription for the optimal aperture in photometry.

As seen in the numerical results presented in this paper, when the aperture radius is too small, valuable information from the outer pixels of the source’s light distribution is neglected, potentially decreasing the estimator’s precision of flux determination. This methodological advancement shifts photometric practices from traditional, often heuristic, aperture choices to a statistically optimized framework that can adapt based on the source characteristics, background conditions, crowding conditions, and instrumental setup. Finally, this approach ultimately enables more accurate and reliable photometric measurements, providing a powerful tool for precision astronomy in studies ranging from exoplanet detection to deep-sky observations.

Acknowledgements

SAET, JFS, and MO acknowledge support from the Advanced Center for Electrical and Electronic Engineering (AC3E), funded by the Chilean ANID Basal Project (AFB240002). Additionally, JFS acknowledges support from Fondecyt Project 1250098. RAM acknowledges support from FONDECYT/ANID grant # 124 0049 and from ANID, Fondo GEMINI, Astrónomo de Soporte GEMINI-ANID grant # 3223 AS0002.

Appendix A Proof of Proposition 1

Equation (22) is a natural consequence of (20) with a Gaussian PSF and considering a set of indices of the form ${\mathcal{J}}_u=\left\{n_{x_c}-n_u,\dots ,n_{x_c},\dots ,n_{x_c}+n_u\right\}$: $$\begin{array}{l}{\displaystyle \sum _{i\in {𝒥}_u}{𝑔}_i^p\left(x_c\right)}\hfill \\ \approx {\displaystyle \sum _{i\in {𝒥}_u}{\varphi }^p\left(x_i-x_c,\sigma \right)\cdot \Delta x^p}\hfill \\ ={\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^p\cdot \Delta x^{p-1}\cdot {\displaystyle \sum _{i\in {𝒥}_u}\exp \left(-p\cdot \frac{{\left(x_i-x_c\right)}^2}{2{\sigma }^2}\right)\cdot \Delta x}\hfill \\ \approx {\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^p\cdot \Delta x^{p-1}\cdot \underset{\Delta x\to 0}{\lim }{\displaystyle \sum _{i\in {𝒥}_u}\exp \left(-\frac{{\left(x_i-x_c\right)}^2}{\frac{2}{p}{\sigma }^2}\right)\cdot \Delta x}\hfill \\ \approx {\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^p\cdot \Delta x^{p-1}\cdot {\displaystyle {\int }_{x_c-u}^{x_c+u}\exp \left(-\frac{{\left(x-x_c\right)}^2}{\frac{2}{p}{\sigma }^2}\right)}dx\hfill \\ ={\left(\frac{\Delta x}{\sqrt{2\pi }\sigma }\right)}^{p-1}\cdot \frac{1}{\sqrt{2\pi }\sigma }\cdot 2{\displaystyle {\int }_0^u\exp \left(-\frac{x^2}{\frac{2}{p}{\sigma }^2}\right)}dx\hfill \\ ={\left(\frac{\Delta x}{\sqrt{2\pi }\sigma }\right)}^{p-1}\cdot \frac{1}{\sqrt{p}}\cdot \frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^{\frac{u}{\sqrt{2}\sigma }\cdot \sqrt{p}}\exp \left(-x^2\right)}dx\hfill \\ =\frac{1}{\sqrt{p}}\cdot {\left(\frac{\Delta x}{\sqrt{2\pi }\sigma }\right)}^{p-1}\cdot \text{erf}\left(\frac{u}{\sqrt{2}\sigma }\sqrt{p}\right)\hfill \\ =\frac{1}{\sqrt{p}}\cdot {\left(\frac{1}{\sqrt{2\pi }}\right)}^{p-1}\cdot {\left(\frac{\Delta x}{\sigma }\right)}^{p-1}\cdot P\left(u\sqrt{p}\right).\hfill \end{array}$$(A.1)

Following similar steps, the approximation of sums of negative powers of ɡ_i(x_c) can be found: $$\begin{array}{l}{\displaystyle \sum _{i\in {𝒥}_u}\frac{1}{{𝑔}_i^p\left(x_c\right)}}\hfill \\ ={\displaystyle \sum _{i\in {𝒥}_u}{𝑔}_i^{-p}\left(x_c\right)}\hfill \\ \approx {\displaystyle \sum _{i\in {𝒥}_u}{\varphi }^p\left(x_i-x_c,\sigma \right)\cdot \Delta x^{-p}}\hfill \\ ={\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^{-p}\cdot \Delta x^{-p-1}{\displaystyle \sum _{i\in {𝒥}_u}\exp \left(p\cdot \frac{{\left(x_i-x_c\right)}^2}{2{\sigma }^2}\right)\cdot \Delta x}\hfill \\ \approx {\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^{-p}\cdot \Delta x^{-p-1}\underset{\Delta x\to 0}{\lim }{\displaystyle \sum _{i\in {𝒥}_u}\exp \left(\frac{{\left(x_i-x_c\right)}^2}{\frac{2}{p}{\sigma }^2}\right)\cdot \Delta x}\hfill \\ \approx {\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^{-p}\cdot \Delta x^{-p-1}{\displaystyle {\int }_{x_c-u}^{x_c+u}\exp \left(\frac{{\left(x-x_c\right)}^2}{\frac{2}{p}{\sigma }^2}\right)}dx\hfill \\ ={\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^{-p}\cdot \Delta x^{-p-1}\cdot 2{\displaystyle {\int }_0^u\exp \left(\frac{x^2}{\frac{2}{p}{\sigma }^2}\right)}dx\hfill \\ ={\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^{-p}\Delta x^{-p-1}\cdot \frac{\sqrt{2}\sigma }{\sqrt{p}}\cdot 2{\displaystyle {\int }_0^{\frac{u}{\sqrt{2}\sigma }\cdot \sqrt{p}}\exp \left(x^2\right)}dx\hfill \\ =\frac{1}{\sqrt{p}}\cdot {\left(\frac{\Delta x}{\sqrt{2\pi }\sigma }\right)}^{-p-1}\cdot \frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^{\frac{u}{\sqrt{2}\sigma }\cdot \sqrt{p}}\exp \left(x^2\right)}dx\hfill \\ =\frac{1}{\sqrt{p}}\cdot {\left(\frac{1}{\sqrt{2\pi }\sigma }\right)}^{-p-1}\cdot \text{erfi}\left(\frac{u}{\sqrt{2}\sigma }\sqrt{p}\right).\hfill \end{array}$$(A.2)

Here the last steps come from the definition of the erfi function: $$\begin{array}{l}\text{erfi\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}=-i\text{\hspace{0.17em}erf\hspace{0.17em}}\left(ix\right)\\ =-i\cdot \frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^{ix}{e}^{-{\upsilon }^2}}d\upsilon \\ =-i\cdot i\cdot \frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^x{e}^{{\upsilon }^2}}d\upsilon \\ =\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^x{e}^{{\upsilon }^2}}d\upsilon .\end{array}$$(A.3)

Appendix B Proof of Proposition 2

The proof follows from the application of the first derivative criterion, we see that $$\begin{array}{c}\frac{\partial S/N\left(u\right)}{\partial u}=\frac{2\sqrt{2}\sigma P^{\prime }\left(\sqrt{2}\sigma u\right)F\left[P\left(\sqrt{2}\sigma u\right)F+4{\sigma }^{2}u\left(B-D\right)\right]}{2{\left[P\left(\sqrt{2}\sigma u\right)F+4{\sigma }^{2}u\left(B-D\right)\right]}^{3/2}}\\ -\frac{P\left(\sqrt{2}\sigma u\right)F\left[\sqrt{2}\sigma P^{\prime }\left(\sqrt{2}\sigma u\right)F+4{\sigma }^{2}u\left(B-D\right)\right]}{2{\left[P\left(\sqrt{2}\sigma u\right)F+4{\sigma }^{2}u\left(B-D\right)\right]}^{3/2}}.\end{array}$$(B.1)

By imposing $\frac{\partial S/N\left(u\right)}{\partial u}{|}_{u_{S/N}^{*}}=0$ and rearranging terms we obtain $$\begin{array}{l}F{P}^{\prime }\left(\sqrt{2}\sigma u_{S/N}^{*}\right)P\left(\sqrt{2}\sigma u_{S/N}^{*}\right)=\hfill \\ -\left[2\sqrt{2}\sigma u_{S/N}^{*}{P}^{\prime }\left(\sqrt{2}\sigma u_{S/N}^{*}\right)-P\left(\sqrt{2}\sigma u_{S/N}^{*}\right)\right]2\sqrt{2}\sigma \left(B-D\right).\hfill \end{array}$$(B.2)

Appendix C Proof of Theorem 2

The proof of Theorem 2 involves three parts. First, the elements of the Fisher matrix are approximated according to the assumption of dominant flux. The second part consists of approximating the joint Fisher information. The third part verifies that the monotonicity conditions hold for both approximations.

Appendix C.1 Approximations for the elements of (27)

Due to the discrete nature of the Fisher information matrix, the objective of this approach comes from the need to obtain a differentiable and continuous version of the terms in this matrix as a function of the aperture radius. Since the flux is more dominant near the central pixel, we assume a dominant flux context for these approximations, which is equivalent to a zero-order Taylor expansion. The zero-order approximations developed in this section follow a path similar to those presented in Mendez et al. (2013), which though naïve, serve as a basic interpretation of how information availability depends on scenario configuration. We have the following: $$\begin{array}{c}{I}_{1,1}\left({𝒥}_u\right)={\displaystyle \sum _{i\in {𝒥}_u}\frac{{𝑔}_i^2\left(x_c\right)}{F\cdot {𝑔}_i\left(x_c\right)+B}}\\ =\frac{1}{F}{\displaystyle \sum _{i\in {𝒥}_u}\frac{{𝑔}_i^2\left(x_c\right)}{{𝑔}_i\left(x_c\right)+\frac{B}{F}}}\\ \approx \frac{1}{F}{\displaystyle \sum _{i\in {𝒥}_u}{𝑔}_i\left(x_c\right)\approx \frac{1}{F}\cdot P\left(u\right)}\end{array}$$(C.1) $$\begin{array}{c}{I}_{1,2}\left({\mathcal{J}}_u\right)=I_{2,1}\left({\mathcal{J}}_u\right)={\displaystyle \sum }_{i\in {\mathcal{J}}_u}\frac{{𝑔}_i\left(x_c\right)}{F\cdot {𝑔}_i\left(x_c\right)+B}\\ =\frac{1}{F}{\displaystyle \sum }_{i\in {\mathcal{J}}_u}\frac{{𝑔}_i\left(x_c\right)}{{𝑔}_i\left(x_c\right)+\frac{B}{F}}\\ \approx \frac{1}{F}{\displaystyle \sum }_{i\in {\mathcal{J}}_u}1\\ =\frac{\left|{\mathcal{J}}_u\right|}{F}\\ =\frac{1}{F}\cdot \frac{2u}{\Delta x}\end{array}$$(C.2) $$\begin{array}{c}{I}_{2,2}\left({𝒥}_u\right)={\displaystyle \sum _{i\in {𝒥}_u}\frac{1}{F\cdot {𝑔}_i\left(x_c\right)+B}}\\ =\frac{1}{F}{\displaystyle \sum _{i\in {𝒥}_u}\frac{1}{{𝑔}_i\left(x_c\right)+\frac{B}{F}}}\\ \approx \frac{1}{F}{\displaystyle \sum _{i\in {𝒥}_u}\frac{1}{{𝑔}_i\left(x_c\right)}}\\ \approx \frac{1}{F}2\pi {\left(\frac{\sigma }{\Delta x}\right)}^2\text{erfi}\left(\frac{u}{\sqrt{2}\sigma }\right).\end{array}$$(C.3)

Appendix C.2 Approximation of the joint Fisher information

Using the approximations for the elements of (27): (C.1), (C.2) and (C.1) we have that $$\begin{array}{c}{I}_F\left({𝒥}_u\right)=I_{1,1}\left({𝒥}_u\right)-\frac{I_{1,2}^2\left({𝒥}_u\right)}{I_{2,2}\left({𝒥}_u\right)}\\ \approx \frac{1}{F}\cdot P\left(u\right)-\frac{{\left(\frac{1}{F}\cdot \frac{2u}{\Delta x}\right)}^2}{\frac{1}{F}2\pi {\left(\frac{\sigma }{\Delta x}\right)}^2\text{erfi}\left(\frac{u}{\sqrt{2}\sigma }\right)}\\ =\frac{1}{F}\left(P\left(u\right)-\frac{2u^2}{\pi {\sigma }^2\text{\hspace{0.17em}erfi}\left(\frac{u}{\sqrt{2}\sigma }\right)}\right).\end{array}$$ $$\begin{array}{c}{I}_B\left({𝒥}_u\right)=I_{2,2}\left({𝒥}_u\right)-\frac{I_{1,2}^2\left({𝒥}_u\right)}{I_{1,1}\left({𝒥}_u\right)}\\ =\frac{1}{F}2\pi {\left(\frac{\sigma }{\Delta x}\right)}^2\text{erfi}\left(\frac{u}{\sqrt{2}\sigma }\right)-\frac{{\left(\frac{1}{F}\cdot \frac{2u}{\Delta x}\right)}^2}{\frac{1}{F}\cdot \left(u\right)}\\ =\frac{1}{F}\left(2\pi {\left(\frac{\sigma }{\Delta x}\right)}^2\text{erfi}\left(\frac{u}{\sqrt{2}\sigma }\right)-\frac{4u^4}{P\left(u\right){\left(\Delta x\right)}^2}\right).\end{array}$$

Appendix C.3 Monotonicity of the joint Fisher information

For the last part of the proof we see that $$\underset{u\to \infty }{\lim }\frac{u^2}{\text{erfi}\left(\frac{u}{\sqrt{2}\sigma }\right)}=0$$(C.4) $$\underset{u\to \infty }{\lim }P\left(u\right)=1,$$(C.5)

which means that: $$\underset{u\to \infty }{\lim }{I}_F\left({𝒥}_u\right)=\frac{1}{F}.$$(C.6)

Appendix D Additional figure

Fig. D.1 Discrepancy analysis of the estimators as a function of the PSF aperture fitting radius in photometry (orange line, left ordinate on each plot, in %). The discrepancy in the source flux (left column) is defined in Eq. (44). The discrepancy in the background (right column) is defined in Eq. (45). The broad black line (right ordinate on each plot) indicates the corresponding S/N value. The vertical black line shows the aperture at which the S/N is maximized. The vertical blue line corresponds to the value of the aperture radius following the Fisher information methodology proposed in Eq. (30), while the green vertical line corresponds to the value of the aperture radius following the approximated methodology in Eq. (34) (in some cases, the blue and green vertical lines overlap). The achieved discrepancy is observed to be close—but not necessarily equal—to the threshold value of 20% due to the discrete nature of the summation. Specifically, ∆maxIF denotes the last value of the normalized discrepancy that remains below the selected threshold δ as the aperture radius increases, which may not exactly match δ. The choice of δ = 20% is arbitrary and is adopted here purely for illustrative purposes. Results are reported for different representative values of F, B and a FWHM of 1.0 arcsec.

References

All Tables

All Figures

Fig. 1 Discrepancy analysis of the joint Fisher information as a function of the PSF aperture fitting radius in photometry (blue line, left ordinate on each plot, in %). The discrepancy in the source flux (left column) is defined as ${\Delta }_{max}{I}_F\frac{I_F\left(\mathcal{N}\right)-I_F\left({\mathcal{J}}_u\right)}{I_F\left(\mathcal{N}\right)}$. The discrepancy in the background (right column) is defined as ${\Delta }_{max}{I}_B\frac{I_B\left(\mathcal{N}\right)-I_B\left({\mathcal{J}}_u\right)}{I_B\left(\mathcal{N}\right)}$. The black line (right ordinate on each plot) indicates the corresponding S/N value. The vertical line shows the aperture at which the S/N is maximized. Results are reported for different representative values of F, B, and a FWHM of 1.0 arcsec.

In the text

	Fig. 2 Discrepancy analysis of the joint Fisher information (in %) as a function of the halo coverage in photometry. The left panel is for the source flux; the right panel is for the background. Results are reported for different representative values of F/B.
In the text