Photometry

Consider techniques for stellar photometry, surface photometry, calibration of broad band and emission line photometry.

In aperture photometry, we simply identify stars and add up all of the light in the surrounding pixels. The light is spread over several pixels in a form given by the point spread function PSF(i, j) . We must also measure the background and subtract its contribution to the sum.

Error estimation: we're going to sum over pixels. First just consider the object aperture:

S = $\displaystyle \sum$ [Star(i, j) + B(i, j)] - NB

$\displaystyle \sigma_{S}^{2}$ = $\displaystyle \sum$ (Star(i, j) + B(i, j) + rn²) + N² $\displaystyle \sigma_{B}^{2}$

$\displaystyle \sigma_{S}^{2}$ = S + $\displaystyle \sum$ (B(i, j) + rn²) = S + N(B + rn²) + N² $\displaystyle \sigma_{B}^{2}$

Now we have to consider how we will estimate B . In the simplest case, we just go away from the star and take then mean level. Generally, we want to minimize effects of any nonuniform background, so we consider an annulus around the star. If this annulus contains N_a pixels, then

$\displaystyle \sigma_{B}^{2}$ = $\displaystyle {\sum [B(i,j) + rn^2]\over N_a^2}$

S/N = $\displaystyle {S\over \sqrt{S + N(B+rn^2)}}$

Using this, one can consider the optimal choice of aperture. As the aperture size increases, you get a increase in total signal, but also an increase in background, so the change in S/N depends on how bright the star is relative to background. The optimal choice of aperture will depend on the brightness of the star: a larger aperture will be better for brighter stars, a smaller aperture for fainter stars. But also recall that we don't want to use a too small aperture or else we won't be able to compare results from different frames because of changes in the PSF.

This leads to the commonly used technique of using small apertures for all of the stars on the frame, but large apertures for a few bright stars. The few bright stars are assumed to have representative PSFs for this frame, so all of the small aperture measurements can be aperture corrected to large aperture measurements without the increase in noise you would get if you actually used a large aperture. Note that you can't go arbitrarily small as you will eventually run into the problem of PSF variations across the frame if for no other reason than any small dependences on pixel centering.

In fact, you can do better for S/N if you use additional information. If you know the shape of the PSF, you can use this information to fit your stellar image, increasing the S/N in the process. Simple linear least squares argument, if you know PSF and position accurately, leads to

$\displaystyle \sigma_{S}^{}$ = S + N_eff(B + rn²)

N_eff = $\displaystyle {1\over \sum PSF(i,j)^2}$

Note that you'll improve S/N, but only if your assumption that your knowledge of the PSF is good is valid. This naturally leads into the next area of stellar photometry in crowding fields, when you are forced into fitting the PSF whether you like it or not.

It is customary to express the observed number of counts in a frame in terms of an instrumental magnitude:

Crowding

Clearly, the above technique breaks down as you have more than a few stars in your frame for two reasons: the stars may have overlapping light, and there may be stars in your sky annuli. This leads to techniques for crowded field photometry, a complicated subject which we'll just review quickly.

In crowded field photometry, the idea to that you have to solve for the brightnesses of overlapping stars while considering the contribution of neighboring stars. To do this requires that you have information about the point spread function. There are several possibilities for how you might consider accounting for neighbors:

Simply, the technique consists of: finding stars, finding the PSF, grouping the stars, and simultaneously fitting for stellar brightnesses and positions. You almost certainly have to do positions because your initial estimates will probably be biased by neighbors. You also might consider fitting for the background as well. We'll consider each one of the steps in order:

Need to consider automation because of completeness issues, not to mention tedium.

Look for peaks; clearly need to consider background noise, so look for peaks above some noise threshold. More sublety: look for peaks which look like they have the right shape to be stars. Matched detector algorithm. Shape parameters, e.g. sharpness, roundness for additional spurious object rejection.

