Laser Remote Magnetometry Using Mesospheric Sodium
Abstract
We have demonstrated a remote magnetometer based on sodium atoms in the Earth's mesosphere, at a 106km distance from our instrument. A 1.33watt laser illuminated the atoms, and the magnetic field was inferred from backscattered light collected by a telescope with a 1.55mdiameter aperture. We theoretically predict a shot noise limited measurement sensitivity of . The measured sensitivity was due to a smaller returned intensity and smaller resonance strength than expected. The value of magnetic field inferred from our measurement is consistent with several models of the Earth's field shape to within a fraction of a percent. Projected improvements in optics, plus the use of advanced lasers or a large telescope, could result in sensitivity.
Key Points
 Remotely measured the geomagnetic field at a 90km altitude, this is an altitude that was previously difficult to monitor
 Measured geomagnetic field was 45,441 nT, comparable to geomagnetic models, with a sensitivity of 162 nT/√Hz
 Several improvements to the present measurements are discussed; the sensitivity could be improved to better than 1 nT/√Hz
Plain Language Summary
The upper end of a spinning top will slowly revolve around a vertical axis, a phenomenon known as precession. The precession frequency of the top is proportional to the strength of the gravitational field in which the top sits. This work uses the precession of atomic spins to analogously measure the strength of a magnetic field. The interaction of light with the atoms is sensitive to this precession and allows existing devices to observe the precessing atoms within an instrument. We have, for the first time, measured the magnetic field 100 km from our instrument using naturally occurring sodium atoms in the Earth's mesosphere. Atomicbased measurements, such as optical magnetometers, are typically performed in wellcontrolled environments. Here we demonstrate the ability to initialize quantum spin states and readout the evolution of the state at a significant standoff distance using a modified laser guidestar. The technique provides unique information about the mesospheric region of the atmosphere that is inaccessible by either aircraft or satellites. The small amount of light collected by our telescope creates a much less sensitive measurement than typical optical magnetometers, but methods are discussed that can improve the performance of this initial setup to reach a sensitivity needed for advancing understanding of geomagnetic field fluctuations such as those due to solar eruptions. This novel measurement technique is also an example of feasibility for other remote quantum measurements in the ambient atmosphere.
1 Introduction/Background
Atoms with an unpaired electron, such as sodium, have ground states with angular momentum and therefore a magnetic moment. Optically pumping a vapor of such atoms creates a net magnetic moment along the direction of the light, which will precess in a magnetic field. The frequency of precession, known as the Larmor frequency, is linearly proportional to the magnitude of the magnetic field. The constant of proportionality is known as the gyromagnetic ratio, and it is known with high precision for sodium and other atoms. Thus, a measurement of the precession frequency of an atom determines the magnitude of the magnetic field in which the atom is immersed. Magnetometry based on this kind of measurement using atoms within the instrument is a wellestablished technology (Budker & Romalis, 2007).
Ninety kilometers above the Earth's surface is a naturally occurring layer of sodium atoms, the residue of meteor ablation. The density of the layer varies with meteor activity but is on the order of 4 × 10^{9} atoms/m^{3} and has a thickness of about 10 km (Happer et al., 1994; McNeil et al., 1995; Plane, 1991). Light at the sodium D_{2}a resonance wavelength of 589.159 nm interacts strongly with these atoms, such that about 4% of the light is resonantly scattered as it passes through the 10kmthick sodium layer. This scattering from mesospheric sodium is what enables the technology of adaptive optics based on “guidestar” lasers. Laser sources with power up to 50 watts, operating exactly at the sodium resonance wavelength, have been built as guidestar lasers for large telescopes incorporating adaptive optics (Denman et al., 2005).
In 2011, Higbie et al. (2011) proposed a system for remotely measuring the magnitude of the geomagnetic field, typically called F by geophysicists, by interrogating sodium atoms in the mesosphere using a pulsed laser at the sodium resonance wavelength. In this measurement, the pulse repetition frequency (PRF) of the laser is adjusted to match the Larmor frequency of the atoms in the magnetic field. If the intensity of the light is in the appropriate range, then the amount of light backscattered by the sodium atoms will be enhanced when the PRF is equal to the Larmor frequency. By tracking the frequency of this resonant enhancement, the magnetic field can be measured remotely.
The existence of a resonance is due to two phenomena. First, optically pumping with circularly polarized photons creates atoms with a net magnetic moment, which precess at the Larmor frequency. Second, there is an enhanced probability for backscatter—that is, scatter back toward the source—from the polarized atoms. Once precessing, they are optimally oriented for enhanced backscatter once per Larmor cycle time. Working at the guidestarlaserequipped Kuiper Telescope (University of Arizona, 2018), one of the facilities of the University of Arizona's Steward Observatory, we have observed this resonance and have thus measured the magnetic field in the sodium layer. The sensitivity of the magnetic field measurement as characterized by an equivalent noise spectral density was . A measurement of magnetic field made by averaging for a period of 1 hr would have a noise bandwidth of (3,600 s)^{−1} = 0.28 mHz. With this bandwidth, the rootmeansquare (RMS) uncertainty due to measurement noise would be 2.7 nT.
The remainder of this paper is divided into four subsections, and then we summarize the work in our conclusion. Section 2 derives equations for the sensitivity of the instrument, based on the parameters of the return signal, and assuming that shot noise in the return signal is the predominant noise source. Section 3 describes the expected nature of the signal returned from the sodium atoms. Section 4 describes our experimental setup and results and makes comparison with the expected results. Section 13 discusses the results and ways how our system could be utilized or enhanced.
Figure 1 is a photograph of the Kuiper Observatory showing the guidestar system in operation, taken from atop Mount Lemmon, at a 6km distance. Figure 2 shows, on the left, the Kuiper Telescope dome on a moonlit night during operation of the guidestar system and, on the right, the 61in.diameter (1.55 m) aperture telescope. Figure 3 shows the Frequency Addition Source of coherent Optical Radiation (FASOR) laser source, which sits on a 4ft by 8ft (1.2 m by 2.4 m) optical table.
2 Description, Results, and Discussion
2.1 Instrument Sensitivity
Figure 4 shows a generic function, R(f), which can be used to describe a resonance peak that sits atop a frequencyindependent background. For simplicity, we use a triangleshaped resonance function; the use of another resonance function would change the sensitivity as expressed by equation 8 by a constant of order 1 but would not change the scaling of that expression. The resonance has a full width at half maximum of Δf, in units of hertz. The vertical axis is dimensionless. The resonance sits atop a background level of one. The height of the resonance above the background is H, most conveniently thought of in units of percent. The peak of the resonance shape function is at f = 0 Hz.
Measurement noise will create changes in S that limit the ability to measure small changes in F. In the following, we find the magnitude of the change in F that is needed to make a change in S equal to the change in S that is due to noise.
The first equivalence is the result of the fact that by definition a change in signal ΔS creates an apparent change in magnetic field of CΔS, and the second equivalence is the result of substituting equations 4 and 5 into 6. A simplifying assumption used here is that H ≪ 1, so that S ≈ S_{0}, and the two are considered equal for purposes of noise calculation.
If the resonance function R(f) is a Lorentzian, and if f is selected to be at the steepest point on the resonance (that is, the inflection point), then a factor of multiplies the right side of the equation. If it is not the case that H ≪ 1, the analysis is somewhat more complicated, since the shot noise would then depend on where you are on the resonance. But for analyzing the magnetometer described in this paper, equation 8 is adequate.
A small value of NES is desirable. Thus, we want a large value of signal S_{0}, which can be accomplished with a more powerful laser or a larger collection telescope, and also the smallest possible value of Δf, the width of resonance peak. We also want a large value for the dimensionless resonance strength, H, so we have a large change in signal with a small change in F (i.e., a large slope in equation 1 yields a sensitive detection of magnetic field changes). Though these parameters are not as readily changed, they do depend on the intensity of the laser at the sodium layer and on other laser properties. The gyromagnetic ratio g is of course a constant.
Equation 8 applies regardless of whether the signal S_{0} is entirely from sodium signal return, or if it is due to some other source, as long as shot noise is the dominant noise. If there is light from a signal other than sodium, that light will not show a resonance, and H will be reduced, more than offsetting the benefit of the increase in S_{0} in equation 8.
2.2 Expected Signal

 Ψ

 is the sodium backscatter cross section, in units of photons (s·sr·atom (W/m^{2}))^{−1};

 A

 the collection area of the telescope, in units of m^{2};

 P

 the timeaveraged power transmitted to the sky, in units of watts;

 N

 the sodium atom density per unit area looking upward, in units of atoms/m^{2};

 η_{atmosphere}

 the oneway efficiency of atmospheric transmission;

 η_{optics}

 the efficiency of the telescope optical train;

 η_{detctor}

 the quantum efficiency of the detector, in units of counts per incident photon;

 θ

 the elevation angle of the telescope, with 90° being straight up, and;

 z

 the vertical distance to the sodium layer in units of meters.
