Methods

PyDIET combines spectral and image models to predict the signal and noise in an observation (Fig. 3). This chapter explains the main quantities and methods behind the calculation.

Fig. 3 Schematic view of the PyDIET simulation workflow.

Signal-to-noise ratio

The signal-to-noise ratio (SNR) is the ratio of the source signal \(F\) to its expected uncertainty \(\sqrt{\var{F}}\). Here, "flux" refers to the number of source photons collected by the detector during the observation:

(1)\[\mathrm{SNR} = \frac{F}{\sqrt{\var{F}}}.\]

The signal \(F=f.t\) follows from the expected photon rate \(f\) and exposure time \(t\). Its uncertainty depends on how the source flux is estimated. The SNR estimator in (1) is biased, but this is usually negligible for exposure-time planning at SNR \(\gtrsim 3\).

Model fitting

If photon noise follows Poisson statistics, readout noise is Gaussian, and noise is independent between pixels, fitting a source model with inverse-variance-weighted least-squares regression approximates the maximum likelihood estimate (MLE) of the flux \(F\). Other effects, such as background-estimation errors and contamination by neighboring objects, also contribute to flux uncertainty. These are difficult to quantify and are omitted from the calculator.

Forced photometry

First, suppose the source position is known, so the only fitted parameter is its amplitude, represented by the flux \(F\). Let \(\vec{\phi}\) be the source image model: the PSF for a point source, or the PSF convolved with a light-distribution model for a resolved source, normalized such that \(\sum_i \phi_i = 1\). The expected source contribution in pixel \(i\) is then \(F\phi_i\).

Before background subtraction, the pixel value may be written as

(2)\[Y_i = N_i + R_i,\]

where

(3)\[N_i \sim \operatorname{Poisson}(B_i + F\phi_i)\]

is the contribution from photons,

(4)\[R_i \sim \mathcal{N}(0,\sigma_r^2)\]

the readout noise, and

(5)\[B_i = \sigma_b^2 t\]

the expected background contribution. After subtracting the expected background, \(I_i = Y_i - B_i\), we have

(6)\[\esp{I_i} = F\phi_i,\]

and

(7)\[V_i(F) = \var{I_i} = \sigma_r^2 + \sigma_b^2 t + F\phi_i.\]

The exact likelihood of \(Y_i\) is the convolution of a Poisson distribution with a Gaussian readout noise distribution. Its Cramér-Rao bound has no simple closed form in general. For this ETC, we use an inverse-variance-weighted Gaussian approximation to the exact Poisson-plus-Gaussian likelihood. The inverse-variance-weighted estimating equation is obtained by minimizing

(8)\[\chi^2(F) = \sum_i \frac{\left(I_i - F\phi_i\right)^2} {V_i(F)},\]

where \(V_i(F)\) is the expected variance used to weight the residuals. The flux estimate satisfies

(9)\[U(F) = \sum_i \frac{\phi_i\left(I_i - F\phi_i\right)}{V_i(F)} = 0.\]

If the weights are kept fixed during one least-squares step, this gives the usual inverse-variance-weighted flux estimator [1]:

(10)\[\hat{F} = \frac{\sum_i \phi_i I_i / V_i}{\sum_i \phi_i^2 / V_i}.\]

In practice \(V_i\) may be updated iteratively using the current value of \(F\).

At the true value of \(F\),

(11)\[\esp{I_i - F\phi_i} = 0,\]

so \(\esp{U(F)} = 0\). The Fisher information associated with this optimal high-count approximation is

(12)\[\begin{split}\begin{aligned} {\cal I}(F) &= \var{U(F)} \\ &= \sum_i \frac{\phi_i^2}{V_i(F)^2} \var{I_i} \\ &= \sum_i \frac{\phi_i^2}{V_i(F)}. \end{aligned}\end{split}\]

Hence

(13)\[{\cal I}(F) \simeq \sum_i \frac{\phi_i^2}{\sigma_r^2 + \sigma_b^2 t + F\phi_i}.\]