Find bright isolated stars. Tabulate the PSF. However, rememeber that you're going to have to use this PSF to estimate brightnesses for other stars, so you'll need to be able to interpolate it accurately. Consider using a functional fit to PSF, possibly carry along residuals as well.

As discussed above, either consider all stars that overlap, or iterate each star individually with improved neighbor subtraction at each iteration.

For all pixels under consideration, compare observed values with first guess values. Use residuals to refine your guesses (nonlinear least squares). Continue iterating until parameters converge.

Multiple iterations of entire procedure to improve PSF by subtracting neighbors, also to find new stars under wings of other stars.

Averaging measurements. Averaging mags not the same as averaging counts. Beware of weighted means if you use estimated errors, because these can be biased.

Atmospheric extinction

If you want to compare brightnesses of stars in different observations, you need to correct for the different absorption of the earth's atmosphere as a function of airmass. Note that if you are just comparing brightnesses within individual frames (differential photometry), correcting for atmospheric absorption isn't necessary, to the extent that the absorption is the same for all objects in the frame.

The extinction coefficient can be determined by making multiple observations of a star at different airmasses. Then you can solve for k and m₀ using least squares. Note that you need to sample a good range of airmasses to get good leverage on the fit, and you should bracket the airmasses of all of your program objects.

When doing broad-band work, however, there's an additional subtlety because of the wavelength dependence. Two different stars might have different extinction coefficients in the same filter because the stars might have very different colors, which has an effect if the filter bandpass is broad. This problem can be solved by using second order extinction coefficients, where the extinction is a function not only of the airmass, but also of the stellar color:

Another important thing to remember is that although we correct for the effects of the Earth's atmosphere using extinction coefficients, the actual response function of an imaging system includes the wavelength-dependent response of the Earth's atmosphere, which varies depending on airmass. If one wants to determine an ``average'' wavelength response of a system, one probably wants to include the effects of the atmosphere at some ``typical'' airmass, and then correct the photometry to this airmass; for example, this is what is done by the SDSS photometric pipeline, where all photometry is corrected to an airmass of 1.3.

Using the above formalism, you have to solve for extinction coefficients plus a magnitude for every star you observe. This requires a fair bit of observing to constrain all of the parameters. Clearly, if you can observe stars of known brightnesses, you will have better constraints on the extinction coefficients. These leads us into the discussion of standard stars.

Standard stars

Standard stars are required so that different observers are able to compare results with each other. The reason this is true is because every observational setup is likely to have different response functions, so the same stars will not be observed to have the same brightnesses (even relative brightnesses!) with each separate setup. Differences in response come from many factors: size and condition of the telescope optics, number and type of optics in the system, bandpass and quality of the filter, response function of the CCD, etc. In addition, it is very difficult to measure all of these responses in any absolute sense, so it is almost impossible to go straight from observations to inferences about physical fluxes in absolute units.

To get around these problems, systems of standard stars have been set up so observers can calibrate any new observations against the known brightnesses of the standard stars. A standard system is just a more-or-less arbitrary choice of one particular instrumental setup, but is generally one which is very stable and which someone has spent a lot of observing time establishing. Generally, the standard system is tied to absolute flux observations of some primary standard star so that observations on the standard systems can be converted to physical fluxes.

Let's consider some examples. First consider two identical systems. Clearly these will give the same observed flux as each other. Now let the systems differ only by a multiplicative factor in throughput. Now the observed fluxes will differ by a constant factor, or a constant additive factor in magnitudes. With observations of standard stars, it is straightforward to find the magnitude difference between the observed magnitude and the standard magnitude to calibrate out this constant factor. The additive constant in magnitudes is called the zeropoint.

A slight divergence is called for here to mention that different packages apply some ``software'' zeropoints to the calculation of instrumental magnitude just to make the numbers (instrumental magnitudes) look ``reasonable'' (i.e., familiar). This is no problem since we determine a zeropoint on top of this to calibrate, but caution is required if talking about instrumental magnitudes measured with different software systems!

