Direct measurements of atomic oxygen in the mesosphere and lower thermosphere using terahertz heterodyne spectroscopy

Atomic oxygen is a main component of the mesosphere and lower thermosphere of the Earth, where it governs photochemistry and energy balance and is a tracer for dynamical motions. However, its concentration is extremely difficult to measure with remote sensing techniques since atomic oxygen has few optically active transitions. Current indirect methods involve photochemical models and the results are not always in agreement, particularly when obtained with different instruments. Here we present direct measurements—independent of photochemical models—of the ground state 3P1 → 3P2 fine-structure transition of atomic oxygen at 4.7448 THz using the German Receiver for Astronomy at Terahertz Frequencies (GREAT) on board the Stratospheric Observatory for Infrared Astronomy (SOFIA). We find that our measurements of the concentration of atomic oxygen agree well with atmospheric models informed by satellite observations. We suggest that this direct observation method may be more accurate than existing indirect methods that rely on photochemical models. Atomic oxygen concentrations in the upper atmosphere can be measured directly with an airborne terahertz heterodyne spectrometer. This approach is probably more accurate than indirect estimates from photochemical models, according to a comparison of the two methods.

A tomic oxygen extends from about 80 km to above 300 km in altitude, but with more than 90% concentrated between 85 and 125 km (Fig. 1a). It plays an important role for the energy balance of the mesosphere and lower thermosphere (MLT), because it participates in exothermic chemical reactions and it contributes to radiative cooling 1,2 . The latter occurs mainly via emission from CO 2 at 15 µm and NO at 5.3 µm. Both molecules are excited by collisions with ground state atomic oxygen. In particular, quenching of CO 2 vibrational levels by collisions with atomic oxygen is important. Therefore, the knowledge about the distribution of atomic oxygen is crucial for retrieval of the kinetic temperature in the MLT from the 15 µm CO 2 radiance 3 . In addition, direct radiative cooling by atomic oxygen occurs via the fine structure transition from the lowest excited state, 3 P 1 , into the ground state, 3 P 2 , at 63.2 µm (4.7448 THz) 4 . This is the dominant cooling mechanism above~250 km. In the MLT the lifetime of atomic oxygen in the 3 P 1 state is several hours. Because of this long lifetime it can be transported over large distances before it releases its energy and therefore it might be used as tracer for the dynamical motions, vertical transport, tides, and winds [5][6][7] . This leads to a strong coupling between the dynamics and the photochemistry in the MLT with the energy being released significantly after and far away from the location, where the UV photon is absorbed. An accurate knowledge of the global distribution of atomic oxygen and its concentration profile as well as diurnal and annual variations are therefore essential for understanding the photochemistry and the energy budget of the MLT.
At the atmospheric conditions of the MLT the fine structure line at 4.7448 THz is thermally Doppler broadened with a line width of ∼12 MHz for emission originating at ∼100 km where the density of atomic oxygen is largest. Due to the increasing temperature, emission from the thermosphere is broader (∼25 MHz at altitudes >300 km). Direct observation of this transition requires airborne or space-borne instruments, because absorption by tropospheric water vapor prohibits observation with groundbased instruments. In addition, terahertz (THz) spectrometers are notoriously complex, in particular when it comes to airborne or space-borne applications. The 4.7-THz line has been observed a few times by rocket-borne instruments 8 . More observations with global coverage have been done with the Cryogenic Infrared Spectrometers and Telescopes for the Atmosphere (CRISTA) that flew on the space shuttle in the 1990s 9,10 . With this spectrometer atomic oxygen densities have been determined at altitudes from 130 to 175 km, which account for about 20% of the total atomic oxygen in the MLT. Below 130 km retrieval was not possible, because of the opacity of the transition and the limb view of CRISTA. Measurements of the fine-structure line have also been made with the Far-InfraRed Spectrometer (FIRS-2), a Fourier transform spectrometer on a high-altitude (∼38 km) balloon 11 . From 1989 to 2003 in total 31 spectra have been obtained and analyzed. But again, the line shape of the transition was not resolved. The radiances observed by FIRS-2 deviate by less than 15% from the computed ones using the Mass Spectrometer Incoherent Scatter (MSIS) model 11,12 .
Atomic oxygen concentrations can be inferred indirectly from measurements of the O( 1 S) green line, the O 2 A-band or the emission from vibrationally excited OH. These methods rely on chemical reaction models and assumptions such as photochemical steady state, quenching rates, radiative lifetimes, and reaction coefficients. The OH emission has been measured with the SABER (Sounding of the Atmosphere using Broadband Emission Radiometry) infrared spectrometer and with the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) 6,13,14 . The O( 1 S) green line was measured by the Wind Imaging Interferometer (WINDII) and SCIAMACHY instruments [15][16][17] 16 . A recent comparison based on an improved photochemical model reveals differences up to 20% between SABER, SCIAMACHY, WINDII, and OSIRIS measurements 19 . The new atomic oxygen concentrations are about 25% smaller than those originally derived 19,20 . Likewise, with an improved photochemical model it was found that the atomic oxygen concentrations derived from OH, O( 1 S), and O 2 A band measurements with SCIAMACHY agree within 15% with each other, which is better than the agreement with SABER data derived from OH measurements 21 .