For \(\sigma_r = 0\), this expression coincides with the exact Fisher information of the Poisson distribution. With non-zero readout noise, it is the usual Poisson+Gaussian approximation. The Cramér-Rao bound gives

(14)\[\var{\hat{F}} \gtrsim \frac{1}{{\cal I}(F)} = \left[ \sum_i \frac{\phi_i^2}{\sigma_r^2 + \sigma_b^2 t + F\phi_i} \right]^{-1}.\]

Substituting into (1) gives an upper bound on the SNR at exposure time \(t\):

(15)\[\begin{split}\begin{aligned} \mathrm{SNR} &\lesssim\, F\, \sqrt{ \sum_i \frac{\phi_i^2}{\sigma_r^2 + \sigma_b^2 t + F\phi_i}} \\ &\lesssim\,f.t\, \sqrt{ \sum_i \frac{\phi_i^2}{\sigma_r^2 + \left(\sigma_b^2 + f\phi_i\right)t}}. \end{aligned}\end{split}\]

\(\var{\hat{F}}\) approaches the Cramér-Rao bound at sufficiently large SNR or in the background-dominated Gaussian regime. At very low SNR, especially when flux is constrained to be non-negative or the source position is fitted simultaneously (see below), the uncertainty distribution may be asymmetric or biased. In this regime, treat (15) as an optimistic upper bound on the SNR rather than a precise prediction. The expression also omits imperfect background subtraction, correlated noise, model mismatch, source blending, and flat-field errors.

Fitting the source position

If the source coordinates are also left as free parameters, the relevant bound is no longer the reciprocal of the flux-only Fisher information. Instead, one must use the inverse of the Fisher matrix \(\boldsymbol{\cal I}\) for \(\vec{\theta} = (F,x_0,y_0)\): in the same approximation,

(16)\[{\cal I}_{\alpha\beta} \simeq \sum_i \frac{1}{V_i} \frac{\partial \mu_i}{\partial \theta_\alpha} \frac{\partial \mu_i}{\partial \theta_\beta},\]

where

(17)\[ \mu_i = B_i + F\phi_i(x_0,y_0),\]

and

(18)\[ \frac{\partial \mu_i}{\partial F} = \phi_i,\]

(19)\[ \frac{\partial \mu_i}{\partial x_0} = F\frac{\partial \phi_i}{\partial x_0},\]

(20)\[ \frac{\partial \mu_i}{\partial y_0} = F\frac{\partial \phi_i}{\partial y_0}.\]

The bound on the flux variance is then

(21)\[ \var{\hat{F}} \geq \left(\boldsymbol{\cal I}^{-1}\right)_{FF}.\]

Equivalently, writing

(22)\[\begin{split} \mathbf{c} = \begin{pmatrix} {\cal I}_{Fx} \\ {\cal I}_{Fy} \end{pmatrix},\end{split}\]

and

(23)\[\begin{split} \mathbf{Q} = \begin{pmatrix} {\cal I}_{xx} & {\cal I}_{xy} \\ {\cal I}_{xy} & {\cal I}_{yy} \end{pmatrix},\end{split}\]

one obtains

(24)\[ \var{\hat{F}} \geq \frac{1}{{\cal I}_{FF} - \mathbf{c}^{\mathsf T}\mathbf{Q}^{-1}\mathbf{c}}.\]

The correction term is positive, so fitting the position can only increase the flux uncertainty. For an isolated source with a centrally symmetric model (and a symmetric fitting domain), the flux-position cross-terms vanish by symmetry: \({\cal I}_{Fx} = {\cal I}_{Fy} = 0\). In that ideal case,

(25)\[ \left(\boldsymbol{\cal I}^{-1}\right)_{FF} = \frac{1}{{\cal I}_{FF}},\]

and the flux SNR bound reduces to (15). In real data, departures from symmetry, masking, undersampling, neighboring sources, or spatially varying noise can make the cross-terms non-zero, in which case fitting the position increases the bound on \(\var{\hat{F}}\).