It should be noted that the backscatter cross section can also be expressed in units of area, according to the equation σ = Ψhυ, where hυ is the photon energy, 3.4 × 10^{−19} J. This yields a value of σ = 8.2 × 10^{−17}m^{2}, an extraordinarily large value, indicating the extraordinary ability of sodium to backscatter light at its resonance wavelength. Equation 9 would then be changed by replacing P with P/hυ. Equation 9 is essentially the light detection and ranging (LIDAR) equation (Wandinger, 2005). When the laser is at an angle other than 90°, the distance to and the thickness of the sodium layer are both increased by 1/ sin (θ). The value z in equation 9 is the vertical component of the distance to the scatterers, not the total distance.
Table 1 gives our estimated values for each parameter of equation 9, and a resultant value for S_{0}, for the experimentally realized case that will be described in section 4.
Parameter  Value 

Ψ  240 photons (s·sr·atom (W/m^{2}))^{−1} 
A  1.8 m^{2} 
P  3.8 W 
N  4 × 10^{13} atoms/m^{2} 
η_{atmosphere}  84% 
η_{optics}  70% 
η_{detector}  0.27 counts/photon 
θ  60° 
z  92 km 
S_{0}, using equation 9  896 ms^{−1} 
S_{0}, 35% duty cycle  314 ms^{−1} 
The first value of S_{0} corresponds to expected return when the laser is operated continuously. The second is what is expected when it is pulsed with a 35% duty cycle, with the time between pulses much shorter than the time spent passing through the sodium layer. In this latter case, the average power in the sodium layer is 1.33 watts, which is 35% of the continuous wave (cw, i.e., not pulsed) output of the laser. The return signal is continuous since many pulses are in the sodium layer at once.
With S_{0} estimated, we still need an estimate of resonance height H and resonance width Δf in order to complete the right side of equation 8 and have an estimate of the noiseequivalent signal in units of .
 Laser wavelength. We assumed a laser centered at the D_{2}a line of sodium, at 589.159 nm.
 Laser linewidth and sidebands. We assumed essentially zero linewidth, with no sidebands. These assumptions are well met for the FASOR we have used.
 Laser intensity at the sodium atoms. We varied this parameter over our range of interest.
 Laser modulation. We assumed that the laser was pulsed, with square pulses at a constant PRF and duty cycle. From model run to run, we varied PRF and duty cycle. Though our modeling suggested a PRF closer to 20% duty cycle might be best (Denman et al., 2012), we found from our onsky measurements that a duty cycle of 35% (pulse length equal to 35% of the PRF cycle) was close to optimal, and for simplicity all results shown in this paper are for a 35% duty cycle.
 Laser polarization. All results are for circular polarization, which is optimal.
 Magnetic field. We assumed 45.3 μT, from National Oceanic and Atmospheric Administrationmaintained National Geophysical Data Center software “Geomag” version 7.0 using the International Association of Geomagnetism and Aeronomy geomagnetic model coefficients predicted for 92km altitude over Tucson, Arizona, United States (Thébault et al., 2015).
 Magnetic field orientation. The optimal orientation for creating a strong Larmor resonance is to propagate the laser, such that at the 92km altitude of the sodium layer, the direction of propagation of the laser is perpendicular to the Earth's magnetic field. Ideally, this direction would have been 9.75° azimuth and 31.03° elevation. However, due to the telescope's English Equatorial mount constraints, see Figure 2 (right), and the restraints placed on us by the Federal Aviation Administration, we were limited experimentally to 25.1° azimuth and 59.7° elevation. Thus, the angle between the geomagnetic field and the laser propagation was close to 60° and that is what we modeled. Our earlier report (Denman et al., 2012) describes how the magnetometer sensitivity becomes worse as this angle decreases from 90°; thus, our sensitivity at this angle is about 2 times worse than if we had propagated optimally. All theoretical results in this paper use this angle of 60° to correspond to the limitations and constraints at the Kuiper Telescope. The optimal 90° angle may likewise only be possible in the tropics; 60° or better should be possible over most midlatitude locations.
 Atom environment. We used mean times between interatomic collisions and spin relaxation appropriate for the high mesosphere.