OK, back to differing instrumental systems. Now let's consider the more realistic case where the shape of the response curve as well as the absolute throughput differs from the standard system. Let's consider a system which has slightly different filters with a different wavelength cutoff. Now you observe a different number of counts from the standard system, and the difference depends on the spectrum of the object you are looking at. Now spectra can differ for lots of reasons, but the biggest effect is just from spectral slope differences. If your filter has a shorter wavelength cutoff (on the red side), then you'll observe less counts than the standard system, and preferentially less counts for a redder star. To first order, this can be calibrated out by solving for an additive constant which depends on the spectrum of the object being looked at. Since we don't in general know the spectrum, we have to parameterize it by something we can observe, and the best choice here is the stellar color, as inferred from observations taken through more than one filter.

The color is generally parameterized by the ratio of the flux at two different wavelengths, or, in magnitudes, the difference between the magnitudes. The two wavelengths should be measured near in wavelength to the wavelength of the filter being corrected; generally, one uses the bandpass being corrected as one of the wavelenghts and an adjacent bandpass as the other. For example, when correcting V magnitudes, people usually use B - V , V - R , or V - I for the color term.

There are two ways to define the color, either in terms of the observational system or in terms of the standard system. The latter is slightly preferred for using least-sqaures (small errors on the independent variable), and also because it allows observations from different nights to be combined. Note that this formulation does not require you to know the colors of your objects a priori, it's just algebra to figure them out as long as you have observations in both filters.

The use of these first-order transformation coefficients is accurate as long as your filter system does not differ much from the standard system, and additionally, that the spectrum of your program objects does not differ significantly from the spectrum of the standard objects. The more these conditions are not met, the less accurate the results. Some additional accuracy in the case of differing systems can be achieved by using higher order transformation coefficients. However, even in this case, it is always important to remember that if the spectrum of the program object differs significantly from the standards, derived fluxes can be significantly in error.

Certainly, you get to a point when the response of one system is so different than the response of another system that no transformation can be determined. In this case, you have two different photometric systems. In fact, there are several different photometric systems at use in astronomy today, and each has advantages and disadvantages.

It is common practice to combine the equations for extinction and transformation coefficients into a single set of equations:

In practice, one observes a set of standards of different colors at different airmasses. Then one uses least squares techniques to solve for the values of k , t , and z which minimize the difference between the observed, transformed, magnitudes and the standard magnitudes. Then one applies these coefficients to measurements of your program objects to derive their magnitudes. Clearly, to do so, you must measure the colors of your program objects by making observations in more than one bandpass since, unlike the standard stars, the colors of your program objects are not known a priori.

Note that for calibration of broad-band observations, we do not require knowledge of the transmission function. Of course, if you have such knowledge, you could derive synthetic transformations, if you also know the transmission functions of the standard system. Beware, however, that getting accurate tranmssion functions can be very difficult.

Least squares fitting

We have encountered the concept of least-squares fitting several times: once to derive extinction and transformation coefficients and zero points, and again when we talked about crowded field photometry. Least squares fitting has a very large number of applications in science in general, as it is a commonly used technique to determine the parameters of some model which one wishes to fit to their data.

The concept of least-squares fitting has a statistical basis when one considers a set of observations which arise from a system which can be described by some mathematical model but where the observations have errors which are normally distributed around the predicted model values. In this case, we can write the probability of observing any particular data point given a model, y(x) :

P(y_i) $\displaystyle \propto$ exp $\displaystyle {-0.5 (y_i - y(x_i))^2\over \sigma_i^2}$

P $\displaystyle \propto$ $\displaystyle \Pi_{{i=1}}^{N}$ exp $\displaystyle {-0.5 (y_i - y(x_i))^2\over \sigma_i^2}$