Aperture photometry

Aperture photometry is another way of estimating the source flux. It is less statistically efficient than fitting an accurate source model, but requires less computation and no profile model. This can be useful for bright, unsaturated stars, whose model-fitting accuracy may be limited by the accuracy of the PSF model. For sources with unknown light profiles, aperture photometry is generally the safest option.

For an aperture \({\cal A}\), the SNR is

(26)\[\mathrm{SNR}({\cal A}) = \frac{f.t\sum_{i\in{\cal A}} \,\phi_i}{\sqrt{ \sum_{i\in{\cal A}} \sigma_r^2 + \left(\sigma_b^2 + f\phi_i\right)t }}.\]

Optimal aperture

The optimal-aperture option chooses the radius \(r=r_\mathrm{opt}\) of the circular aperture \({\cal A}(r)\) that maximizes the expected aperture signal-to-noise ratio \(\mathrm{SNR}({\cal A})\) in (26):

(27)\[ r_\mathrm{opt} = \underset{r}{\operatorname{arg\,max}}\,\mathrm{SNR}({\cal A}(r)).\]

The optimization is performed numerically over the allowed aperture radii, from zero to the radius of the circular simulation mask.

The optimal aperture depends on the exposure time and noise regime. PyDIET therefore recomputes the optimal aperture every time the signal-to-noise ratio is evaluated.

For stars, the optimal aperture depends on the seeing, sky background, source brightness, and exposure time. For galaxies, the same calculation is applied to the seeing-convolved Sérsic model (see Galaxy model). The optimal aperture depends on the Sérsic radius, Sérsic index, seeing, sky background, source brightness, and exposure time.

Aperture enclosing 96% of the flux

The large-aperture option is defined as the circular aperture enclosing 96% of the model source flux.

PyDIET computes the enclosed model-flux fraction

(28)\[f(r^2) = \sum_{r_i \leq r} \phi_i,\]

where \(r_i\) is the distance of pixel \(i\) from the source center on the internal raster grid.

The aperture radius is obtained by solving

(29)\[ f(r_{96}^2) = 0.96.\]

PyDIET solves (29) numerically with a Brent root finder. The search interval extends from the center of the image to the edge of the circular simulation mask.

The 96% aperture is close to a total-flux aperture. It captures most of the model flux, but it is not generally the aperture with the highest signal-to-noise ratio, because the outer parts of the aperture may contain more noise than useful source signal. In background-dominated conditions, the 96% flux aperture is usually larger than the optimal aperture. The outer wings of the source profile contain little signal but contribute sky noise, so including them lowers the signal-to-noise ratio. In source-dominated or very high signal-to-noise regimes, it may become smaller, because the penalty for including more of the source profile is reduced.

Exposure time

To find the exposure time needed for a target SNR, PyDIET solves (15) or (26) by iterative root finding, using Brent's method [2].

Source models

PyDIET provides models for point sources, galaxies, and uniform extended sources. Each describes how source photons are distributed over the detector pixels.

For compact and galaxy-like sources, the source model is represented by a discrete image \(\phi_i\), normalized as

(30)\[\sum_i \phi_i = 1,\]

where \(i\) runs over image pixels (or oversampled sub-pixels, depending on image sampling). If the total detected source rate is \(f\), in photons per second, the expected source contribution in pixel \(i\) during an exposure time \(t\) is therefore

(31)\[f.t.\phi_i.\]

For flat extended sources, the model is instead expressed as a surface brightness: the source rate is then interpreted per unit angular area rather than as the total rate of a finite object.

Point-source model

A point source is modeled as an unresolved object observed through the PSF. PyDIET uses a circular Moffat profile [3],

(32)\[M(r) = \left(1 + \frac{r^2}{\alpha^2}\right)^{-\beta},\]

where \(r\) is the angular distance from the source center, \(\alpha\) the Moffat scale radius, and \(\beta\) the parameter controlling the wings. The Moffat scale radius is set by the requested FWHM. Since \(M(\mathrm{FWHM}/2) = M(0)/2\), one obtains