Holzlöhner et al. (2010) used almost the same model of sodium atomic physics for their paper on expected guidestar return, and that paper is a good reference if further understanding of sodium modeling is desired.
Figures 5a and 5b show results of model calculations of the resonance height H_{Na}, in units of percent, and the sodium resonance width Δf, in units of kilohertz, as a function of the timeaveraged intensity of light at the sodium atoms. All atoms are assumed to see the same intensity, as would be the case with an ideal “tophat” beam shape. (Here it is assumed that all the signal S_{0} is due to sodium, and for this case there is no distinction between H and H_{Na}. Later, we will define the difference.) The atmospheric pressure is an important variable. At lower pressure, the mean time between collisions of sodium atoms with air molecules is longer, and the light has more time to act on the sodium atoms before their angular momentum states are randomized by collisions, so the same degree of synchronized angular momentum can be achieved at lower light intensity. Over the ~10km thickness of the sodium layer, atmospheric pressure varies by a factor of 6. Figures 5a and 5b each show three cases: for atoms at the center of the sodium layer at 92 km, and for atoms above 25% and 75% of the sodium layer, at 89 km and 95 km. The atmospheric pressure at these three altitudes is 0.14, 0.23, and 0.084 Pa. The parameters in the model corresponding to atmospheric pressure are the rates of velocity changing collisions and spin relaxation, which were assumed to be linear functions of pressure.
It can be seen from inspecting Figures 5a and 5b, in light of equation 8, that there is an optimum intensity for sensitive magnetometry. To minimize the sensitivity, as is desired, the ratio Δf/H must be minimized. Figure 5c was generated using equation 8, with g = 6.99812 Hz/nT, taking S_{0} for 35% duty cycle from Table 1 and taking H and Δf from Figures 5a and 5b. At the altitude of 92 km, we calculate a sensitivity of at approximately 0.5 W/m^{2}. It will be seen that experimentally observed values are almost an order of magnitude worse, primarily because values of both S_{0} and H turned out lower than expected.
2.3 Experimental System
2.3.1 Transmitter
We call our light source a FASOR. This source generates light at the wavelength of 589.159 nm by frequency summing the output of two singlefrequency lasers, one at 1,319 nm and the other at 1,064 nm. (All available coherent sources at 589 nm are based on nonlinear conversion from infrared lasers, since no practical laser at 589 nm is currently available. Dye lasers at 589 nm have been built but have many problems.) The output of the FASOR is up to 10 watts of diffraction limited, singlefrequency coherent light, easily tunable to exactly the sodium D_{2}a resonance wavelength of 589.159 nm. This FASOR was designed and built originally for use as a concept prototype guidestar laser (Bienfang et al., 2003; Denman et al., 2004). It now belongs to the Steward Observatory and is installed in the Kuiper Telescope dome, atop Mount Bigelow, in the Catalina Mountains north of Tucson, Arizona. Figure 6 is a schematic diagram of the measurement system.
This light is passed through an acoustooptic modulator, model 350853 from Gooch & Housego. This device can deflect up to 80% of the light, with a rise and fall time of 300 ns. The deflected, firstorder output beam from the modulator was utilized so that the offstate power was 0 watts. After modulation, the polarization state of the light could be adjusted into any arbitrary state by a quarter waveplate and a half waveplate, in combination. A beam expanding telescope expands the slightly elliptical beam so that it has 1/e^{2} radii along the two transverse axes of 38 × 26 mm. Not shown in the diagram is a flat mirror that directs the beam to the sky. We used a shear plate to optimize the telescope for collimation of the light. We measured both the power and the polarization state of the light after the final flat mirror. The polarization was adjusted to circular. The total power exiting to the sky was 3.8 watts, when the acoustooptic modulator was continuously deflecting the beam. When modulated with a 35% duty cycle, average power was 1.33 watts.
2.3.2 Receiver
We used the 61in. (1.55 m) Kuiper Telescope as the receive telescope to observe the beam in the sodium layer. This telescope has an effective aperture of 1.8 m^{2}, taking into account the light blocked by the 0.41mdiameter secondary mirror.
The photodetector was a Hamamatsu Multipixel Photon Counter (MPPC), model 12662150. (Other vendors call this device a silicon photomultiplier or an avalanche photodiode array.) This module contains a semiconductor detector in the form of a 1mm square. The quantum efficiency is specified as 27%. The dark signal is specified to be below 40 femtowatts, equivalent to a dark count rate of 32 ms^{−1} at 589 nm. The internal avalanche gain of this device is large enough that the detection system is limited by the shot noise in the dark signal, not by noise in the electronics.
A narrowband optical filter, a Chroma model ET592/21x, was used to eliminate light that was not near the sodium wavelength. It had a pass band of 21 nm centered at 592 nm and 97% transmission at 589 nm. Without this filter, signal due to background light was significant; with the filter in place, we could operate under a full moon with no significant signal from background light.
2.3.3 Signal Acquisition System
In an idealized system, free of any noise other than shot noise in the detected signal, we could map out the Larmor resonance by slowly scanning the PRF of the transmitted light. Once mapped out, we could then tune the PRF to the point of steepest slope on the side of the resonance. Changes in signal detected would then correspond to changes in magnetic field F, as per equation 4.
In reality, changes in return signal due to scintillation induced by wind turbulence (the atmospheric effect that causes stars to twinkle) (Dravins et al., 1997; Ryan et al., 1997), laser power fluctuations, beam pointing fluctuations, and other lowfrequency disturbances lead to changes in return signal significant compared to the resonance height H. The technique of phasesensitive or lockin amplifier detection allows the measurement to be moved to a higher frequency, where shot noise is the dominant noise. The lockin amplifier we used was model 7230 from AMETEK Signal Recovery.
To implement phasesensitive detection, the PRF of the system is dithered. This allows the PRF to be moved quickly on and off the resonance or to be moved quickly from one side of the resonance to the other. The lockin detector looks for a signal at the dither frequency, but with a phase that takes into account the delay due to the ~200km round trip of the light.
In practice, a dither frequency of a few hundred hertz was adequate to achieve shot noiselimited detection (Denman et al., 2012). We used a dither frequency of 314.159 Hz. The frequency synthesizer used was the Rigol model DG1062Z, which had all of the features needed to create the frequencydithered pulse train with a few parts per billion frequency stability, while also providing a reference signal at the dither frequency for the lockin detector. Thus, both the dither oscillator and the pulse generator of Figure 6 are included in the Rigol device. The accuracy of the magnetic field measurement will be directly limited by the accuracy of the frequency synthesizer.
2.3.4 Characterizing the Beam in the Mesosphere
Figure 7 shows the beam in the mesosphere, as recorded by a chargecoupled device (CCD) camera in the same location as that used by the Hamamatsu MPPC. The wedge of light coming in from the top is light scattered from air due to Rayleigh scattering. The spot centered in the box is the sodium signal. The separation between the sodium return and the Rayleigh return is due to the fact that the axis of the transmit telescope is offset by about 2 m from the axis of the receive telescope, so that the highdensity region of the atmosphere where Rayleigh scattering is significant is seen offset from the sodium spot. The transmit axis is launched above the receive telescope, as is seen in the image in Figure 7.
The box in Figure 7 corresponds to the size of our 1mmsquare detector at an auxiliary focal plane (BigGuider path) of the telescope. The plate scale at the focal plane is 20 arcseconds/mm (University of Arizona, 2018) or 97 μrad/mm. At a range distance of 106 km, the side of the box corresponds to 106 km × 97 μrad = 10.3 m. A fit of the beam in the box gives a 1/e^{2} radius of 2.6 m. Since atmospheric effects will increase the apparent size above the actual size, this is a good upper limit on our estimate of the beam size in the sodium layer.
A lower limit on beam size can be calculated by assuming that the propagation from the transmit telescope to the mesosphere at 106 km is diffraction limited, so that propagation is according to standard Gaussian beam equations. Starting at 38 × 26 mm, this leads to a spot at the 106km range with 1/e^{2} radius of 0.52 × 0.76 m.
These limits on beam size allow us to calculate beam intensity in the mesosphere. Assuming 3.8 watts launched, and 84% atmospheric transmission, the intensity on the beam axis, given by 2Pη_{atmosphere}/(πω_{x}ω_{y}) where ω_{x} and ω_{y} are the beam radii in the two dimensions and other variables are as in equation 9, must be in the range from 0.3 to 5.1 watt/m^{2}. Note that this is the spatial and temporal peak intensity—not the timeaveraged, or spatially averaged, intensity. Our duty cycle when collecting signal for magnetometry is 35%, so timeaveraged intensity, on beam axis, is 35% of what we have just stated—that is, between 0.1 and 1.8 W/m^{2}. This can be compared to the range from 0.35 W/m^{2} to 1.05 W/m^{2} where optimum sensitivity is expected, according to Figure 5.
2.3.5 Identifying the Sodium Return Signal
On 25 March 2016 we gathered return signal data. The first data set, shown in Figure 8, was done at the low PRF of 700 Hz, with pulses of duration 0.29 ms. At this low PRF, the signal from the sodium layer can be separated from closer sources by time of flight. (At the higher PRF we used for magnetometry, this is not possible, since the distance between adjacent pulses is short compared to the thickness of the sodium layer, so that multiple pulses are in the sodium layer simultaneously, and the return signal is continuous.) Figure 8 shows return data as a function of time after the beginning of pulse transmission, averaged over multiple cycles to improve signal to noise.
The signal before 0.6 ms is not due to sodium but, instead, is Rayleigh atmospheric scattering. Given that the image of Figure 7 shows good separation between sodium and Rayleigh, that is, little Rayleigh in the drawn box, we were surprised by the significant amount of nonsodium light detected. The ratio of sodium signal to total signal (ignoring dark signal), as found by integrating the data of Figure 8, is 53%.
The sodium signal centered at 0.85 ms has a level, after dark count subtraction, of 250 ms^{−1}. This is disappointing compared with the predicted value from Table 1 of 896 ms^{−1}. The low level of signal may be due to a lower than expected value of sodium atom density N, which can vary by a factor of 2. It is also likely that there are unexpected losses in the telescope optics that reduce η_{optics} below the 70% value of Table 1. The signal after 1.1 ms is dark signal. The dark count is estimated from the data of Figure 8 to be 17 ms^{−1}, comfortably below the 32 ms^{−1} derived from the detector specification.
2.3.6 Measurement of Larmor Resonance
The data set of Figure 9 shows the magnetic resonance. For this data set, the PRF was dithered by ±15 kHz from its average frequency f_{avg}. Thus, the PRF is never at this average; it is either 15 kHz above or 15 kHz below. Over a period of 90 min, we slowly swept f_{avg} from 283 kHz to 372 kHz. The dither frequency—the frequency at which the PRF makes a full cycle from f_{avg} + 15 kHz to f_{avg} − 15 kHz and back—was 314.159 Hz.
The signals plotted in Figure 9 are the modulation seen on the return, at the dither frequency, in phase and 90° out of phase (inquadrature) from the delayed dither frequency. We adjusted the lockin phase delay to compensate for the time of flight, with the goal of putting all signal into the inphase channel. The phase correction used was 82.7°, which at 314.159 Hz is equal to (82.7°/360°)/314.159 Hz = 0.731 ms or a roundtrip distance of 219 km. For our telescope angle of 59.7°, and site elevation of 2.5 km, this means phase was optimized for atoms at an altitude of 219 km × sin (59.7°)/2 + 2.5 km = 97.1 km.
When the return signal is larger during the first half of the dither cycle, corresponding to the time when light pulsed at f_{avg} + 15 kHz is returning to the telescope, the inphase value from the lockin is positive; conversely, when pulses with PRF f_{avg} − 15 kHz create an increase in return, the inphase signal is negative. Thus, the positive signal peaks 15 kHz below the Larmor frequency, and its negative image appears 15 kHz above the Larmor frequency. Thus, the data of Figure 9 show the underlying resonance but convolved with a function that creates one peak shifted −15 kHz and another peak inverted and shifted +15 kHz.
If lockin phase were adjusted ideally, the inquadrature signal should be pure noise. However, since the signal has a range of delay times, there is no single phase value that compensates for all ranges; both the nearer and farther sodium atoms contribute some signal to the inquadrature channel. The average count rate S_{0} during the acquisition of the data of Figure 9, due to all sources (sodium return, Rayleigh return, and dark signal), was 151 ms^{−1}.
The inphase data of Figure 9 can be deconvolved and normalized to yield the resonance height H_{Na} and the resonance width Δf. Figure 10 shows the result of deconvolving the data of Figure 9 and normalizing by the amount of DC sodium return signal observed far from resonance. (By DC signal, we mean signal observed directly, rather than using lockin detection that only can observe an AC signal.) The signal due to Rayleigh and to dark counts were subtracted off, so the data of Figure 10 contain only sodium signal and thus can be compared with the predictions of Figures 5a and 5b.
Figure 10 shows that the resonance height, which we have called H_{Na} to remind us that the reference level is the offresonance signal due to sodium, is about 5.9%. The steepest part of the resonance, which is where you want to operate to be highly sensitive to a magnetic field change, can fit to a triangle function with a full width at half maximum of 4 kHz. (A triangle is a rather simple function that clearly does not fit the data over a wide range of frequencies. But looking back at equation 4, observe that all that really matters for magnetometry is the maximum slope of the resonance. The use of a triangle is a simple way to estimate the peak slope of the data and to get an effective value of Δf that will result in that peak slope.)
The value of H_{Na} is interesting for comparison to the model results of Figure 5b, but it is not the value of H needed for use in equation 8. For that, we need the change relative to all signal, including that from dark signal and from Rayleigh signal. This yields a value of H = 2.8%. With values of S_{0} = 151 ms^{−1}, H = 2.8% and Δf = 4 kHz, and the known value of g, we can use equation 8 to calculate an expected system sensitivity of .
Figure 11 shows how we determined both S_{0}, the average count rate of the detector, and V_{DC}, the DC voltage needed for the normalization that led to Figure 10. Figure 11 is the direct output of the detector, recorded for a period of 0.2 ms, under conditions identical to those of Figure 9. The individual pulses created by detected photons can be seen.
2.3.7 Magnetic Field Measurement
The Larmor resonance peak frequency f_{L}, and thus the magnetic field magnitude F, can be estimated from the data of Figure 10. Creating Figure 10 required a long, slow scan of frequency, with most of the frequencies during the scan being far off resonance, where there is little magnetic signal. If one's goal is to quickly measure small changes in magnetic field, a better approach to the estimation of f_{L} is to reduce the size of the frequency dither, so that the full amount of frequency dither is Δf. (That is, a dither of ±Δf/2 from the center frequency f_{avg}.) If f_{avg} is set at f_{L}, then the lockin detector will produce a signal of zero, since the signals at f_{L} + Δf and at f_{L} − Δf are the same. (We are assuming a symmetric resonance function. Asymmetry will introduce a slight offset.) Any change in f_{L} − f_{avg} due to change in magnetic field or due to change in f_{avg} will produce the largest possible change in signal, since both f_{L} + Δf and f_{L} − Δf are set at the inflection points, where the slope is maximum. This change in signal is exactly that expressed by equation 4.
Figure 12 shows data taken in this manner. The lockin signal was recorded, for 45 min, from 10:11 to 10:56 UTC on 25 March 2016. The value of f_{avg} was set at 318 kHz for most of this data run. The dither value was set to ±2.5 kHz for all of the run, with a dither frequency of 314.159 Hz. From 10:31 until 10:33, f_{avg} was shifted to 316 kHz. From 10:54 until the end of the run at 10:56 f_{avg} was set to 320 kHz.
The two, 2min long, 2kHz frequency changes created exactly the signal that would be expected if a magnetic change occurred, with an absolute value of 2 kHz/g = 286 nT. Thus, Figure 12 provides a direct way to measure the sensitivity of the magnetometer.
The two 2kHz signal injections, equivalent to 286nT changes in F, created a signal, estimated from Figure 12, of 1.75 μV, with the sign of the voltage change the opposite of the sign of the frequency change, and the same as the sign of the spoofed magnetic field change. Thus, the calibration factor C is 286 nT/1.75 μV = 163 nT/μV. Ignoring the 4 min of injected signal, the RMS fluctuation of the data of Figure 12 is 0.13 μV, corresponding to 21 nT. The equivalent noise bandwidth of the signal of Figure 12 is 0.0167 Hz, as determined by the lockin amplifier output filter time constant and slope (a 10s time constant with a 12dB/octave slope) and reported by the equivalent noise bandwidth (or ENBW) command. These settings were applied prior to collection of the signal data and are used by the lockin amplifier to apply to the envelope of the frequency response function of its digital finite impulse response output filters. Thus, the amplitude spectral density of the background signal of Figure 12 is .
We have no reason to believe that any actual fluctuation of magnetic field caused any of the signal in Figure 12. The 25 March 2016 was magnetically quiet. The field magnitude F as measured by the United States Geological Survey Tucson Magnetic Observatory (see Acknowledgments) stayed within a range of width 1.02 nT over the 45 min of Figure 12, small compared to the 21nT RMS noise seen in that figure. Thus, it is likely that all of the background noise is due to the measurement process, and thus that value is a fair characterization of the sensitivity of our instrument. However, it is possible that field fluctuations at 92 km exceed those at ground level, and that some real signal contributes to the noise.
Figure 13 shows the amplitude spectral density of the background noise of the Figure 12 data, that is, the data outside the two 2min injected signals. It shows the high variability of nonaveraged spectral data, but nothing makes us believe that it is not simply lowpassed white noise.
The magnetic field F corresponding to a Larmor frequency f_{L} = 318 kHz according to equation 3 is 45,441 nT. Ignoring the offset spikes, the average signal during the 45 min was −0.005 μV, equivalent to a magnetic field change of −0.8 nT, which is negligible.
2.3.8 Expected Magnetic Field
The United States Geological Survey operates a magnetic observatory conveniently located for our work. The Tucson Magnetic Observatory is at location 32.175°N, 110.734°W, at an elevation of 946 m. Our magnetometer's laser system transmitter and the Kuiper Telescope receiver were at 32.416°N, 110.735°W, at an elevation of 2,510 m. Our beam was launched with an elevation of 59.7° (90°equals vertical) and an azimuth of 25.1° (east of due north). Table 2 gives the locations where our beam intersected the sodium layer at three altitudes. Also included is the measured or calculated magnetic field from three different models at each location. The frequency value in the last column is calculated using equation 3, using the value of F from the CHAOS6 model (Finlay et al., 2016).
Altitude(km)  Latitude(deg)  Longitude(deg)  F(nT) IGRFa  F(nT) WMMb  F(nT)CHAOS6c  F(nT) measured  f_{L} (kHz)  

