Aerosol Measurements with the FRAM Telescope

Precision stellar photometry using a telescope equipped with a CCD camera is an obvious way to measure the total aerosol content of the atmosphere as the apparent brightness of every star is affected by scattering. Achieving high precision in the vertical aerosol optical depth (at the level of 0.01) presents a series of interesting challenges. Using 3.5 years of data taken by the FRAM instrument at the Pierre Auger Observatory, we have developed a set of methods and tools to overcome most of these challenges. We use a wide-field camera and measure stars over a large span in airmass to eliminate the need for absolute calibration of the instrument. The main issues for data processing include camera calibration, source identification in curved field, catalog deficiencies, automated aperture photometry in rich fields with lens distortion and corrections for star color. In the next step, we model the airmass-dependence of the extinction and subtract the Rayleigh component of scattering, using laboratory measurements of spectral sensitivity of the device. In this contribution, we focus on the caveats and solutions found during the development of the methods, as well as several issues yet to be solved. Finally, future outlooks, such as the possibility for precision measurements of wavelength dependence of the extinction are discussed.


Motivation
The Pierre Auger Observatory [1] is the world's largest detector of ultra-high energy cosmic rays (UHECR).The extensive air showers (EAS) initiated by the UHECR are observed using an array of surface detectors, which register the secondary particles that reach the ground, and 27 optical telescopes which, during sufficiently clear and dark nights, observe the fluorescence light emitted during the passage of the EAS through the atmosphere.The later method is particularly sensitive to the transmission properties of the atmosphere, as showers are observed to a distance of tens of kilometers from the telescopes.To monitor those properties, many devices, such as cloud cameras, elastic and Raman LIDARs, calibration lasers, all-sky cameras and the FRAM telescope are installed on the site.The FRAM telescope consists of narrow-field (23×15 arcmin 2 ) and wide-field (7 • ×7 • ) optical setups comounted on a commercial astronomical mount and operated as a fully autonomous robotic observatory.Currently the wide-field camera is mostly used for atmospheric monitoring [2].
The main purpose of the FRAM telescope in its current setup is part of the so-called rapid monitoring program [3], which focuses on targeted measurements of atmospheric properties triggered by the observation of spee-mail: ebr@fzu.czFull author list: http://www.auger.org/archive/authors_2016_09.htmlcific EAS in the fluorescence telescopes (this method is often called Shoot-the-Shower, or StS, by the Auger Collaboration).The main goal of this program is the observation of showers with anomalous longitudinal light profiles -for each candidate for such shower, it must be first determined whether the measured fluorescence profile is affected by clouds or other inhomogeneities in the atmosphere.For this purpose, the FRAM telescope has an advantage over the laser-based measurements, because it does not emit any light and thus does not interfere with the operation of the fluorescence telescopes.When a trigger is received for an interesting shower, the wide-field camera of the FRAM telescope is used to measure the extinction of starlight along the observed path of the shower on the sky, taking several (depending on the orientation of the shower, but at least six) images across the whole field of the fluorescence telescopes, that is between 30 and 2 degrees in altitude above the horizon [4].In principle, the difference between measured brightness of a star and that from a catalog can be written as where Z is the zeropoint (a calibration constant of the system), A is airmass (amount of air traversed by the light relative to a star in zenith), k is the extinction coefficient (1.086 × the vertical optical depth) and f (...) is in general a correction function of some properties of the star and the optical system.Local deviations from this formula are an indication of the presence of clouds or inhomogeneous aerosol layers.
A curious reader could ask, what are the values of the extinction coefficient and whether this method could be used to determine this value and thus the vertical aerosol optical depth (VAOD).The answer is a qualified yeswhile extracting the value from the fit is easy, reaching a precision good enough that the data become relevant for atmospheric calibration requires significant effort.In this paper, we describe this effort and the current status of the method, focusing on aerosol measurements.Two main things should be kept in mind: firstly that this is still a work in progress and secondly that the method has specific limitations, in particular in the fact that we always measure the integral extinction across the whole atmosphere and cannot determine the height profile of the aerosols.