(33)\[\alpha^2 = \frac{\mathrm{FWHM}^2}{4\left(2^{1/\beta} - 1\right)}.\]

The Moffat \(\beta\) parameter is set for each instrument in the data configuration file. The default value is \(\beta = 3.2\), corresponding to the best-fitting value for MegaCam observations.

The Moffat model is rasterized on the internal image grid, truncated inside the circular simulation domain, and normalized numerically:

(34)\[\phi_i = \frac{M(r_i)} {\sum_j M(r_j)},\]

giving the final rasterized point-source image model.

Galaxy model

A galaxy is modeled as a circular Sérsic profile [4, 5, 6] convolved with the seeing PSF. Before convolution, the intrinsic surface-brightness profile is

(35)\[S(r) = \exp\left\{-b_n\left[\left(\frac{r}{R_\mathrm{e}}\right)^{1/n} - 1\right]\right\},\]

where \(R_\mathrm{e}\) is the effective radius and \(n\) the Sérsic index.

The effective radius \(R_\mathrm{e}\) is the radius enclosing half of the intrinsic Sérsic-model flux. The Sérsic index controls the concentration of the profile. Values close to \(n=1\) correspond to disk-like profiles, whereas larger values produce more centrally concentrated profiles with more extended wings, such as in ellipticals or galaxy bulges. For numerical stability, the current implementation limits the Sérsic index to \(0.36 \leq n \leq 10\).

The observed galaxy model is obtained by convolving the intrinsic Sérsic profile with the PSF model,

(36)\[G(r) = (S * P)(r),\]

where \(*\) denotes convolution. The convolution is performed on the rasterized model using Fourier transforms.

The final galaxy image model is then truncated inside the circular simulation domain and normalized numerically:

(37)\[\phi_i = \frac{G_i} {\sum_j G_j}.\]

The total source rate \(f\) is therefore the total detected photon rate of the galaxy.

Flat extended-source model

The extended-source model represents a flat source with constant surface brightness. It is not normalized to unit total flux, because the source is not treated as a finite object. Instead, the input brightness is interpreted per unit angular area.

Let \(a_i\) be the angular area of pixel \(i\), in arcsec². For a flat source with detected surface photon rate \(f_\Omega\), in photons s^-1 arcsec^-2, the expected source contribution in pixel \(i\) during an exposure time \(t\) is

(38)\[f_\Omega . t . a_i.\]

For a regular detector grid, all pixels have the same angular area \(a_\mathrm{pix}\). The model image is therefore simply rasterized as

(39)\[\phi_i = a_\mathrm{pix}\]

inside a large disk, and zero outside it. Convolving a constant surface-brightness field by a normalized PSF gives the same constant surface brightness. Therefore, the seeing parameter has no direct effect on the mathematical model of an extended (flat) source.

Synthetic photometry

PyDIET converts source brightnesses and sky surface brightnesses into detected count rates before evaluating the signal-to-noise ratio. This conversion is performed by integrating reference spectra or sky spectra through the selected instrumental and atmospheric responses.

Transmission

Let \(T_\mathrm{inst}(\lambda)\) be the instrumental response, including the selected mirror state, telescope and instrument optics, filter transmission, and detector response. Let \(T_\mathrm{atm}(\lambda; X)\) be the atmospheric transmission at airmass \(X\), and \(\tau\) the additional transparency factor entered by the user. The full source transmission used for source count rates is

(40)\[T_\mathrm{full}(\lambda) = T_\mathrm{inst}(\lambda) T_\mathrm{atm}(\lambda; X) \tau .\]

The effective collecting area is

(41)\[A_\mathrm{eff} = A_\mathrm{tel} - A_\mathrm{obs},\]

where \(A_\mathrm{tel}\) is the telescope collecting area and \(A_\mathrm{obs}\) is the obstructed area, which is often instrument-dependent (especially for prime-focus cameras).

Magnitude zero-point