Results and discussion
Direct observations of atomic oxygen are potentially more accurate. We have analyzed its emission at 4.7 THz measured with the GREAT spectrometer on SOFIA, the Stratospheric Observatory for Infrared Astronomy 22,23 . This is a by-product of the astronomical observations in the same frequency band. In fact, atomic oxygen has been observed in the atmosphere of Mars using GREAT 24 . Due to Doppler shift the line from the astronomical object is often significantly shifted (up to ∼3 GHz) from the atmospheric line. This allows analyzing the spectra of the atomic oxygen line in the Earth atmosphere for those astronomical objects, which are significantly Doppler shifted and not too much broadened, ensuring that the lines do not interfere.
Heterodyne spectroscopy has a number of differences compared to grating spectroscopy and Fourier-transform spectroscopy regarding the measurement of atomic oxygen in the atmosphere. First of all it provides a much higher spectral resolution. In the case of GREAT the spectral resolution of ∼6 MHz (full width at half maximum (FWHM), see "Methods" section) is limited by the line width of the local oscillator (LO) 25 . This is significantly smaller than the atomic oxygen line width. In comparison, CRISTA has a spectral resolution of ∼4.2 GHz 8 and FIRS-2 has 120 MHz spectral resolution 11 . The resolving power of GREAT allows spectrally resolving the atomic oxygen line (Fig. 1b). Because of the high spectral resolution the problem of saturation of the atomic oxygen line can be overcome in the vertical sounding observation geometry of GREAT (see "Methods" section). This occurs in the center of the line while the line wings are not saturated and contain information about the atomic oxygen concentration at higher altitudes. The second advantage is the high sensitivity of the GREAT spectrometer. Its singlesideband noise temperature is 2200 K. This corresponds to a noise equivalent power (NEP) of 3 × 10 −20 W within an integration time of 1 s. With this sensitivity GREAT is capable of measuring one spectrum with an integration time of 19 s. For comparison, 31 separate spectra have been measured with FIRS-2 each requiring between 11 and 16 min of integration time. It should be noted that GREAT is a rather large instrument, which, in its present form, cannot be implemented on a balloon or on a satellite. However, recent improvements in heterodyne technology, in particular high-frequency Schottky diode mixers 26 , compact quantum-cascade lasers as local oscillators 27 , and low-power CMOS (Complementary Metal-Oxide Semiconductor) backend spectrometers 28 have initiated the development of balloon-borne and satellite instruments 29,30 .
The spectra, which are presented here, have been measured on 14 January 2015 between 1:15 and 4:15 am (flight ID 2015/01/14). The flight trajectory is shown in Fig. 1c. Essentially it is in southnorth direction above the Pacific Ocean west of northern Mexico and California between 27°and 48°north latitude and 134°and 128°west longitude. The flight direction is from south to north and the telescope is viewing west. Three typical profiles obtained during this flight are shown in Fig. 2. They have been measured at three different elevations, as well as at different locations along the flight track and different altitudes. The measurement uncertainty in intensity is ±15%. It is caused by uncertainties of the absorption by water vapor and calibration uncertainties (see "Methods" section). At the lowest elevation of 29.1°the profile is the broadest with signs of saturation as the flat top indicates. With increasing elevation it becomes narrower and at the largest elevation of 50.6°saturation almost disappears, because the path through the atmosphere becomes shorter and less atomic oxygen is probed. It should be noted that most of the signal in the wings of the profile is from altitudes around 100 km. The red solid lines are radiative transfer calculations based on atomic oxygen column densities and temperature profiles from NRLMSISE-00 (US Naval Research Laboratory Mass Spectrometer and Incoherent Scatter Radar Exosphere − 2000) 31 . Each line profile is a convolution with the 6 MHz spectral resolution of the GREAT spectrometer. The elevation of the telescope as well as the absorption by the atmosphere between SOFIA and the MLT are also taken into account (see "Methods" section). The residuals (modeled minus measured data) are up to 15% and the largest deviations occur in the wings of the profiles, indicating that the atomic oxygen concentration is somewhat different from the model in particular at higher altitudes.
In order to analyze the accuracy of the model we have calculated the emission profiles assuming various deviations from the atomic oxygen concentration and temperature profiles of NRLMSISE-00 (Fig. 3). First, we have varied the concentration profile from 120 to 80% of the nominal concentration at all altitudes. The residuals show that the best agreement is obtained with the nominal concentrations from NRLMSISE-00. It should be noted that the residuals become more pronounced in the wings of the profile while the peak remains unchanged because of saturation. This underlines the importance of the high spectral resolution. When the temperature, T, is changed by ±5% while the concentration profile is kept nominal the emission becomes much stronger (T = +5%) or much weaker (T = −5%), because more or less atomic oxygen is excited. An important difference to the residuals for different atomic oxygen concentrations is that in this case the residual at the peak of the line changes significantly. In Fig. 3c calculated profiles are shown which are obtained by shifting the concentration profile by +5 km or −5 km relative to the nominal profile while keeping the temperature profile unchanged. Again, the deviation from the measured profile is pronounced. When the concentration profile is shifted by −5 km a single peak occurs and saturation of the emission is reduced. A double-peak emission profile appears when the concentration profile is shifted upwards, because the wings of the emission profile are less attenuated. They probe higher atmospheric layers and more atomic oxygen contributes to the signal.
Because NRLMSISE-00 is not a physics-based model, an agreement between this model and the GREAT measurements might be accidental. Therefore we compare the GREAT data with atomic oxygen profiles obtained with the SABER instrument. There is a set of nine SABER profiles, which was measured close to the geolocation and time of the GREAT measurements (see "Methods" section). This data has been measured along the flight track shown in Fig. 1c  Since there are no measured data available for altitudes above 100 km we have combined SABER and NRLMSISE-00 atomic oxygen and temperature profiles. Below 100 km these consist of the SABER profile while above 100 km the NRLMSISE-00 profiles were taken. With these concentration and temperature profiles the expected atomic oxygen emission is calculated and compared with the measured emission (position B in Fig. 1c). The agreement between the GREAT and the SABER/   NRLMSISE-00 profiles is similar as with the pure NRLMSISE-00 profiles (Fig. 4a). The differences are within the uncertainty of the GREAT measurement. Also, day-to-day variations of atomic oxygen and temperature as well as uncertainties of the SABER data and NRLMSISE-00 may contribute to the differences. An analysis of the altitude profiles is shown in Fig. 4b. Along with the combined SABER/NRLMSISE-00 profiles it displays a range (±1σ) of atomic oxygen concentration profiles, which are compatible with the GREAT spectra. At the peak this is an uncertainty of approx. ±40%. The range of profiles is derived by varying the atomic oxygen concentration profile, the temperature profile and the relative positon of both starting with the best fit NRLMSISE-00 profile. The variations of these parameters were determined, which provide residuals of less than the measurement uncertainty of ±15%. The combined SABER/NRLMSISE-00 profile falls into this range. We now turn to the analysis of the flight path in Fig. 1c with 254 spectra. For each spectrum the radiance has been determined by integration of the emission profile and plotted as a function of flight time (Fig. 5). All data are corrected for absorption by atmospheric water vapor, which is determined by fitting the maximum of the saturated emission profiles (see "Methods" section). In contrast to Fig. 2a single value for the precipitable water vapor (pwv) content has been applied to all spectra (see "Methods" section). The calculated radiances take into account the atomic oxygen concentration and temperature profiles from NRLMSISE-00 and the elevation of the telescope. It should be noted that along the flight path the atomic oxygen concentration as well as the temperature predicted by NRLMSISE-00 change by less than ±5% at all altitudes between 50 and 400 km. If only the region with the highest concentration of atomic oxygen is considered (85-125 km) the changes are less than ±2%. The measured radiances agree with the calculated ones within the measurement uncertainty of ±15%. The radiance increases by about 20%, because the elevation of the telescope decreases leading to a longer path through the atmosphere and as a consequence to broader emission profiles and increased radiances (cf. Fig. 2). The shapes of the measured emission profiles along the flight path of SOFIA agree well with the predictions by NRLMSISE-00 and the corresponding atomic oxygen concentration profiles are well within the uncertainty of the GREAT measurement shown in Fig. 4b. It is worth noting that out of the 31 FIRS-2 radiances 26 are in the range of our measurement between 1.5 and 2.2 nW cm −2 sr −111 .
The first spectrally resolved, direct measurement of atomic oxygen in the MLT through its fine-structure transition at 4.7 THz is presented. This is a very promising alternative to established indirect methods which rely on of photochemical  The measured data are corrected for absorption by atmospheric water vapor using the same pwv value for all data. The straight line has been calculated with the NRLMSISE-00 model taking into account the elevation of the telescope. Data obtained from 254 spectra are plotted. The radiance increases until t = 5000 s and from t = 6000 s onwards, because the elevation of the telescope decreased from ∼51°to ∼32°(until 5000 s) and from ∼42°to ∼29°(from 6000-10,000 s). The interruption between 5000 and 6000 s corresponds to a reorientation of the telescope. b Normalized radiance (black dots, measured radiance divided by calculated radiance) during flight. The measurement uncertainty indicated by the grey bars is ±15%. models. It is an important step towards a conclusive understanding of the photochemistry, energy balance, and dynamics of the Earth atmosphere. This will be substantiated once the large amount of atomic oxygen data, which has been obtained with SOFIA, is completely analyzed. With the recently installed sevenpixel upGREAT heterodyne spectrometer even more data will be available. In particular, the fine-structure transition at 2.06 THz can be measured simultaneously, enabling measurements of wind speed and temperature profile. Beyond that, with the current progress in THz technology balloon-borne and space-borne heterodyne spectrometers become feasible. Combining such a THz spectrometer with optical instruments similar to SABER or SCIAMACHY will be even more advantageous for the determination of atomic oxygen in the MLT.