Remotely measured value  45,441  318.0  
Beam crosses bottom of sodium layer  87  32.820  −110.511  45,560  45,551  45,580  319.0  
Beam crosses center of sodium layer  92  32.844  −110.498  45,464  45,456  45,484  318.3  
Beam crosses top of sodium layer  97  32.868  −110.485  45,369  45,361  45,389  317.6  
USGS Tucson Magnetic Observatory  0.95  32.175  −110.734  47,109  47,098  47,121  47,404  329.8 
The Tucson Magnetic Observatory measured the value of F = 47, 404 nT during our 45min measurement. The estimated values are from the International Geomagnetic Reference Field, (Thébault et al., 2015), the World Magnetic Model (Chulliat et al., 2015; National Oceanic and Atmospheric Administration (NOAA), 2015), and the CHAOS6 model (Finlay et al., 2016). The F value derived from our measurement, 45,441 nT, is consistent with the estimates. Our measured value of f_{L} = 318 kHz corresponds to an altitude within but somewhat above the center of the sodium layer. This makes sense in light of the fact that the low pressure of high altitude increases the sodium resonance H, as is seen by observing Figure 5a.
A model prediction of the magnetic field in the mesosphere is difficult because there are no current observations that can be integrated into a model. Hopefully, our described technique and future observations using it can help calibrate these models. All current models use some combination of observations from magnetic observatories at the Earth's surface and satellite observations from the Swarm Constellation (a European Space Agency mission to study the Earth's magnetic field), currently in polar orbits at 450 and 530 km. This absence of mesospheric observations also makes the calculation of the model uncertainty difficult. The time period of our observation was a very quiet one geomagnetically, with a magnetic disturbance storm index Dst of −2 to −3 nT and a solar flux activity index F10.7 value of 85.1 solar flux unit (1 solar flux unit = 10^{−22} W·m^{−2}·Hz^{−1}; National Oceanic and Atmospheric Administration, 2016; WDC Kyoto Observatory, 2016). This corresponds to expected external fields, and corresponding induced fields, at mesospheric altitudes of about 10 nT (Sabaka et al., 2004). The unmodeled field associated with shorter wavelengths of the crustal field is of order 50 nT in the mesosphere (Lesur et al., 2015). The World Magnetic Model (Chulliat et al., 2015) provides a onesigma estimate of 152 nT for commission and omission errors.
2.4 Discussion of Results
2.4.1 Sensitivity
In previous sections, we have provided three estimates of instrument sensitivity, as measured in .
The first estimate is based on calculated system parameters. This estimate of used a value of photon return S_{0} based on values from Table 1 and used values of H and Δf based on the sodium modeling summarized by Figure 5, with the assumption of optimized beam intensity and with all signal being sodium signal. Equation 8 then provides the estimate of sensitivity.
The second estimate, , is based on observed system parameters and also uses equation 8. It is based on measured values of S_{0} and values of H and Δf taken from the fit to the data of Figure 10.
The third estimate, , is based on actual measurement of both signal and noise, based on the data presented in Figure 12.
Case  Total count rate, S_{0}  Na count rate, S_{0, Na}  H_{Na}  H  Δf  Noiseequivalent signal, using equation 8 

