I guess I'm a bit confused about the definition of sigma_sky...
Let me first explain what I have in mind so we can check if we're
actually saying the same thing or not.
The error from the sky is the sum in quadrature of two contributions:
the first (let's call it e_1) is just the total uncorrelated noise
from all individual pixels in the aperture and the second (e_bkg) is
the uncertainty in the actual average background level in the
aperture. e_sky=sqrt(e_1^2+e_bkg^2)
e_1=sqrt[N_pix*(sigma_1*Omega_pix)^2]=sigma_1*Omega_pix*sqrt(N_pix),
where sigma_1 is just the rms of the SB in individual pixels.
e_bkg could be directly evaluated (in principle) by measuring the
average <SB>_i, i=1..N_ap, within N_ap apertures of the same size of
the galaxy. Let sigma_ap be the rms of the <SB>_i, then
e_bkg=N_pix*Omega_pix*sigma_ap

So when I sum the two contributions I get:
e_sky=Omega_pix*sqrt(N_pix*sigma_1^2+N_pix^2*sigma_ap^2)

One thing I don't understand is whether your sigma_sky is the same as
my sigma_1, as it seems to be the case if I compare the first term in
our two equations respectively.
Moreover, if I make the assumption that sigma_ap=sigma_1/sqrt(N_sky)
(and sigma_1=sigma_sky), my equation becomes:
e_sky=Omega_pix*sqrt(N_pix*sigma_sky^2+N_pix^2*sigma_sky^2/N_sky)=
        =sigma_sky*Omega_pix*sqrt(N_pix+N_pix^2/N_sky)
which is exactly the same as yours.
However, sigma_ap=sigma_1/sqrt(N_sky) only if large scale fluctuations
are negligible. Strictly speaking, this assumption holds only if all
background pixels are gaussian distributed with
sigma=sigma_1=sigma_sky, such that the error on the mean value of
N_sky pixels is sigma_1/sqrt(N_sky). This would imply that by
increasing the number of sky pixels you can reduce this contribution
as much as you like, but in reality you have correlated fluctuations
on the scale of the aperture size. This is a solid limit to the
accuracy with which you can infer the actual level of the background
below the galaxy from a measurement of the surrounding background.
This is reflected in the equation I've written above, where the term
N_pix^2*sigma_ap^2 is independent of the number of sky pixels you use.
Does it make sense to you?
Hope this clarifies things!

*** Thanks for taking the time to teach me these details!  I ultimately decided to go with my original approach since I have used generous values for sigma_sky.