Methods
The emission of the atomic oxygen fine structure line at 4.7448 THz is observed with GREAT, the German Receiver for Astronomy at Terahertz Frequencies. GREAT is regularly operated on board of the Stratospheric Observatory for Infrared Astronomy (SOFIA) 22 , which is a Boeing 747SP with a 2.5-m diameter telescope in the rear. The elevation of the telescope can be varied from about 20°to 60°and its azimuth is determined by the aircraft heading. At 100 km altitude the area of the sampled volume is about 5 m wide and 5 km long and the distance between two sampled volumes is approximately 15.5 km. This is determined by the spatial resolution of the 2.5-m diameter telescope of SOFIA, which has a diffraction limited beam width of 6 arcsec (FWHM), by the speed of the aircraft, and the integration time for one spectrum, which is 19 s.
GREAT is a heterodyne spectrometer with two frequency channels, one of them at 4.7 THz. The observing altitude of SOFIA is 12-14 km, meaning that the emission from atomic oxygen in the MLT will be present in nearly all of the astronomical atomic oxygen observations. Due to the velocity of astronomical objects relative to the Earth the atomic oxygen emission from these objects is usually Doppler shifted relative to the line in the Earth's atmosphere. Thus the astronomical data nearly always contain the atomic oxygen line originating from the MLT. The GREAT heterodyne spectrometer relies on a hot-electron bolometer as mixer and a quantum-cascade laser (QCL) as local oscillator (LO) 25,32 . It has single-sideband noise temperature of 2200 K enabling a measurement time of 19 s per spectrum and an intermediate bandwidth ranging from 0.2 to 2.5 GHz. The backend spectrometer is a digital fast Fourier transform spectrometer 33 . Its channel spacing is 76.3 kHz, more than sufficient for the present observations, where we, in fact, used channel binning down to effectively 0.763 MHz channels. The spectral resolution is ultimately limited by the emission linewidth of the QCL LO, which is Gaussian-shaped with 6 MHz FWHM. The spectral resolution of the GREAT instrument has been determined prior to the flight in the laboratory be measuring reference spectra of CH 3 OH at a pressure of 1 hPa. From the analysis of these spectra a spectral resolution of 6 MHz (FWHM) was determined. This linewidth is taken into account when comparing the measured radiances with our radiative transfer model.
For radiometric calibration the GREAT spectrometer is equipped with two blackbody calibration sources. One has a temperature T hot of 294 K while the other one is cooled providing a temperature T cold = 149 K. A measurement cycle for the atmospheric atomic oxygen line is as follows: First, the signals C hot and C cold from the hot calibration source and the cold calibration source, respectively, are measured. Then the signal of the sky, C sky , with no astronomical source in the fieldof-view of the telescope is measured. The calibrated signal, S cal , is determined according to Here, B(T hot , ν) and B(T cold , ν) are the radiances of the hot and cold calibration source. For analysis of the spectra we have implemented a radiative transfer code. In this code the atmospheric radiance is modelled by evaluating the following radiative transfer equation with R int being the spectrally integrated radiance (nW cm −2 sr −1 ). T is the temperature at altitude z and ν is the frequency. B ν (T(z)) is the source function at altitude z. As has been shown in previous studies, local thermodynamic equilibrium (LTE) determines the population of atomic oxygen up to at least 350 km 34 . This is in agreement with FIRS-2 measurements, which show that more than 80% of the atomic oxygen radiance originates from below~140 km and its distribution is well modelled assuming LTE 11 . Therefore B ν (T) is the Planck blackbody function. The monochromatic transmission τ ν is given by Here Here ν is the frequency, Q is the electronic partition function and E L is the energy of the lower state of the transition. E L is zero because the 3 P 2 is the ground state.
To calculate the radiance the atmosphere is divided into horizontally homogeneous spherical layers with a thickness d = 1 km starting at an altitude of 400 km and ending at 50 km. Absorption and scattering have been neglected. In each layer the transmission τ ν is calculated assuming an average temperature within the layer, which is taken from the NRLMSISE-00 temperature profile 31 . In each layer a portion of the incoming radiation is absorbed which is taken into account by an attenuation factor exp(−α ν l) (α ν ː absorption coefficient at frequency ν, l: absorption length). The radiation which is thermally emitted by the layer is added. The absorption path length takes into account the elevation φ of the telescope.
Our radiative transfer code evaluates the radiance across the atomic oxygen line with a step size of 0.763 MHz. This ensures that the line profile is properly sampled with at least 20 spectral bins. The monochromatic radiances are integrated in order to obtain the total, integrated radiance (cf. Eq. 2). The integration is done from −35 MHz to +35 MHz around the line center frequency. We have validated our radiative transfer code by calculating the radiances reported in the ref. 11 . As input we used the concentration profiles and temperature profiles provided by NRLMSISE-00. The integrated radiances calculated with our code agree within 1% with the radiances calculated by the radiative transfer code of ref. 11 , which in turn agrees within 1% with the Full Transfer by Ordinary Line-by-Line Methods (FUTBOLIN) radiative transfer code 35 and the Monochromatic Radiative Transfer Algorithm (MRTA) 36 . Figure 6 displays a comparison of the lines shapes simulated with the radiative transfer code assuming a spectrometer with a spectral resolution of 6 MHz (such as GREAT), and an instrument which is not capable to spectrally resolve the atomic oxygen emission (spectral resolution of 120 MHz such as FIRS-2). The spectra are calculated based on atmospheric layers up to a maximum altitude. The highresolution spectra start to saturate in the center of the line when atmospheric layers above an altitude of approximately 110 km are taken into account (Fig. 6a). However, the wings of the emission line become wider when atmospheric layers from higher altitudes are included in the calculation. This demonstrates that atomic oxygen lines observed with a high-resolution instrument such as GREAT contain information from all altitudes. Figure 6b displays the same spectra but normalized in order to better visualize the wings of the profiles, which becomes wider the more atmospheric layers are taken into account. In contrast the line shape observed with a low-resolution instrument is entirely determined by the resolving power of the instrument and independent of the atmospheric layers, which are included in the calculation (Fig. 6c).
In order to compare the calculated radiance with the measured radiance the absorption by water vapor in the tropopause and lower stratosphere above SOFIA has to be taken into account 37 . It depends on the precipitable water vapor (pwv) in the atmosphere above SOFIA. From fitting the line profiles in Fig. 2 we obtain pwv values of 14.2 µm (elevation: 50.6°), 11.1 µm (38.3°), and 10.1 µm (29.1°). This corresponds to zenith pwv values of 11.0 µm, 6.9 µm, and 4.9 µm, respectively. The pwv overburden of SOFIA towards zenith can be calculated using the ATRAN (Atmospheric TRANsmission) model 38 . This yields 9.6 µm (50.6°), 7.3 µm (38.3°), and 5.5 µm (29.1°), which is in good agreement with the fitted values. Figure 7 displays the changes of the emission profile when the pwv is changed by ±1 µm. This corresponds to a change of transmission of ±7%. It should be noted that the transmission affects the whole spectrum by the same factor. The uncertainty of the pwv is taken into account in the total uncertainty (see below).
For obtaining the radiances in Fig. 5 the pwv is deduced from fitting the saturated peak signals of all emission profiles during the flight and taking into account the elevation of the telescope. It yields an average value of 10.9 ± 1.5 µm pwv. This is larger than the deviation shown in Fig. 7 and corresponds to change of transmission by ±12.5%. The uncertainty is estimated from comparing the ATRAN calculation with the fitted pwv data. The single pwv value of 10.9 µm has been used for all radiances in Fig. 5. The deviation of the measured radiances from the calculated ones is larger in the first ∼3000 s of the flight, because the deviation to the pwv of the fit in Fig. 2 is largest. In principle, it is possible to obtain a much better agreement between the measured and calculated radiances by fitting each spectrum with a specific pwv.
The total measurement uncertainty is determined by two parameters, namely the pwv and the reproducibility of the calibration procedure. The 1.5 µm uncertainty of the pwv corresponds to an uncertainty of up to 12.5% of the radiance with the smaller uncertainty at smaller telescope elevation. Another source of uncertainty is the calibration procedure with the hot and cold loads. This includes the temperatures of the loads as well as gain drifts in the system and standing waves, which can change between calibration cycles. This uncertainty is 9%. The total uncertainty of 15% is obtained from the root-sum-square of the pwv uncertainty and the calibration procedure uncertainty.
In Fig. 8a nine SABER atomic oxygen concentration profiles are shown. These spectra have been measured between 1:15 and 4:15 am on 14 January 2015 along the track shown in Fig. 1c. For comparison with the GREAT measurements we have computed the average profile of these nine profiles. This average profile is shown in Figs. 1a and 4b. The nighttime SABER measurements of the whole year 2015 (similar geolocation as GREAT: 20°−40°north latitude, 120°−150°west longitude) are visualized in Fig. 8b. The 1σ variation is less than 30% of the mean at altitudes below 95 km and it increases above.