S_{0}, H, and Δf from models  314 ms^{−1}  314 ms^{−1}  20%  20%  11 kHz  
S_{0}, H, and Δf as measured  151 ms^{−1}  71 ms^{−1}  5.9%  2.8%  4 kHz  
Effective S_{0}, H, and Δf to fit actual  102 ms^{−1}  48 ms^{−1}  5.9%  2.8%  7.3 kHz 
where S_{0} is the total received signal count rate, including Rayleigh scatter and dark counts, and S_{0, Na} is the count rate from sodium scattering alone.
Our modelbased calculation is based on the count rate S_{0} of 314 ms^{−1} taken from Table 1, with all of S_{0} assumed to be due to sodium scattering. In fact, as measured, only 71 ms^{−1} was produced by sodium; 63 ms^{−1} was due to Rayleigh, and 17 ms^{−1} was dark counts. Since the Rayleigh and dark counts contribute to the background without increasing the resonance, they reduce the useful value of H to 2.8%, according to equation 10.
 At 71 ms^{−1}, the sodium signal was 22% of what we expected, which could be due to natural variations in the density of the sodium layer.
 We did not expect a significant Rayleigh signal to overlap and interfere with our sodium signal; in fact, it was almost as big as the sodium signal, thus reducing by more than half the effective value of H.
 The intensity was probably as low as 0.1 watt/m^{2}, outside the optimal range of 0.35 watt/m^{2} to 1 watt/m^{2}. This reduces the value of both H and Δf, especially at lower altitudes.