$\displaystyle \chi^{2}_{}$ $\displaystyle \equiv$ $\displaystyle \sum_{{i=1}}^{N}$ $\displaystyle {(y_i - y(x_i))^2\over \sigma_i^2}$

This is the standard least-squares equation: one wishes to minimize the sum of the square of the differences between observation and model by adjusting the model parameters. We define this quantity as ``chi-squared'', or $\chi^{2}_{}$ . For a good fit, one would expect that on average, points deviate from the model by roughly $\sigma_{i}^{}$ , so one would expect a $\chi^{2}_{}$ for a set of N observations to approach N if the model is a good model and the errors have been estimated properly. In fact, it is possible given an observed value of $\chi^{2}_{}$ to compute the probability that this value would be obtained if the model is correct; this allows one to judge the quality or likelihood that the model which minimizes $\chi^{2}_{}$ is actually the correct model. One often discussed the reduced $\chi^{2}_{}$ quantity, $\chi_{\nu}^{2}$ , which is just $\chi^{2}_{}$ divided by N + M , where N is the total number of points, and M is the number of free parameters in the fits; the latter is there because any set of observations which does not have more data points than the number of free parameters can generally be fit perfectly even for an incorrect model. The total N + M is called the degrees-of-freedom of the fit.

How does one do least-squares in practice? Basically, one wishes to minimize $\chi^{2}_{}$ with respect to one or more parameters of a model. Let us first consider models which can be written in the form:

y(x_i) = $\displaystyle \sum_{{i=1}}^{M}$ a_kX_k(x_i)

$\displaystyle \chi^{2}_{}$ = $\displaystyle \sum_{{i=1}}^{N}$ y_i - $\displaystyle \sum_{{k=1}}^{M}$ a_kX_k(x_i) $\displaystyle \sigma_{i}^{}$ $\displaystyle \left.\vphantom{ y_i - \sum_{k=1}^M a_k X_k(x_i)\over \sigma_i}\right]^{2}_{}$

$\displaystyle \sum_{{i=1}}^{N}$ 2y_i - $\displaystyle \sum_{{k=1}}^{M}$ a_kX_k(x_i) $\displaystyle \sigma_{i}^{}$ $\displaystyle \left.\vphantom{ y_i - \sum_{k=1}^M a_k X_k(x_i)\over \sigma_i}\right]$ (- X_m(x_i)) = 0

$\displaystyle \sum_{{i=1}}^{N}$ $\displaystyle {X_m \sum_{k=1}^M a_k X_k(x_i) \over \sigma_i^2}$ = $\displaystyle \sum_{{i=1}}^{N}$ y_iX_m(x_i)

$\displaystyle \alpha_{{km}}^{}$ $\displaystyle \equiv$ $\displaystyle \sum_{{i=1}}^{N}$ $\displaystyle {X_m(x_i) X_k(x_i) \over \sigma_i^2}$

$\displaystyle \beta_{{k}}^{}$ $\displaystyle \equiv$ $\displaystyle \sum_{{i=1}}^{N}$ $\displaystyle {y_i X_k(x_i) \over \sigma_i^2}$

$\displaystyle \alpha_{{km}}^{}$ a_m = $\displaystyle \beta_{k}^{}$

a_m = $\displaystyle \alpha^{{-1}}_{}$ $\displaystyle \beta$

Associated with solving for the ``best'' parameters, one often wishes to compute the errors associated with the fit parameters. This is discussed in Numerical Recipes (chapter 15), with the result:

$\displaystyle \sigma^{2}_{}$ (a_m) = $\displaystyle \alpha^{{-1}}_{{mm}}$

One can also consider the situation where a model is non-linear in the parameters, a_m , i.e., the derivatives with respect to each parameter may depend on the parameter value. This leads to non-linear least squares techniques. These are more complex, as one can have situations where there are many minima in $\chi^{2}_{}$ and one needs to find the global minimum rather than a local one. Such problems require a starting guess of a reasonable solution and then iteration towards the best solution. The crowded field photometry problem falls into this category because the model is nonlinear in the position parameters.