Photometry and extinction
First we introduce and clarify some quantities.The light of a star reaches the Earth as an apparent flux J in Wm −2 or similar units.The apparent brightness m is defined in logarithmic units -magnitudes -with respect to a reference flux by m = −2.5 log 10 (J/J ref ). ( The physical flux is detected by a detector, in our case a CCD chip, where it is registered as a current and digitized, resulting in an instrumental flux J inst in arbitrary units.The CCD camera is a highly linear detector (after data reduction) and thus for some calibration constant C. For apparent brightness, we get where Z is called the zeropoint and is the primary calibration constant of any photometric device.
In the atmosphere, light suffers exponential extinction where the optical depth τ can be easily translated into a difference in apparent brightness because τ = ln(J space /J ground ) = 0.921Δm.
For a given direction, the optical depth can be related to the vertical optical depth (VOD) by where the dimensionless quantity A is called airmass and quantifies the relative amount of air traversed by the light with respect to a star at the zenith.For a flat Earth it would have a 1/cos dependence on zenith angle θ, for curved Earth the following approximation is valid up to 89 • : We further define the extinction coefficient k so that it corresponds to the change of apparent brightness by one VOD For a given star, we thus expect the observed magnitude to be related to the catalog one as As the optical depth and extinction coefficient are numerically very close, we use them somewhat interchangeably in the text, however taking care to include the correct factor in actual calculations.In reality, the extinction has several components.In the blue visible light where FRAM mostly operates, the dominant contribution is the Rayleigh scattering off air molecules, while the molecular absorption is negligible (Fig. 1).When the Rayleigh contribution is subtracted from the VOD, the rest is thus assumed to be due to aerosols and is referred to as the vertical aerosol optical depth, VAOD.
In principle, we could extract the extinction coefficient from a single observation of a single star, providing that we know the zeropoint precisely enough.This was the idea of the photoelectric photometer installed on FRAM until 2010, but the technical problems related to its robotic use proved insurmountable.For CCD photometry, such approach is not feasible because of the variations in the zeropoint due to many effects, if one wants to achieve precision higher than of the order 0.1 in optical depth, which is quite coarse considering that 0.1 is the usual VAOD cut in analyses at the Auger Observatory.Therefore, our goal must be to measure the VAOD at the level of at least 0.01.
To do so requires to measure m inst − m cat as a function of A and extract both k and Z simultaneously.We first tried to follow some reference stars during the night as they naturally pass through a range of airmasses, but we found out that probably both k and Z vary considerably during the time needed.Moreover this approach makes it difficult to deal with systematic effects caused by star color and its inevitably changing position on the CCD chip.This had led us to the wide-field approach, where we measure a lot of stars in different airmasses within a short time window and we deal with systematics and outliers within a complex fitting procedure.

Data selection
As we already mentioned, the motivation for the aerosol measurements came from the rapid monitoring program.In this program, we take pictures along the observed sky path of an interesting shower.To fully control the atmospheric influence, we focus not only on the directions where the shower was observed in the fluoresence telescopes, but scan the whole path of the shower across their field of view to make sure that the non-detection of the shower in some parts is not due to a cloud.Thus every shower triggers a scan in altitude covering a swath of sky between 2 and 30 degrees above the horizon, which is potentially a huge span of airmasses -the airmass at 30 degrees altitude is 2 whereas it is roughly 40 at horizon.The useful span is however limited by our ability to detect stars dimmed by very high extinction and also the fact that near to the horizon we are looking to very large distances, increasing the probability to encounter a distant cloud.Such altitude scans provide a quite homogeneous dataset for aerosol measurements but are limited to periods when the rapid monitoring program was working (and interesting showers were detected).Since the beginning of 2016, we started to take additional altitude scans at frequent interval during all nights (even during full Moon when the fluorescence telescopes do not operate at all).
Only those scans that do not show any cloud contamination at all (apart from very low altitudes) are useful for the aerosol analysis, because the presence of clouds makes it very difficult to fit the airmass dependence of extinction.The present work is based on scans taken between March 2013 and April 2016, from which 657 clear scans were selected.These were picked by hand using the plots generated for the rapid monitoring program and then further cleared taking into account the quality of fit, removing any suspicious artifacts and making sure that the camera temperature was properly set for all the images.Since we started the dedicated observations at the beginning of 2016, there is a wealth of new data (likely at least doubling the dataset) awaiting analysis.All the data shown here are based on images taken in the B filter with 30-second exposures with the wide-field camera of the FRAM telescope, which is a Moravian Instruments G4-16000 CCD with a 36mm × 36 mm chip mounted on a Nikkor 300/2.8lens.This setup provides a rectangular field of view 7 • × 7 • .

Data reduction
The CCD chip converts incoming photons to electrons and those are then read out for each pixel and the signal is amplified and digitized in an A/D converter.Thus the signal for each pixel (in ADU) is in principle linearly dependent on the incoming light into the pixel.In reality, the signal is a combination of the light of the star, the night-sky background and other scattered light, the dark current of the chip, and the noise of the electronics.The electronics contribution is very small for Moravian Instruments cameras and also the dark current (the signal recorder with zero light, mostly the thermal current of the semiconducting chip and thus dependent on temperature) is also very small for most of the pixels when cooled to −20 • , with the exception of typically a thousand of "hot" pixels that can even become saturated by the thermal current when the exposure is sufficiently long.The dark current can be dealt with by taking a series of exposures with the shutter closed ("dark frames") and subtracting the average or median of those exposures from the light images, or, if it is sufficiently low, by detecting the hot pixels in the light images and removing them.
Most of the background light is removed by the aperture photometry method (see later), but the dark-current correction plays an important role when we take into account that the sensitivity of the lens-camera system is not uniform across the field of view, because of the different sensitivity of each pixel (a small effect), the vigneting of the lens (a large effect, up to 50 %, Fig. 2) and due to dust and other contamination on various optical surfaces.Normally this is dealt with using the "flat field", an image of a uniformly lit surface.Thus for a CCD image the following holds schematically and thus one would like to retrieve the original image from the one taken by We observe that the dark current slowly increases during the night, even when the temperature is reported to be constant (probably due to the fact that the sensor is not directly on the chip and the chip heats up when used).A single dark image is noisy and thus many have to be taken to create a master dark image used in calibration.This is normally done in the evening or in the morning, but such master dark image either undercompensates or overcompensates the dark current and that affects the flat subtraction.While this is a small effect, we are currently developing methods to deal with it, such as interpolation between evening and morning dark frames.The flat field correction is never perfect anyway, because for a wide-field camera, the evenly lit surface is hard to realize and thus we use the dawn and dusk sky, which has its own drawbacks.We correct for the flat-field deficiencies that are mostly radial, using the data.
After the dark and flat corrections, the flux for each star is evaluated using aperture photometry.In this method, a circle of a fixed size is centered on each detected light source in the image and the light in all the pixels is summed.Another larger circle is then taken concentrically and the signal in the annulus is used to estimate the background to be subtracted from the previous value (Fig. 3).The aperture needs to be large enough so that most of the point spread function is contained in the aperture -taking an aperture too small leads to nonlinearity of the measurement, because for bright stars the amount of light lost outside the aperture becomes more significant.On the other hand, taking the aperture too large means a danger of contamination with the light of another nearby star.We avoid this problem by rejecting any pair of detected sources that is closer together than twice the aperture, but this leads to significant losses of observable stars in some of in star-rich areas of the Milky Way.
Another interesting issue with aperture photometry arises when the shapes of the images of stars are not constant across the field of view.The coma and astigmatism that appear for stars farther from the center of the image are actually not a big issue, because their dependence on the position on the chip is radial and thus mimic the effect of the vignetting which we correct for using the data anyway.A bigger problem was found by visually inspecting the images, when it became apparent that the stars have the tendency to be less sharp at one side of the image.We attribute this behavior to mechanical flexing between the lens and the camera and have since improved the mechanical support accordingly.We have found that the increase in the FWHM of the stars from one side of the image to the other correlates with the decrease of their observed brightness (Fig. 4).This effect cannot be corrected for in the data because it can appear in the same direction as does the atmospheric extinction that we are trying to fit.

Fitting stars with a model
When the light flux from the stars has been measured, we need to identify those in the catalog.The wide field of view of the camera is distorted and thus a simple linear identification algorithm fails.While tools to identify stars in curved fields in principle exist, we were not able to make them work properly for our data.Instead of writing our own method, we have resorted to a simpler solution: we cut the images into smaller sub-images and within those the distortion is small enough so that the linear identification works.We have tested various configuration, monitoring the distances between the measured and catalog sources and we have concluded that cutting the image based on a 4 × 4 grid of equally sized tiles and merging the middle 2 × 2 into one works very well (Fig. 5).This method is slightly problematic in cloudy weather, when some of the tiles may not have enough stars to properly identify the star fields, but that is not relevant to the aerosol analysis.
The comparison with the catalog is further complicated by different star colors.In astronomy, it is customary to measure star brightness in one of pre-defined "filters", which nowadays mean a complete specification of a bandpass.The most widespread set of filters is the Johnson UBVRI system, to which FRAM with its physical color filters is roughly compliant.However the Tycho2 catalog that we use (as the only homogeneous all-sky catalog for the relevant range of star magnitudes available at the time when we started the analysis) uses B and V filters that are significantly different from the Johnson standard.We fit the dependence of the brightness observed by FRAM on the catalog values given by Tycho using our data and we have found that B FRAM ≈ 0.7B Tycho + 0.3V Tycho (13) where this formula holds only for stars with B − V < 0.8, because at higher values of B − V (more red stars) another population appears.This is likely a real physical effect of star populations in the Galaxy.We are currently investigating the option to use a new all-sky photometric catalog APASS from AAVSO which may make the color calibration easier and also allow us to extend the measurements to the R filter where the Tycho2 catalog has no data.We will come to the problems caused by the stellar spectra later.Putting all the effects we have mentioned together, we can formulate the following model for the observed brightness of a star: where B and V are magnitudes of the star in the catalog in the respective filters, r is the distance of the star from the center of the respective image and M, c 1 , c 2 , R 1 , R 2 are parameters of the system (optics+CCD+data reduction).
The last term represents the dependence of the extinction on the star color (and we will revisit it later).
We fit this model to many scans at once.For each scan (which comprises several images) we allow a different pair of Z and k (thus the index i) while the other parameters are the same for all scans -even k C , because we expect the color-dependence to be dominated by the Rayleigh scattering, which is constant.We take only stars brighter than 9.5 mag because of the unreliability of the Tycho2 catalog below this value and only stars fainter than 6.5 mag because of saturation.The time needed for the multi-dimensional fit increases with the number of parameters; we have found that we can comfortably fit about 200 scans at once.We perform this procedure in two iterations, cutting obvious outliers after the first iteration.We have tested that the multi-dimensional fitting algorithm is reliable by fitting only Z and k shower by shower using the system parameters from the global fit.The parameter M comes out very close to one, showing that our assertions of linearity are justified (this however only works for apertures that are much larger than how big the stars look at the image).
The resulting fit has RMS of residua 0.07 magnitudes (Fig. 6) and there are no apparent trends when the scans are cut at A < 8 (we can now do it for even higher airmasses, thanks to a correction that we discuss later).The zeropoint obviously changes with time and there was a period at the beginning when it was rather fluctuating -we believe that a large effect on the zeropoint changes is in the focusing (again, due to the aperture photometry) and we are planning to implement a temperature compensation for the focusing, which we have already tested on the narrow-field camera.The time series of zeropoints clearly shows the time when the optics was cleaned, as the environment is very dusty.However the effects of the cleaning were reversed rather quickly (Fig. 7).

Rayleigh subtraction
To isolate the aerosol information, we must subtract the effects of the molecular atmosphere from the measured extinction.As we have noted, molecular absorption is negligible in the B filter (but very important in the V and R fil-ters) and thus all we need to calculate is the Rayleigh scattering.To do this, we need to characterize both the spectrum of the stars and the spectral response of the FRAM system to form an overall spectral response R(λ) (which in general depends on the spectrum of a given star).We considered taking every star to be a black-body radiator of a given temperature, but calculations with realistic stellar spectra have proven this to be a poor approximation (Fig. 8); however, while realistic stellar spectra depend on a number of parameters, the calculated Rayleigh optical depth still correlates well with B − V, at least for stars that are not too red (B − V < 1), as shown in the left panel of Fig. 9.The spectral properties of the system have been measured in field, but to apparently insufficient precision as for example the rather slight change in the extrapolation of the spectral transmission of the lens shown in red in the right panel of Fig. 9 (with results in pink on the left panel) changes the Rayleigh OD prediction at B − V = 0 by 0.01.More precise laboratory measurements are underway.
The Rayleigh scattering as a function of wavelength was then calculated according to [5].For each scan we take the total atmospheric depth from the GDAS model used at the Auger Observatory [6].Here it is important to note that simply taking the convolution of the wavelengthdependent Rayleigh VOD with the overall spectral re-  sponse R(λ) is technically wrong, because exponentiation and integration obviously do not commute.Instead, one has to calculate the reduction in the flux at ground for the light across the system passband as The −0.0012A 2 term is determined numerically and needs to be added to (6) if we want to properly describe stars at A > 8 (Fig. 10) and in typical atmospheric conditions τ FRAM ≈ 0.208 for a star with B − V = 0 with the assumptions on spectral properties of the system and black-body stellar spectra as described above.
The description of stellar spectra can be tested by comparing the color-dependence of the extinction k C fitted from data (0.015) with the predicted value for Rayleigh scattering, because the extinction is mostly Rayleighdominated and the effect of aerosols on k C is small.Using black-body approximation, the predicted value is too small (0.008), using realistic stellar spectra, the value is slightly too large (0.022); further improvement in this area is clearly needed. 1

Discussion
The presented method allows us to extract VAOD values for scans taken without interference of clouds.In general, we obtain very small but positive values starting from almost zero indicating that the method is working correctly.We also observe that the cleaner nights have the tendency to be stable while dirtier nights show large temporal variations.Most importantly, there is no visible correlation between Z and VAOD, indicating that the fitting procedure does indeed remove this degeneracy and that the variations in VAOD are not artifacts related to the variations of the zeropoint.
The overall precision of the measurement is complex to estimate at the moment.The statistical errors in determination of k from the fits are 0.01 and less.The biggest worry here is whether there could not be a systematic effect with the same trend as the extinction that would influence the fitted value, but the consistency of the fits and the small RMS indicate otherwise.The largest uncertainty then comes from the Rayleigh subtraction as explained above, leading to a possible shift of the VAOD values at the order of 0.01-0.02depending on how the laboratory measurements of the optics turn out.Moreover, the current calculation of the Rayleigh scattering now overestimates the k C coefficient: 0.022 is predicted using the spectral properties of the spare lens, while the fitted value is roughly 0.015.Using the theoretical value instead of the fitted one shifts the VAOD by 0.003.

Conclusions and outlook
We have shown that we can use the data on apparent brightness of a large number of stars distributed across a wide range of zenith angles gathered in the Shoot-the-Shower program to extract the integral VAOD values to the top of the atmosphere.The methods of data processing, selection of suitable stars and corrections for various instrumental effects lead to a model that fits the observed difference between catalog and measured values with an RMS of 0.07 magnitudes.The model parameters can be then turned into VAOD values using theoretically calculated values for Rayleigh scattering.There are still important questions regarding the use of realistic stellar spectra and the need for detailed knowledge of spectral response of the FRAM components -however these will be addressed with precise laboratory measurements.
The methods to measure the aerosol content of the atmosphere developed with the Auger FRAM have already found application in atmospheric monitoring for the planned Cherenkov Telescope Array (CTA) observatory.See [7] for more details.
Eventually, we would like to establish FRAM as an alternative source of information about the aerosol content of the atmosphere -even though FRAM can provide only the integral value to the top of the atmosphere, it has the advantage of being able to take many measurements independently of the measurements by the fluoresence telescopes and with no interference to their operations.To this end, we have started (since January 2016) to take altitude scans similar to the Shoot-the-Shower data, independently of any StS triggers, during any night when stars are visible above the FRAM site.We are also regularly taking measurements in various azimuth directions in order to study the spatial uniformity of the aerosols and measurements in different filters (B, V and R) in order to measure the Angstrom coefficient every night.The wavelength dependence of the VAOD is usually assumed to be a power law and the mean value of the Angstrom coefficient α has been measured to be 0.7 at the Auger Observatory [8].However, already these data indicate that the immediate values vary between 0 and 2. Recent data from "nearby" AERONET stations [9] indicate that changes between 0 and 2 can occur during a single day and thus measurement of its value with high temporal resolution is interesting.

Figure 2 .
Figure 2.An example of a photo in a rich star field, not corrected for flat field, shows the significant vigneting in corners that must be removed using a flat field image or by a fit to the stellar data.

Figure 3 .
Figure 3. Aperture photometry: the inner circle is used to calculate the light flux and the annulus defines the background.

Figure 4 .
Figure 4.The effect of the varied star FWHM caused by mechanical problems between camera and lens (right) on the deviation of the measured brightness of those stars from the model (left).

Figure 5 .
Figure 5. Distances between detected sources and catalog positions of stars in pixels (color coded) when the source identification was done after splitting into 3 × 3 tiles (left) and 4 × 4 tiles with inner 2 × 2 merged (right).

Figure 6 .
Figure 6.Residua of the global fit as a function of various variables.

Figure 7 .
Figure 7. Zeropoint as a function of scan number.Consecutive points with the same color belong to the same night, but there might be large gaps between nights, such as the 196-day gap marked in the bottom row caused by problems with the Shoot-the-Shower trigger system.

Figure 8 .Figure 9 .Figure 10 .
Figure 8.Comparison of the overall spectral response R(λ) calculated using black-body radiation (blue) and realistic stellar spectra (green) for stars with surface temperature of 4000 K (left) and 10 000 K (right).