We must also discuss the discrepancy between the sensitivity estimated based on observed system parameters and the sensitivity derived directly from the data of Figure 12, the two differing by a factor of 2.2. It could be that the noise is higher than expected, or the signal could be lower than expected. (It could also be that we have measured some real geomagnetic noise, but we will assume that this is not the case.)
One way to characterize excess noise in a detector is as a reduced quantum efficiency, leading to a reduced effective count rate for purposes of calculating noise. Looking to equation 8, observe that a reduced value of S_{0} will increase (that is, degrade) the sensitivity in inverse proportion to S_{0}^{1/2}. Thus, the effective value of S_{0} presented in the bottom row of Table 3 is reduced from the value in the row above by a factor of (1/1.21) 2 = 68%. This factor can be thought of as multiplying the quantum efficiency of the detector to create a reduced, effective quantum efficiency of 27% × 68% = 18%.
One likely source of this excess noise is the considerable variance in the size of the pulses created by photons in the semiconductor detector, due to the nature of the device (Adamo et al., 2014; Hamamatsu, 2013; Hamamatsu, 2016). This can be observed in Figure 11. One cause of pulse height variation is double pulsing, caused by the fact that the electron avalanche creates a small amount of light, which can create a second, simultaneous avalanche.
We can calculate an expected value of ΔV based on the parameters extracted from Figures 10 and 11. With V_{DC} = 228 μV, H = 2.8%, Δf_{injected} = 2 kHz, and Δf = 4 kHz, equation 11 yields ΔV = 3.2 μV.
When compared to the observed value of 1.75 μV taken from Figure 12, the actual signal is lower than expected. This reduction in signal raises (degrades) the sensitivity by a factor of 3.2/1.75 = 1.83. Multiplying the increased noise factor of 1.21 by the reduced signal factor of 1.83, we get 2.2, closely matching the observed difference between calculated and observed sensitivity.
The reduced signal may be due to the fact that the value of Δf = 4 kHz, taken from the triangle fit in Figure 10, is optimistic. Such a low value of Δf is at the very limit of what our sodium model predicts, even at low intensity and high altitude. A value of Δf = 7.3 kHz would cause equation 12 to predict the observed value of ΔV. This value of Δf, along with the reduced value of S_{0}, is used in the last row of Table 3, predicting when plugged into equation 8, close to the observed sensitivity of .
We believe that the data in general support the theoretical work we have presented, and that our sensitivity is worse than expected primarily because the intensity in the mesosphere was lower than expected, the signal collection efficiency was worse than expected, and the discrimination against Rayleigh scattered light was worse than expected. Also, the noise properties of the detector can be characterized by a reduced quantum efficiency.
All of these issues relate to the optimization of optics—the launch optics, collection optics, and detector. In a better optimized and better characterized system, we expect that actual sensitivity would much more closely match theory.
2.4.2 Improving Sensitivity and Resolution
The sensitivity, measured in , is important for improving the system's ability to measure a timevarying field. A system with a sensitivity of would provide a sensitivity approaching that of typical observatorygrade magnetometers but at a previously inaccessible altitude.
There are many ways in which our system could be improved. Two of the clearest are an improvement in the quality of the launched beam, so that we could be working at the optimum beam intensity, and an improvement of the efficiency of the collection optics. In both respects, our results are well below the state of the art, due to time and cost constraints, and it should be possible to make improvements.
In order to optimize the launch optics, we would need to interferometrically analyze the launch telescope and find if the wavefront were close enough to flat to create a diffractionlimited beam. Then we would replace components shown to be inadequate. In order to improve the efficiency of the collection optics, we would need to improve the throughput of each component in the receive telescope. We would also add an aperture or other mechanism to exclude Rayleigh light. We would begin by measuring total throughput, by observing the signal from a star of known magnitude. All these improvements were attempted, but we were unable to get useful results due to weather delays and the time and cost limits of this project.
Another improvement would be to use a stateoftheart guidestar laser. Such lasers, which are commercially available, have a power of 20 watts and are optimized for high photon return. (An important contributor to increased photon return is the addition of a second laser wavelength, with about 10% of total power, resonant with the D_{2}b line of sodium.) This wavelength repumps the sodium and increases Ψ to a value of 450 photons (s·sr·atom (W/m^{2}))^{−1} substantially above the value of 240 as shown in Table 1, which is valid for our singlewavelength laser. Holzlöhner et al. (2010) include theory and calculations of the effect of sodium atom repump. The expected return from such lasers is 36 × 10^{6} photons/s/m^{2} of incident light at the telescope aperture when the laser is operating continuously at its rated power of 20 watts (Drummond et al., 2007). These lasers are becoming more common at the large, advanced astronomical observatories, which have collection areas of 50 m^{2} or more.
A further improvement in sensitivity is theoretically possible if the laser used has a broadened linewidth. A laser with a linewidth roughly equal to the Dopplerbroadened sodium linewidth, that is, about 500 MHz, will interact with all of the sodium atoms, instead of only one velocity class. For a narrow linewidth laser, any velocitychanging collision by a sodium atom results in the atom being lost from the set of precessing atoms interacting with the laser beam, even if the angular momentum of the atom is unchanged by the collision. The linewidth Δf of the Larmor resonance is proportional to the rate at which atoms are lost from this set. A broad linewidth laser continues to interact with an atom after a change in its velocity, if the angular momentum state is unchanged.
It turns out that collisions with O_{2} change angular momentum, but collisions with N_{2} do not (Holzlöhner et al., 2010; Milonni et al., 1999). The atmosphere is about 20% O_{2}. Thus, even with a broad laser linewidth, about 20% of all collisions will remove the atom from the set that is precessing in phase with the laser. The Larmor linewidth assuming that only O_{2} collisions contribute is expected to be about 20% of what it would be if all collisions contributed to linewidth.
A pulsed laser with broadened linewidth has been proposed and patented (Kane et al., 2014; U.S. Patent No. 9,001,853, 2015). The proposed design has the additional advantage of greater power in a pulsed mode. A disadvantage of the broadened linewidth is that the optimum mesospheric intensity is increased, meaning that the beam has to be smaller in the mesosphere or the laser increased in power.
Table 4 summarizes three design options for remote magnetometers with improved sensitivity. All assume improved launch beam quality and collection efficiency and assume that the beam intensity in the mesosphere is optimized. A detector with 0.27 counts per incident photon noiseeffective quantum efficiency is assumed, such as a photomultiplier tube or an improved MPPC. An optimized 20watt guidestar system with Ψ = 450 photons (s·sr·atom (W/m^{2}))^{−1} is assumed throughout Table 4, other values are as from Table 1.
System  Current with improved optics and laser  Case 1 with broadened linewidth laser  Case 1 with large receive telescope 