Photometric systems

What are some of the considerations for choosing a photometric system? One general approach is to consider a broadband system to give a low resolution approximation to the spectra of objects, regardless of the details of what objects you might be looking at.

For broadband systems, the most common has been the UBVRI system, originally established by Johnson in the 1960s. However, even the UBVRI system has undergone some revisions, and Johnson's original system is not quite the one used predominantly today.

Choice of photometric system for studies of stars: determination of temperatures, reddenings, surface gravities, metallicites.

What goes into a choice of filters to use often involves a tradeoff of better spectral information vs. S/N considerations. Note that S/N considerations must include both the errors in the bandpasses being considered, and also the degree to which different objects are separated in color space. (possible problem to determine best filter choice given some sensitivities, astrophysical colors of stars).

Emission line photometry

Emission line objects: often emission line objects have underlying continuum. Distinction about what you want to measure: fluxes vs. equivalent widths. In either case need to measure continuum. Use of neighboring filters (beware of neighboring lines) and/or multiple filter widths. Optimal choice of filter widths and exposure time ratios.

Calibration issues: for emission line sources, total filter transmission isn't relevant quantity, rather transmission at line location is. Two possibilities for calibration:

Surface photometry

What if we have extended (resolved) objects? Can talk about either integrated or surface brightnesses. If integrated brightnesses, use some sort of aperture, remember aperture corrections (though you probably can't estimate them very well, and they probably don't matter too much).

For surface brightnesses, just look at fluxes per area on sky. Remember that there will be some mixing from PSF/seeing; topic of extracting most spatial information will be discussed later. Counts/pixel go directly to mag/square arcsec using stellar calibration, size of pixels (beware of distortion).

If object is irregular, not much more you can do. Hard to observe faint irregular objects accurately because of S/N. However, if object has some degree of regularity (e.g., galaxies), one can average over regions to increase S/N substantially. Generally, many galaxies can be fairly well parameterized by elliptical isophotes, so if you can determine ellipse parameters, you can average along a given isophote to increase S/N by $\sim$ $\sqrt{{N_{pix}}}$ .

Note that if all isophotes are concentric and have same ellipticity, one could use elliptical aperture photometry. But many galaxies have twisting isophotes, so often one needs to solve for varying ellipticity and position angle as a function of semimajor axis.

Several methods for doing this have been presented by Kent (ApJ 266, 562) and Lauer (ApJ 311, 34). They are similar in that they solve for elliptical isophotes using nonlinear least squares. Both start with some a priori guess of isophote parameters. Kent method solves for Fourier coefficients which parametrize (observed-guess) ellipse; ellipse parameters modified to minimize A₁, B₁, A₂, B₂ :

I = I₀ + A₁cosE + B₁sinE + A₂cos2E + B₂sin2E + ^...

y = a₀sinE(1 - $\displaystyle \epsilon_{0}^{}$ )

Lauer describes galaxy as having power law surface brightness profile with constant ellipticity/position angle between some a priori specified semimajor axes.

G(x, y) = $\displaystyle \mu_{k}^{}$ ( $\displaystyle {r\over d_k}$ )

$\displaystyle \gamma$ = log( $\displaystyle {\mu_{k+1}\over \mu_k}$ )/log( $\displaystyle {d_{k+1}\over d_k}$ )

d_k = a_k(1 - $\displaystyle \epsilon_{k}^{}$ )[sin²( $\displaystyle \phi_{k}^{}$ - $\displaystyle \theta_{k}^{}$ ) + (1 - $\displaystyle \epsilon_{k}^{}$ )²cos²( $\displaystyle \phi_{k}^{}$ - $\displaystyle \theta_{k}^{}$ )]^-1/2

Reminder: surface photometry, integrated photometry and variable pixel area effects.