For a given filter and instrument configuration, PyDIET computes a photometric zero-point from the count rate of a reference source.

For AB magnitudes, the reference spectrum is the constant-flux-density \(m_\mathrm{AB}=0\) spectrum. For Vega magnitudes, the reference spectrum is the Vega spectrum. In both cases, the reference count rate is computed by integrating the reference spectrum through the relevant transmission curve.

For the full atmospheric response, the reference count rate is

(42)\[ f_0 = \int R_{0,\lambda}(\lambda) T_\mathrm{full}(\lambda) A_\mathrm{eff}{\lambda \over hc}\,d\lambda ,\]

where \(R_{0,\lambda}\) is the chosen reference spectrum. PyDIET evaluates this integral using countrate() from the Synphot observation module.

If the detector gain is \(g\), in electrons per ADU, the corresponding count rate in ADU per second is

(43)\[ C_0 = {f_0 \over g}.\]

The magnitude zero-point is then

(44)\[ \mathrm{ZP} = 2.5 \log_{10} C_0 .\]

Equivalently, a source producing a measured count rate \(C\), in ADU per second, has magnitude

(45)\[ m = \mathrm{ZP} - 2.5 \log_{10} C .\]

PyDIET reports two related zero-points:

• The full AB zero-point, including atmosphere and transparency. This is the zero-point corresponding to \(T_\mathrm{full}\).

• The instrumental AB zero-point (accessible through the web API), computed from \(T_\mathrm{inst}\) only, without atmospheric transmission. This value is cached for each mirror/filter configuration and is used in particular to express the sky background as an AB surface brightness.

For a source magnitude \(m\), the detected source count rate is

(46)\[ f = f_0 10^{-0.4m} .\]

For flux-density units, PyDIET converts the input value to the equivalent multiplicative factor relative to the AB reference spectrum. For integrated physical fluxes, the conversion uses the pivot wavelength and equivalent rectangular bandwidth of the full transmission curve.

Sky background

The sky background enters the SNR model as a count rate per detector pixel. PyDIET supports two cases:

• a predefined sky model selected from tabulated sky spectra;

• a user-specified sky brightness.

For the predefined sky modes, PyDIET selects a tabulated sky-emission spectrum according to:

• the sky state, for example dark, grey, or bright;

• the solar-activity level;

• the observing airmass.

The sky spectrum is interpolated linearly between the available tabulated airmass values. The result is a sky surface-brightness spectrum,

(47)\[ I_\mathrm{sky}(\lambda; X, s, q),\]

where \(X\) is the airmass, \(s\) the sky state, and \(q\) the solar-activity state.

This sky-emission spectrum is already a surface brightness at the observing site. It is therefore propagated only through the instrumental response, leaving out the transmission of the atmosphere. The sky count rate per square arcsecond is

(48)\[ B_\Omega = \int I_\mathrm{sky}(\lambda) T_\mathrm{inst}(\lambda) A_\mathrm{eff} {\lambda \over hc}\, d\lambda + B_\mathrm{therm} .\]

The term \(B_\mathrm{therm}\) is the instrumental thermal-emission contribution for the selected mirror and filter configuration. The current calculation includes thermal emission from the telescope, optics, and filter; detector emission is not included.

When the sky background is specified manually, PyDIET converts the user input into a detector count rate per pixel.

PyDIET reports the sky background both as a count rate per pixel (accessible through the web API) and as an AB surface brightness.

If the detector pixel scale is \(p_x \times p_y\), in arcseconds per pixel, the value passed to the image model as the background rate per pixel is then

(49)\[ B_\mathrm{pix} = B_\Omega p_x p_y .\]

The reported AB sky surface brightness is derived from the instrumental AB zero-point:

(50)\[ \mu_\mathrm{sky,AB} = \mathrm{ZP}_\mathrm{inst,AB} - 2.5 \log_{10}\left( {B_\Omega \over g}\right).\]

The instrumental zero-point is used here because the sky-emission spectrum is already the sky brightness at the observing site; it should not be attenuated by the atmosphere a second time.