Telescope aperture  1.8 m^{2}  1.8 m^{2}  50 m^{2} 
Laser average power  7 watts  20 watts  7 watts 
Laser linewidth  <1 MHz + D_{2}b  500 MHz + D_{2}b  <1 MHz + D_{2}b 
Detected signal, effective  3,095 ms^{−1}  8,842 ms^{−1}  85,964 ms^{−1} 
Resonance height H  20%  20%  20% 
Resonance width Δf  11 kHz  2.2 kHz  11 kHz 
Sensitivity using equation 8 
The first option would make use of a stateofthe art cw singlefrequency laser, with a telescope similar to what we used. Laser power is 35% of 20 watts, because the average power launched from the 20watt laser is reduced by the 35% duty cycle of the pulse format. The second option assumes that the laser is a broad linewidth, naturally pulsed design, so that no duty cycle penalty is incurred. The third option assumes that the stateoftheart cw laser is used with a large telescope, the type that is typically equipped with a guidestar laser.
Fan et al. (2016) have made an interesting and theoretically sound proposal for avoiding the dutycycle penalty. Instead of pulsing the laser near the Larmor frequency, they would reverse the sense of the circular polarization at that frequency. The only drawback is the need for fast polarization modulation of the highpowered 589nm beam.
2.4.3 Spatial Resolution
According to Table 2, the value of f_{L} changes by 1.3 kHz along the path of the laser beam through the layer. The sodium atoms at lower altitude have a higher value of f_{L} and a larger linewidth Δf.
For Z = 106 km, L = 10 km, and x = 10 m, Y = 106 m. With our current spot size of 2.6 m, these values would lead to four distinct resolvable spots where the magnetic field could be measured independently. Of course, the return signal from each spot would be only onequarter of the total signal, increasing sensitivity by a factor of 2 at each range.
Another way to get altitude resolved measurements would be to lower the PRF to a subharmonic of the Larmor frequency. The spin relaxation time (250 μs) is much longer than 10 times the Larmor period (32 μs). A PRF of 31.8 kHz would produce only one pulse in the sodium layer at a time, and the altitude detected could be determined by timegating the detector.
Either of the two methods could eliminate the undesirable Rayleigh backscatter photon flux.
3 Conclusions
We have demonstrated a remote magnetometer that has measured the magnetic field in the Earth's mesosphere, using sodium atoms at a 106km distance from our instrument. A pulsed laser with average transmitted power of 1.33 watts illuminated the sodium atoms, and the field was inferred from backscattered light collected by a telescope with a 1.55mdiameter aperture. The magnetic field was measured with a sensitivity characterized by an equivalent noise spectral density of . At this sensitivity, a measurement based on 1 hr of data averaging would have an expected RMS error due to noise of 2.7 nT. The value of magnetic field inferred from our measurement, 45,441 nT, is consistent with a value estimated based on the Earth's modeled field shape to within a fraction of a percent.
Projected improvements in optics could lead to sensitivity of , and the use of advanced lasers or a large telescope could approach sensitivity. All experimental and theoretical sensitivity values are based on a 60° angle between the laser beam axis and the magnetic field vector; at the optimal 90° angle sensitivity would be improved by about a factor of 2. But we were limited to a 60° angle at the Kuiper Telescope, so we consistently use that value. The optimal 90° angle may only be possible in the tropics; 60° or better should be possible over most midlatitude locations.
This technique works by tracking a resonance in the return signal that appears when the sodium atoms are illuminated by a circularly polarized laser pulsing at the frequency of the sodium atom precession, typically in the range 200–350 kHz. For the singlefrequency laser we used the height and width of the resonance reach an optimum when the average intensity of the pulsed laser light at the atoms is about 0.5 W/m^{2}. Our system was below this optimum, closer to 0.1 W/m^{2}.
Because neither satellites nor balloons can reach this altitude, this technique creates a unique window into an interesting zone of the geomagnetic field, on the lower edge of the ionosphere. Only rockets have been able to reach this zone previously and have only made transient measurements, not sustained measurements needed to yield the geomagnetic background noise (Sesiano & Cloutier, 1976; Volland, 1984).
Sodium layers exist in other planets and moons in our solar system. It should be possible to use this technique from a satellite in orbit, provided that the orbit is low enough to create a reasonably small beam in the sodium layer without the use of large optics. The fact that this technique is ultimately a measurement of frequency means that very long averaging times can be used, since frequency references can be extraordinarily stable. Thus, it is possible to imagine useful data resulting even at extremely low levels of return signal or during the daytime with higher background.
Acknowledgments
The results presented in this paper rely on data collected at the Tucson Magnetic Observatory. We thank the U.S. Geological Survey for supporting its operation and INTERMAGNET for promoting high standards of magnetic observatory practice (www.intermagnet.org), and the National Geomagnetism Program (observatory location Tucson, AZ): http://geomag.usgs.gov/monitoring/observatories/tucson/, http://geomag.usgs.gov/map/#realtime, and http://geomag.usgs.gov/plots/ We would also like to acknowledge the University of Arizona, Steward Observatory, for the use of the Kuiper Telescope and the 589nm FASOR coherent light source. Special thanks go to Jim Grantham and Steve Bland of the Stewart Observatory Mount Operations Team for making some last minute unscheduled telescope configuration changes, and finally Duncan Reed for weathering several cold nights airplane spotting. This material is based upon work supported by the U.S. Navy, Office of Naval Research  Code 321 under contract N0001411C0314 and N0001414C0110. The Final Report for the earlier contract provided estimates for both the sensitivity of the laserbased remote measurement under various conditions, and for the correlation that would be observed between highaltitude data and surface data. The authors feel the most pertinent figures in this paper are Figures 9 and 12. The original raw data for these figures are available for download in two Excel spreadsheets at https://wp.optics.arizona.edu/adaptiveoptics/wpcontent/uploads/sites/68/2018/03/figures.zip. Data for other figures discussed in the analysis can be derived from these data. We also thank the two reviewers for thoroughly reading our paper and making a few suggested changes to help clarify some fine details.