We have an ensemble of images f_i(\bm{x}) covering the same region, indexed by i. The PSF for the images are denoted as p_i(\bm{x}). The image before the PSF convolution is denoted as I(\bm{x}). In Fourier space, the observed image can be written as \tilde{f}_i(\bm{k})= \tilde{I}(\bm{k}) \tilde{p}_i(\bm{k}) + \tilde{n}_i(\bm{k})\,, where \tilde{n}(\bm{k}) is the Fourier transform of noise. We assume that the oversampled images (f_i) are registered and placed on the same grids. With the assumption of oversampling, we write signal as continuous function. In addition, we assume noise is homogeneous and uncorrelated: \langle n_i(\bm{x}) n_i(\bm{x'})\rangle = \sigma_i^2\delta_\text{D}(\bm{x}-\bm{x}'), where \delta_\text{D} is the Dirac delta.
The Kaiser coadding algorithm first smooth each observed image with its own PSF rotated by \pi, weighted by \frac{1}{\sigma_i^2} and sum over the images to create a matched-filtered image that is optimal for source detection: \tilde{S}(\bm{k}) = \sum_i \frac{1}{\sigma_i^2} \tilde{f}_i(\bm{k}) \tilde{p}^*_i(\bm{k}) = \sum_i \frac{1}{\sigma_i^2}\tilde{I}(\bm{k}) |\tilde{p}_i(\bm{k})|^2 + \sum_i \frac{1}{\sigma_i^2}\tilde{n}_i(\bm{k})\tilde{p}^*_i(\bm{k})\,. The noise on the matched-filtered image is \tilde{n}_s (\bm{k}) = \sum_i \frac{1}{\sigma_i^2}\tilde{n}_i(\bm{k})\tilde{p}_i(\bm{k})\,. The power spectrum of the noise field is \langle \tilde{n}^*_s (\bm{k}) \tilde{n}_s (\bm{k}) \rangle = \sum_i \frac{(2\pi)^2}{\sigma_i^2} |\tilde{p}_i(\bm{k})|^2 \,. To make sharp coadded images, the matched-filtered image is deconvolved to whiten the noise: \begin{equation} \begin{split} \notag \tilde{T}(\bm{k}) &= \frac{\sum_i \frac{1}{\sigma_i^2}|\tilde{p}_i(\bm{k})|^2} {\sqrt{\sum_i \frac{1}{\sigma_i^2}|\tilde{p}_i(\bm{k})|^2}} \tilde{I}(\bm{k})\\ &= \sqrt{\sum_i \frac{1}{\sigma_i^2}|\tilde{p}_i(\bm{k})|^2} \, \tilde{I}(\bm{k}) \end{split} \end{equation} Finally, the image can be normalized to the original magnitude zero point \tilde{T}(\bm{k}) = \frac{\sqrt{\sum_i \frac{1}{\sigma_i^2}|\tilde{p}_i(\bm{k})|^2}} {\sqrt{\sum_i \frac{1}{\sigma_i^2}}} \tilde{I}(\bm{k}), and the resulting uncorrelated, homogeneous noise has variance \sigma_T in configuration space, where \begin{equation} \frac{1}{\sigma_T^2} = \sum_i \frac{1}{\sigma_i^2}\,. \end{equation} The final PSF is \begin{equation} \tilde{p_T} = \frac{\sqrt{\sum_i \frac{1}{\sigma_i^2}|\tilde{p}_i(\bm{k})|^2}} {\sqrt{\sum_i \frac{1}{\sigma_i^2}}}\,. \end{equation} Using the PSF and noise on coadd image, we can derive the statistical error on the flux measurement of point source is \begin{equation} \sigma_F^2 = \frac{(2\pi)^2\sum_i \frac{1}{\sigma_i^2} \int \mathrm{d}^2k\,|\tilde{p}_i(\bm{k})|^2} {(\sum_i \frac{1}{\sigma_i^2})^2}\,. \end{equation} Assuming that the flux of the point source is F, we can derive the signal-to-noise ratio (\mathrm{S/N}) of the flux measurement: \begin{equation} \begin{split} \mathrm{S/N} &= F \sqrt{\sum_i \frac{1}{\sigma_i^2 (2\pi)^2} \int \mathrm{d}^2k |\tilde{p}_i(\bm{k})|^2}\\ & = F \sqrt{\sum_i \frac{1}{\sigma_i^2 } \int \mathrm{d}^2x |\tilde{p}_i(\bm{x})|^2}\,. \end{split} \end{equation}