Given a relative flat-fielding error of sigma_flat [%], the relative 
photometric error in total source flux scales as:

sigma(S_tot)/S_tot = sigma_flat * sqrt[ 1 + 2*(B/S)^2 ] / sqrt(Np),

where B/S is the background to source flux ratio on the same spatial scale.
The "2" in front of the (B/S)^2 is due to the error we make when applying
the flat calibration (on pixels containing "S+B"), then another error we make
when estimating the background on a different location to subtract from 
"S+B". You can argue this second error is negligible because we usually 
median or average pixels within a large annulus to get a more accurate 
estimate of B, but it's best to include this error for conservatism.

-----------
Derivation:

Start with estimate of source flux in a single pixel, following flat-fielding
and background subtraction:
S_pix_est = (1 + e1)*(S_pix + B_pix) - (1 + e2)*B_pix,

where:
S_pix, B_pix = per-pixel source and background true signals.
e1 = fractional error in flatfield per-pixel on top of source.
e2 = fractional error in flatfield per-pixel where background
     is measured (another location).

We assume e1,e2 are drawn from the same noise distribution, ~ N(0, var_flat)
where var_flat = (sigma_flat)^2. We also assume that e1 and e2 are uncorrelated.

Simplifying,
S_pix_est = S_pix + e1*S_pix + (e1 + e2)*B_pix.

Take variances where errors in true values are zero; we also neglect Poisson
noise and other noise contributions:
var(S_pix_est) = var(e1)*S_pix^2 + 2*var_flat*B_pix^2,

Divide through by square of source pixel signal (S_pix) and take square root
of both sides:
sigma(S_pix)/S_pix = sqrt[ sigma_flat^2 +  2*sigma_flat^2 * (B_pix/S_pix)^2 ].

Simplifying,
sigma(S_pix)/S_pix = sigma_flat * sqrt[ 1 + 2 * (B_pix/S_pix)^2 ]  ----- (1)
This is the fractional error for a single pixel containing source signal.

Sigma in total source signal (a point source covered by effectively ~ Np pixels):
sigma(S_tot) = sqrt(Np) * sigma(S_pix)  ----- (2)

Total source signal:
S_tot = Np * S_pix  ----- (3)

Combining (2) and (3), the fractional error in total source signal is:
sigma(S_tot)/S_tot = [sigma(S_pix) / S_pix] / sqrt(Np)  ----- (4)

Plugging (1) into (4), we have the final result:
sigma(S_tot)/S_tot  =  sigma_flat * sqrt[ 1 + 2 * (B_pix/S_pix)^2 ] / sqrt(Np).

---
F. Masci
07 / 19 / 21