Full profile incoherent scatter analysis at Jicamarca

Incoherent scatter data from a hybrid longpulse/double-pulse experiment at Jicamarca are analyzed using a full-profile analysis similar to the one implemented by Holt et al.(1992). In this case, plasma density, electron and ion temperatures, and light ion composition profiles in the topside are estimated simultaneously. Full-profile analysis is crucial at Jicamarca, since the long correlation time of the incoherent scatter signal at 50 MHz invalidates conventional gated analysis. Results for a 24 h interval in April of 2006 are presented, covering altitudes through 1600 km with 10 min time resolution, and compared with results from the NRL ionospheric model SAMI2. The analysis provides the first comprehensive assessment of ionospheric conditions over Jicamarca at sunrise as well as the first 24-h record of helium ion layers. Possible refinements to the experiment and the algorithm are discussed.


Introduction
The Jicamarca Radio Observatory has design features that distinguish it from other incoherent scatter radars (ISRs) including its 50 MHz operating frequency. Jicamarca's long wavelength make it relatively immune to finite Debye-length effects that would otherwise limit its usefulness as an ISR at high altitudes. Ionospheric profiles have been measured at altitudes exceeding 6000 km in the past (Farley, 1991). Backscatter from neutral density irregularities is stronger at VHF than at higher frequencies, and Jicamarca is the only ISR able to observe echoes continuously from the troposphere through the mesopause (Woodman and Guillén, 1974;Maekawa et al., 1993). So-called "type 1" and "type 2" elec-Correspondence to: D. L. Hysell (dlh37@cornell.edu) trojet echoes have comparable strengths at 50 MHz, making it a suitable frequency for studying both (Balsley and Farley, 1971). The frequency is also optimal for studying echoes from plasma irregularities due to equatorial spread F, which would be weaker at higher frequencies but suffer more from refraction and scintillation at lower frequencies.
Disadvantages of 50 MHz compared to higher frequencies include a higher sky noise temperature that varies from about 4000-40 000 K, depending on sidereal time, and limits the sensitivity of ISR experiments. Furthermore, Jicamarca's frequency implies long correlation times (of the order of 0.5 ms) for incoherent scatter echoes observed obliquely to the geomagnetic field. Relatively long pulses must therefore be used to probe the ionosphere with satisfactory sensitivity and frequency resolution.
The data arising from long-pulse ISR experiments are voltage autocorrelation measurements or "lag products". The lag products are related to the autocorrelation function of the density fluctuations in the scattering volume, but the effective size of that volume is different for different lags (see for example Turunen, 1986). The volume for the shortest lag (the zero lag) is determined by the pulse length, while the volume for the longest lag is determined by the filter impulse response time, which is typically much shorter. Ad hoc summation rules can be used to degrade the resolution of some of the lag products in an attempt to equalize them, but uniform resolution cannot be obtained in practice. Furthermore, poorly-resolved lag products are unacceptable. Incoherent scatter theory is nonlinear, superposition does not hold, and parameter estimates based on spatially averaged autocorrelation functions can be unrepresentative.
For most of its history, Jicamarca employed a double pulse technique for incoherent scatter experiments (Farley, 1969b;Pingree, 1990). Pairs of short pulses with opposite circular polarizations were transmitted, and lag product matrices were assembled from the backscatter, one lag at a time. The size of the scattering volume was set by the pulse lengths Published by Copernicus Publications on behalf of the European Geosciences Union. rather than the lag spacing, and problems associated with inhomogeneity were thereby avoided. The double pulse experiment was also relatively immune to clutter from the ground, the equatorial electrojet, and other sources. Faraday rotation measurements using phase information from the orthogonal polarizations yielded an independent measurement of absolute electron density for calibration purposes (Farley, 1969a). However, the double pulse experiment was very slow (in terms of the rate of accumulation of information pertaining to the probability density function of the desired parameters as defined by Lehtinen, 1986) and did not fully exploit the average power limitations of the transmitters. Given reasonable incoherent integration times, the experiment could not yield useful ionospheric profiles at protonospheric heights. Nevertheless, it was a useful experiment and remains so for studying the valley and F regions.
More recently, a hybrid experiment combining the double pulses with an alternating code was implemented at Jicamarca (Hysell, 2000). The alternating code permits the estimation of the lag product matrix with range resolution set by the bit length of the code rather than the pulse length (Lehtinen, 1986;Lehtinen and Häggström, 1987;. The range resolution of all the lags is uniform except for the zero lag. The alternating code thus resolves many analysis problems associated with inhomogeneity. The alternating code is also much faster than the double pulse, providing multiple lag product estimates from each pulse and fully utilizing the duty cycles of the transmitters. Given reasonable incoherent integration times, the hybrid experiment yields useful ionospheric profiles well into the protonosphere. Error analysis for the alternating code is simpler than for a long-pulse experiment since the error covariance matrix of the former is diagonally dominant . A price must be paid in the form of increased self-clutter, however. Another promising strategy for coping with radar range ambiguity was suggested by Lehtinen (1986), implemented by Holt et al. (1992), and further investigated by . Rather than manipulating the lag products from a long-pulse experiment, constructing approximations of autocorrelation functions, and fitting the results for ionospheric parameters one range at a time (so-called "gated" analysis), these investigators proposed the simultaneous estimation of complete parameter profiles directly from the lag product matrices as they are measured. The technique is termed "full profile analysis" and belongs to a class of problems in statistical inverse theory. Inverse theory is applied widely in other branches of geoscience including seismology, geodesy, hydrology, crystalogrophy, and geology as well as in radar, radio astronomy, and inverse scattering. It seeks to find model parameters (ionospheric state variables in this case) that are consistent with whatever data are available within statistical confidence limits as well as with other a priori expectations.
In this paper, we investigate the application of full profile analysis to long-pulse incoherent scatter measurements at Jicamarca. We are interested in recovering plasma density, electron and ion temperature, and ion composition profiles in the F region and topside. The paper begins with a description of the hybrid double-pulse/long-pulse radar experiment at Jicamarca. After reviewing the forward problem that specifies how the measurements depend on the sought-after parameters, we describe the inversion algorithm which differs from that of Holt et al. (1992) in some respects. Finally, experimental results are presented and compared with simulation results from the NRL ionospheric model SAMI2, and potential improvements are discussed.

Experiment description
The experiment is performed using the entire main antenna at Jicamarca in its "on axis" or uniform phase configuration. The angle the main antenna beam makes with respect to the geomagnetic field on this position is slightly less than two degrees off perpendicular near the F peak and decreases with altitude. Left-and right-circularly polarized signals are used for transmission and reception. A separate transmitter excites either polarization through hybrid networks with a peak power of about 1 MW each.
The interpulse period is 40 ms, during which time three pulse sequences are emitted. The first is a 1.6 ms (240 km) long pulse, which is transmitted on a single circular polarization. Following that, two double pulse pairs with individual pulse lengths of 100 µs (15 km) are transmitted. The pulse separation is n times 200 µs (30 km), where n takes on all integers from 0 through 10. Opposite circular polarizations are used for the first and second pulse in the pair. The phase of the second pulse is inverted in the second pair to help suppress crosstalk in the final analysis.
Reception is performed using digital receivers sampling at 50 µs intervals. Additional boxcar filtering of the data reduces the sample rate to once in 100 µs. The interval between the long pulse and the first double pulse pair transmission is roughly 26 ms or about 3900 km. Noise estimates are drawn from longest range delays.
Clutter from the ground, the equatorial electrojet, and 150 km echoes render the long-pulse data useless in range gates below about 450 km. (This is true for all the pointing positions available at the observatory.) Ionospheric parameters below this range are estimated entirely from the double pulse data. The double pulse is relatively immune to clutter, since only one or two lags of the lag profile matrix for any given range gate are typically affected. These are easily identified and discarded. Error analysis for the double pulse experiment is elementary since errors for different lags of the autocorrelation function are statistically independent (Farley, 1969b;Pingree, 1990). The double pulses also provide an independent estimate of absolute electron density through Faraday rotation. Electron density profiles inferred from Faraday rotation are used to normalize double-pulse Ann. Geophys., 26, 59-75, 2008 www.ann-geophys.net/26/59/2008/ power profiles which, in turn, are used to normalize power profiles recovered from full profile analysis of long-pulse data. While crosstalk between the polarizations can lead to various kinds of errors, this can be removed to first order using phase flipping, as mentioned above. Whereas double-pulse data can be analyzed using conventional gated analysis, long-pulse data must be analyzed using full-profile analysis and inverse methodology for reasons outlined in the introduction. In practice, this entails solving the forward model iteratively until a cost function composed of the model prediction error (essentially the chi-squared statistic) and a combination of model metrics is minimized. The following section specifies the forward model, which includes both the so-called radar ambiguity function and incoherent scatter theory. A more general treatment of the complete radar instrument function can be found in Woodman (1991) and .

Parameter estimation and the forward model
Free electrons in the ionosphere can be viewed as a uniform background, which does not scatter radio waves, with superimposed electron density fluctuations N, which do. A radar signal radiated at time t−2r/c and scattered off density fluctuations at time t−r/c and radar range r will be received at the radar at time t. The relative phase of the transmitted and received signals will be k · r radians, where k is the scattering wavevector (the difference of the incident and scattered wavevectors) and r is the scattering coordinate. In light of the first Born approximation, summing over the radar scattering volume gives the total signal applied to the receiver x(t) = d 3 r e ik·r s(t − 2r/c) N(r, t − r/c) Following Farley (1969b), we exclude constants associated with the transmitter power, antenna gain, electron scattering cross section, and target range that are not relevant to the current analysis. Here, s(t) is the transmitted pulse shape which defines the range extent of the scattering volume. If the impulse response function of the receiver is h(t), then the corresponding receiver output at the sample time t s will be The receiver output can be regarded as a zero-mean Gaussian random variable, with information contained in the second-order statistics: y(t s2 )y * (t s1 ) = dt 1 dt 2 d 3 r 1 d 3 r 2 e ik·(r 2 −r 1 ) · N (r 2 , t 2 − r 2 /c) N * (r 1 , t 1 − r 1 /c) The following variable substitutions will be expedient.
Because electron density fluctuations in thermal equilibrium have a short correlation length, we can take r 2 −r 1 =r ∼ 0 everywhere in this expression except in the exponent. This yields We can now perform the d 3 r integral, noting that the term in the angle brackets is the spatio-temporal autocorrelation function of the density irregularities and that the integral in question is then its spatial Fourier transform. The result is the temporal autocorrelation function for electron density irregularities with a wavevector k or ρ(k, t2, t1; r). This term is stationary and depends only on τ ≡t 2 −t 1 . Consequently. another change of variables can now be made: Finally, defining the one-and two-dimensional radar ambiguity functions (Lehtinen, 1986;Nygrén, 1996) leads to the following compact result y(t s2 )y * (t s1 ) = dτ dr ρ(k, τ ; r)W ts1,ts2 (τ, r) The quantity on the left is the lag product matrix, which constitutes the measurements for the experiment with associated noise and statistical uncertainty. The ambiguity functions are known. The autocorrelation function is related to the lag products through a two-dimensional convolution with through incoherent scatter theory. Estimating the state parameters from lag product measurements is a problem in inverse theory. Figure 1 shows the sixteen ambiguity functions plotted side-by-side describing a long-pulse experiment with a pulse width 16 times the receiver sampling interval. In this case, the pulse shape s(t) was made to have a linear downward slope associated with a linear transmitter power droop.

Error analysis
The data used for the analysis are the real parts of the lag product matrices formed from the long-pulse experiment. The imaginary parts of all the lag products are small given the largest line-of-sight drifts observed at Jicamarca. Denoting ρ 12 as the lag product formed by summing the product of voltage samples labeled '1' and '2', the elements of the data covariance matrix bridging ρ 12 and ρ 34 can be shown to be (see appendix) where k is the number of statistically independent samples. Note that the sample indices can be repeated (i.e. (2) provides the autocovariances as well as the covariances).
To the extent that noise and self-clutter contribute to any of the lag products, they also contribute to the error covariance matrix in accordance with (2). No special calculations or partitioning of the lag products into signal, self-clutter, and noise components need be performed for the error estimator to hold in this kind of analysis. This represents a significant bookkeeping simplification compared to traditional gated error analysis (Hysell, 2000). When noise estimates are subtracted from the lag product matrices, the subtraction must however take place after C d is calculated.
It is worth comparing the error performance of a longpulse experiment (using (2), no summation rules) and a randomized 16-bit alternating code experiment (using the formulas developed in Hysell (2000)). A comparison based on a model incoherent scatter autocorrelation function roughly representative of the equatorial topside is shown in Figure 2. the two-dimensional ambiguity function. The autocorrelation function itself is related to the state parameters of the ionosphere (density, temperature, composition, drifts, etc.) through incoherent scatter theory. Estimating the state parameters from lag product measurements is a problem in inverse theory. Figure 1 shows the sixteen ambiguity functions plotted side-by-side describing a long-pulse experiment with a pulse width 16 times the receiver sampling interval. In this case, the pulse shape s(t) was made to have a linear downward slope associated with a linear transmitter power droop.

Error analysis
The data used for the analysis are the real parts of the lag product matrices formed from the long-pulse experiment. The imaginary parts of all the lag products are small given the largest line-of-sight drifts observed at Jicamarca. Denoting ρ 12 as the lag product formed by summing the product of voltage samples labeled "1" and "2", the elements of the data covariance matrix bridging ρ 12 and ρ 34 can be shown to be (see appendix) where k is the number of statistically independent samples. Note that the sample indices can be repeated (i.e. (2) provides the autocovariances as well as the covariances).
To the extent that noise and self-clutter contribute to any of the lag products, they also contribute to the error covariance matrix in accordance with (2). No special calculations or partitioning of the lag products into signal, self-clutter, and noise components need be performed for the error estimator to hold in this kind of analysis. This represents a significant bookkeeping simplification compared to traditional gated error analysis (Hysell, 2000). When noise estimates are subtracted from the lag product matrices, the subtraction must however take place after C d is calculated.
It is worth comparing the error performance of a longpulse experiment (using (2), no summation rules) and a randomized 16-bit alternating code experiment (using the formulas developed in Hysell, 2000). A comparison based on a model incoherent scatter autocorrelation function roughly representative of the equatorial topside is shown in Fig. 2. For simplicity, the model ACF is made purely real and positive definite. The comparison is made in the intermediate signal-to-noise ratio (SNR) limit, where we take the SNR for the zero lag long-pulse estimator to be unity, typical for the most distant ranges considered here. Given this value, selfclutter and noise contribute about equally to the covariances of the alternating code whereas the long pulse is noise dominated in all but the longest lags.
Given 16 samples over the pulse length, 136 distinct lag products can be formed. Figure 2 depicts the corresponding relative error covariances in a grayscale format spanning 40 dB of dynamic range. The same absolute scale is used for both plots, with the scale maximum set to correspond to the terms on the main diagonal of the alternating code plot. Normalization is done by dividing each element of the matrices by the product of the maximum possible values of the two lag products involved (post noise power removal) in that element. Since zero lag estimates are not provided by the alternating code experiment, the relevant rows and columns of the matrix are blanked.
Errors in the long-pulse experiment are more highly correlated than in the alternating code experiment. If Fig. 2 was replotted to show only the upper 20 dB of the scale, almost everything but the main diagonal of the alternating code plot would vanish. Along the main diagonals, however the alternating code autocovariances are larger than the long-pulse autocovariances by as much as 20 dB in the terms representing short lags. These terms are of primary interest in the topside. These terms suffer less self-clutter from a long pulse than from an alternating-coded pulse. If Fig. 2 was replotted to show only the upper 20 dB of the scale, almost everything in the long-pulse plot would vanish.
We have not performed a sensitivity analysis and cannot fully assess the relative merits of long pulses and alternatingcoded pulses on the full-profile algorithm, which could utilize data from either. Because the eigenvalues of the longpulse error covariance matrix are all smaller than the terms of the main diagonal of the alternating code matrix and that the median eigenvalue is 20 dB smaller, we undertook the present effort using long-pulse data. The merits of this approach are still uncertain, however, since the long-pulse data inverse problem is apt to be more poorly conditioned than alternating-code problem (see below).
Throughout most of the remaining analysis, we neglect the off-diagonal components of the error covariance. Huuskonen Ann. Geophys., 26, 59-75, 2008 www.ann-geophys.net/26/59/2008/ and  found that this can lead to a factor of two underestimate in the errors of the plasma parameters in a gated-analysis long-pulse experiment. It is not at all clear what effect it has on the full-profile analysis. The problem is beyond the scope of the present study but requires further investigation.

Incoherent scatter theory
The theory that relates the autocorrelation function ρ(k, τ ; r) to the state variables of the ionospheric plasma (electron density, electron and ion temperatures, ion composition, drifts, etc.) is well known and need not be reviewed here (see for example Dougherty and Farley (1960); Farley et al. (1961); Dougherty and Farley (1963); Swartz and Farley (1979); Fejer (1960aFejer ( ,b, 1961; Salpeter (1960Salpeter ( , 1961; Hagfors (1961); Woodman (1967). Application of the theory at Jicamarca is complicated by the small magnetic dip angle, however. A considerable body of work exploring the theory at small aspect angles has emerged in recent years. Here, we highlight findings relevant to the our forward model. It is possible to measure line-of-sight plasma drifts very accurately at Jicamarca when the beam is pointed perpendicular to the geomagnetic field and the incoherent scatter correlation time becomes very long, permitting pulse-to-pulse analysis (Woodman and Hagfors, 1969). However, Kudeki et al. (1999) first noted that the incoherent scatter spectrum at small aspect angles is even smaller than the standard theory predicts. (See also Kudeki and Milla, 2006.) At the same time, it had long been known that the standard theory failed to predict the correct electron-ion temperature ratio for finite but small aspect angles, of the order of a few degrees or less (Pingree, 1990). Sulzer and Gonzalez (1999) were able to account for both of these phenomena (qualitatively for the former and quantitatively for the latter) by including the effects of electron Coulomb collisions in the calculation of the electron admittance. Their calculations were based on Monte Carlo simulations. Summarizing these findings, Aponte et al. (2001) went on to show that the new Coulomb collision corrections bring temperature ratio observations made a few degrees off perpendicular at Jicamarca into compliance with expectation. Those corrections, which have been made available in tabular form (M. Sulzer, personal communication) are incorporated in the present analysis.
A second issue involves the relationship between plasma density, the radar scattering cross-section, and power profile measurements. At small aspect angles, the incoherent scatter cross section depends on the electron and ion temperatures and the plasma density in a complicated way not predicted by the approximation for large aspect angles given by Farley (1966) employed widely at other ISR facilities. Milla and Kudeki (2006) have examined the effect in detail, including for the first time the effects of electron Coulomb collisions in their calculations. Extending the work by Woodman (2004) through incoherent scatter theory. Estimating the state parameters from lag product measurements is a problem in inverse theory. Figure 1 shows the sixteen ambiguity functions plotted side-by-side describing a long-pulse experiment with a pulse width 16 times the receiver sampling interval. In this case, the pulse shape s(t) was made to have a linear downward slope associated with a linear transmitter power droop.

Error analysis
The data used for the analysis are the real parts of the lag product matrices formed from the long-pulse experiment. The imaginary parts of all the lag products are small given the largest line-of-sight drifts observed at Jicamarca. Denoting ρ 12 as the lag product formed by summing the product of voltage samples labeled '1' and '2', the elements of the data covariance matrix bridging ρ 12 and ρ 34 can be shown to be (see appendix) where k is the number of statistically independent samples. Note that the sample indices can be repeated (i.e. (2) provides the autocovariances as well as the covariances).
To the extent that noise and self-clutter contribute to any of the lag products, they also contribute to the error covariance matrix in accordance with (2). No special calculations or partitioning of the lag products into signal, self-clutter, and noise components need be performed for the error estimator bookkeeping simplification compared to traditional gated error analysis (Hysell, 2000). When noise estimates are subtracted from the lag product matrices, the subtraction must however take place after C d is calculated.
It is worth comparing the error performance of a longpulse experiment (using (2), no summation rules) and a randomized 16-bit alternating code experiment (using the formulas developed in Hysell (2000)). A comparison based on a model incoherent scatter autocorrelation function roughly they incorporate the effects of Coulomb collisions analytically with the introduction of a heuristic operator in the socalled Gordeyev integral which, when evaluated numerically, reproduces the results of Sulzer and Gonzalez (1999 Jicamarca power profiles to density profiles, finding good agreement with Faraday rotation measurements. We adopt the same methodology here, making use of Faraday rotation only for absolute density calibration.

Inverse methodology and regularization
Full profile analysis is an example of a discrete inverse problem of the form G(m)=d, where d is a column vector representing known data, m is a column vector representing unknown model parameters, and G is the known theory relating them. In this instance, the data are lag profile estimates, the model parameters are plasma density, temperature and composition specified at various altitudes, and the theory is incoherent scatter theory combined with (1) above. If the model is linear or its estimate m est is sufficiently close to the solution that the theory can be linearized about m est , then the problem becomes simply Gm est =d, where d ∈ R n , m ∈ R m , and G ∈ R n×m . The apparent simplicity of even this problem is deceptive, and its solution in the presence of noisy data almost inevitably necessitates recourse to statistical inverse methods. For an introduction to inverse methodology, consult for example Menke (1984), Tarantola (1987), Kaipio and Somersalo (2004), and Aster et al. (2005).
We are here interested in the overdetermined or mixeddetermined problem, with n>m≥ rank(G). In such problems, the column space of G does not span R n , and there exists a data nullspace (space occupied by nontrivial solutions of the problem G t x=0) from which components of d cannot be drawn. In general, this makes it impossible for any model to solve the inverse problem exactly, and the problem becomes one of finding a model estimate which is in some sense best. The least-squares solution is the one that minimizes the L2 norm of the model prediction error ||Gm−d|| 2 2 . If the data are normally distributed and the data covariance matrix C d is known, then the weighted least squares estimator can be used to find the model that minimizes (Gm−d) t C −1 d (Gm−d). In the absence of any other information, the weighted least squares solution is the most likely solution in a frequentist sense. The estimator that minimizes the χ 2 statistic can be expressed as Despite their widespread use, linear and nonlinear least squares strategies suffer from important deficiencies that render them suboptimal for many problems. In the case of mixed-determined problems for which m > rank(G), there exists a model nullspace (space occupied by nontrivial solutions of the problem Gx=0), rendering the least-squares solution non-unique. Moreover, we do not expect the solution to the inverse problem to have uniformly small χ 2 . The χ 2 statistic is itself a random variable with expectation governed by the number of degrees of freedom in the system. In both the under-and mixed-determined problems, there will generally exist a large number of solutions with acceptable values of χ 2 . While being statistically indistinguishable, these solutions can be markedly different, particularly if the system is poorly conditioned and unstable. Full profile analysis represents one such system. In simply minimizing χ 2 , one is likely to end up merely "fitting the noise." The importance of conditioning can be appreciated by considering the generalized inverseĜ which solves the problem m est =Ĝd. The Moore-Penrose generalized inverse can be constructed using singular value decomposition (SVD) or generalized singular value decomposition (GSVD) ( Van-Loan, 1985). In SVD, there is the factorization where U and V are constructed from the column and row space of G, respectively, and S is a diagonal square matrix composed of the nonzero singular values of G. The factorization works for over-, mixed-, under, and even-determined problems, and results in the least squares estimate in the first case. Performing the decomposition reveals that m est is formed from the row space of G, with components weighted by the reciprocal of the associated nonzero singular values. Any noise in the data propagates through to the model with the same weighting. Problems with the smallest singular value much smaller than the largest are termed poorly conditioned since they will tend to exhibit excessive noise amplification. Model solutions to such problems may tend to oscillate wildly even though predicting the data accurately. The characteristic is endemic to inverse problems and generally arises when attempting to measure more parameters than the underlying measurement technique warrants. Such problems are often approached casually using simple interpolative strategies, often with poor results.
Regularization attempts to remedy the problem by attenuating the components of m est most subject to noise amplification. This can be done using SVD or GSVD by multiplying the singular values of G (or rather their reciprocals) by "filter factors", imposing a floor on their magnitude. Different filter factors bring about different results, typically minimizing the L2 norm of a model metric of the form ||Lm|| 2 2 . Choosing the identity for L minimizes the size of the model in a strategy termed "minimum length" and "model simplicity". It is also referred to as Tikhonov regularization. Other forms for L can be used to minimize the first or second derivative of m globally. This is referred to as generalized Tikhonov regularization.
In what follows, we will consider regularization strategies based on model simplicity. The following strategies can be shown to be equivalent for the case L=I : 1. SVD using filter factors of the form f i =s 2 i /(s 2 i + α 2 ), where s i are the singular values and α is the so-called regularization parameter. The data covariance matrix can be incorporated by transforming and scaling G and d.
The model estimator that accomplishes this is m est =(G t C −1 d G+α 2 C −1 m ) −1 G t C −1 d d This strategy is termed "weighted damped least squares".

Conjugate gradient weighted least squares minimization
with an initial guess m est =0 and with early iteration termination consistent with some finite α.
4. Augmented least squares, seeking the least-squares solution to the augmented minimization problem:

Solving (for m est ) the characteristic equation
the result being the weighted damped least squares estimator above.
In addition to the frequentist perspective, inverse problems are often solved with Bayesian methodologies, in which the model parameters are themselves treated as random variables and apriori probabilities are explicitly formulated and included. It can be shown that the minimum length strategy is equivalent to the maximum a posterior probability (MAP) solution of a Bayesian strategy where the model parameters are assumed a priori to be normally distributed (Tarantola, 1987).
Just as regularization can be incorporated into linear least squares minimization algorithms through augmentation of G and d as illustrated in the fourth item above, so too can it be incorporated in nonlinear least squares algorithms through analogous matrix augmentation of the Levenberg Marquardt problem. In either case, the regularization parameter must be determined. While guides for selecting α exist (including generalized cross-validation and the L-curve method), the subjectivity associated with the choice is a weakness of regularization that remains an topic of research and debate.

Full profile algorithm
The algorithm is based on the one designed by Holt et al. (1992) with a few innovations. The most important of these are described below. Additional complexity stems from our attempt to infer temperature and composition profiles in addition to plasma density. Other factors peculiar to Jicamarca also enter.

Noise and debris removal
Sky noise power is estimated from a weighted ranking of the power estimates for the sample rasters of all the receivers. It is then subtracted from the double-pulse and long-pulse lag products. It is important to note that both the zero and adjacent first lags of the long-pulse data receive contributions from sky noise. Only the zero lags of the double-pulse data receive contributions.
Space debris poses special problems for Jicamarca since everything in low Earth orbit eventually passes over the equator. Most of the debris is in the Rayleigh regime at 50 MHz, meaning that the debris backscatter gain is always small but finite. The balance of the debris is observed between about 700-1000 km range. It is usually but not always slowly fading and persists for the beam transit time, typically a few seconds.
The most effective strategy for debris removal we have found involves ranking lag profile estimates computed from small datasets. We typically compute long-pulse lag products from 128 adjacent pulses representing about 5 s of experiment time. Sixteen sequential estimates made this way are stored and ranked by the size of their real parts. The two smallest and two largest of these are then discarded, with the remaining lag products used for analysis. The "outliers" are self selecting and so improve the statistics of the experiment even having been rejected.

Double pulse data processing
Double-pulse data are processed according to the procedure described in Pingree (1990). The use of digital receivers obviates the need for many of the data correction algorithms used in the past associated with receiver bias, gain imbalance, and orthogonality imbalance. Crosstalk is still corrected to first order with phase flipping. Lags of the ACF contaminated by ground clutter and coherent echoes are rejected on the basis of the relative strength of the echoes seen on the two circular polarizations. The ACFs are analyzed using conventional nonlinear least-squares fitting. We generally fit for electron and ion temperature but assume a purely O + plasma. The results are power and temperature profiles along with an absolute density profile derived from Faraday rotation. Useful data are obtained starting at about 180 km altitude, whereas clutter dominates below. These are the only useful data we have below about 450 km altitude.

Power profile recovery
Each of the measured long-pulse lag products represents the convolution of the radar range-time ambiguity function and the autocorrelation function of the scattering medium. The zero lag is no exception. However, since the ACF is nearly uniform over the small delays τ entering into (1), we can regard the zero lag as the one-dimensional correlation of a www.ann-geophys.net/26/59/2008/ Ann. Geophys., 26, 59-75, 2008 "true" power profile with the emitted radar pulse shape to a good approximation. Here, "true" means the profile that would be measured ideally if short pulses were used. We can therefore recover the former by deconvolving the latter from the zero lag profile. The idea is the same as that employed by for Arecibo long-pulse analysis except that we apply it to the zero lag only (R. Nikoukar, personal communication, 2007), the one-dimensional approximation failing to hold generally for other lags in our case. Deconvolution is performed using Tikhonov regularization. This is accomplished with augmented linear least squares minimization, the fourth enumerated strategy above. Here, L is the identity matrix, and the regularization parameter α is set to a constant value that suppresses oscillations in all of our profiles. Minimization is performed using a non-negative least squares (NNLS) routine that prevents the occurrence of negative results anywhere in the deconvolved profile (Lawson and Hanson, 1987). This strategy enhances the stability and robustness of the methods and permits smaller regularization parameters than would otherwise be possible. It is noteworthy that the droop in transmitter power that occurs during long pulse transmission produces a forward model that is better conditioned than would be a uniform pulse.
Coherent scatter from altitudes up to about 180 km render long-pulse samples below about 180+240=420 km useless for incoherent scatter analysis. In practice, finite receiver recovery time makes this figure more like 450 km. We therefore restrict the full-profile analysis to lag products corresponding to ranges greater than or equal to 450 km. However, forward modeling of the lag products starting at 450 km requires an estimate of the "true" power profile at altitudes below 450 km, something that cannot be recovered very accurately from long-pulse data deconvolution. We therefore extend the "true" power profile downward by using the double pulse power profile at lower altitudes. The strategy works well except when the double-pulse data are contaminated by coherent scatter from equatorial spread F or when the backscatter is too weak for accurate power estimation.
An approximate electron density profile can be constructed by multiplying the true power profile by the square of the radar range and scaling the result for best agreement with electron density profiles measured using Faraday rotation in the least-squares sense. In the event of unequal electron and ion temperatures, the profile must be further corrected as described earlier in the incoherent scatter section of the paper.

Forward model
The "true" power profile completes the forward model, since the product of this profile with the normalized (max. value of unity) theoretical incoherent scatter autocorrelation function gives ρ(k, τ ; r) in (1). Given specifications for the ionospheric density, temperature, and composition profiles along with the radar range-time ambiguity function, (1) can be used to predict the absolute lag product matrix. In practice, the integrals in (1) are evaluated discretely using a trapezoidal rule. The proxy electron density profile is a sufficiently good approximation for n e for all forward modeling intents and purposes and need only be corrected at the conclusion of the analysis.
For the forward model, electron temperature, the electron/ion temperature ratio, and the hydrogen and helium ion fractions are specified coarsely at 20 uniformly-spaced ranges spanning the solution space. These become the knots used for cubic B-spline interpolation (De-Boor, 1978). The interpolated functions thus created are twice-differentiable continuous and provide a means of specifying the ionospheric parameter profiles with arbitrarily fine resolution for forward modeling. Computing the B-spline coefficients is fast and efficient, and experimentation has shown that the topside profiles measured with conventional gated analysis at Jicamarca can be well represented with 20 or fewer knots.

Regularization
The model parameters, the 4×20 knot values, are calculated using augmented nonlinear least squares minimization following a nonlinear generalization of the fourth regularization strategy enumerated above. The minimization algorithm is conventional Levenberg-Marquardt. The parameters are initialized with values resulting from an extensive grid search when a large family of representative trial profiles are generated and a χ 2 statistic calculated for each. It is generally possible to find initial profiles this way for which χ 2 /m<5, m being the number of parameters. Levenberg Marquardt iteration generally yields final values for χ 2 /m ∼2-3. Smaller values can be produced only at the cost of producing very irregular parameter profiles.
Augmentation in this case means introducing six additional cost functions that help to regularize the output of the algorithm. Each of these is multiplied by an adjustable weight. The weights have been set so as to inhibit results that are inconsistent with general expectations for the ionosphere. The cost functions are enumerated below.
The first three cost functions pertain to the global range averages of the second derivatives of the electron temperature, ion temperature, and hydrogen ion fractions, respectively. Suppressing these quantities inhibits spurious oscillations in Ann. Geophys., 26, 59-75, 2008 www.ann-geophys.net/26/59/2008/ the profiles associated with poor conditioning. Note that the net hydrogen fraction in our algorithm is defined by A similar function is used to express the helium ion fraction. The functions restrict the ion fractions to values between 0 and 1. It remains possible for the combined hydrogen and helium ion fractions to exceed unity, and so a fourth cost function was introduced to inhibit this. A fifth cost function is a very weak prohibition of helium ions overall. This is meant to prevent the algorithm from producing spurious He + layers that might otherwise come about as a consequence of any systemic errors in data acquisition and processing. Helium ions indicated by the full profile analysis must be pointed to by clear signatures in the lag product data. Finally, a sixth cost function is a prohibition against electron-to-ion temperature ratios less than unity. This is a condition that we expect never to occur in nature, and the cost function prevents the full profile algorithm from exploring that region of solution space.

Parameter estimation and error analysis
The algorithm iterates until meeting convergence criteria involving the residuals of the characteristic equation, the rate of change of the model parameters, and the rate of change of the data residuals. Confidence limits are placed on the model results according to the usual error propagation formula for the model covariance matrix: where J is the Jacobian matrix evaluated for the model solution vector. The same formula can be used to place confidence limits on the electron density profile calculated through one-dimensional deconvolution; in that case, the problem is linear, and G for deconvolution replaces J . Full profile analysis by its nature produces smooth profiles with variations often much smaller than their associated confidence limits. For that reason, the algorithm can be interpreted as an efficient interpolation mechanism designed to produce smooth profiles as demanded by Occam's Razor. This reflects the fact that the parameter errors are highly correlated; errors in any one parameter in a profile are shared by its neighbors. Smoothness should not be confused with confidence, and computed confidence limits should be taken seriously.

Experimental results
Twenty-four hours of hybrid double/long-pulse data were taken during 19-20 April 2006, at Jicamarca and processed using the full-profile algorithm. In the analysis, only the diagonal components of the error covariance matrices were used. While the formalism surrounding the calculation and incorporation of the complete covariance matrix is straightforward, the storage and computational demands are significant, and sparse matrix solvers are required in practice. Sparse matrix math has yet to be incorporated in our algorithm.
Processed data representing 10 min incoherent integration are shown in Figs. 2-5 which follow. Each figure is composed of five panels. From left to right, the first panel shows lag products for different range gates computed from double pulse data. Vertical lines are drawn through lag products deemed to be contaminated by coherent scatter or ground clutter. Measurements are represented by points with associated error bars. Theoretical fits to the data are represented by solid curves.
The second panel in each figure shows lag products computed from long-pulse data. Points with error bars and solid curves once again represent measurements and theoretical fits, respectively. The plots are normalized against the zero lag value for each range delay. Whereas the double pulse data were processed using conventional gated analysis, the longpulse data were processed using the full profile technique.
The third panel depicts electron density versus altitude. Three curves are shown. Above 450 km (below 600 km), electron densities derived from long-pulse (double pulse) data are plotted. The error bars associated with the respective analyses are also plotted. In addition, absolute electron density derived from Faraday rotation is indicated by the curve in green. This curve was used to normalize the double-pulse power measurement which was, in turn, used to normalize the long-pulse measurement where the two overlap.
At the magnetic equator, sources of clutter including meteor echoes, the equatorial electrojet, and the 150-km echoes generally contaminate incoherent scatter echoes below about 180 km during the day and 140 km at night. It is seldom appropriate to interpret echoes below these altitudes in terms of incoherent scatter theory. At night, coherent scatter from equatorial spread F contaminates the echoes at higher altitudes. Bottomside scattering layers are virtually always present for a few hours after sunset except during summer solstice and are another source of clutter and contamination. Full profile analysis is especially sensitive, since clutter at any altitude can invalidate the analysis at all altitudes.
The fourth panel in each figure shows electron and ion temperatures as red and black plotter symbols, respectively. These are essentially identical at night, in which case the symbols overlap. Error bars are shown. Lastly, the fifth panel shows light ion fraction. The black and blue plotter symbols represent hydrogen and helium ions, respectively. Below 450 km, a pure oxygen ionosphere is assumed, and no composition estimates are plotted. profile, and distinctly different topside density scale heights below and above the 50% H + transition height. The temperature profile is essentially isothermal, as expected for this time forward. Nothing about the full profile analysis constrains the topside temperatures to be uniform in altitude, however. That they are uniform and continuous with the temperatures derived from gated double pulse analysis below 450 km argues for the efficacy of the new analysis. Note that Jicamarca was running at reduced power for the first few hours of this experiment. This resulted in lag profiles with more scatter than at subsequent times. The layers in the power profile between 200-300 km altitude are weak bottom-type coherent scattering layers that obscure the shape of the bottomside, which typically has a gradient scale length of about 10 km after sunset. Figure 4 show postmidnight data. Now, the F peak resides at only about 250 km altitude, and the O + /H + transition has fallen to about 650 km. The H + profile varies more rapidly and linearly and with altitude than before. There is very little O + above 800 km now and for the remainder of the night. It is evident that the hydrogen ion fraction below 450 km al-titude is finite, but the double pulse data are not sufficiently robust to permit composition estimation, and the temperature estimates immediately below 450 km are suspect as a result.
Postsunrise conditions are illustrated by Fig. 5. Plasma density and composition are just beginning to respond to sunlight. However, electron and ion temperatures are substantially elevated, and the electron-ion temperature ratio is also discerningly greater than unity everywhere below about 600 km for reasons mentioned below. The accuracy of the double pulse temperature estimates is poor above the F peak, as indicated by the scatter in the fit temperatures just below 450 km, but the full profile results appear to agree with the double pulse data below 400 km. Previous analyses at Jicamarca were insufficiently robust to permit the estimation of unequal electron and ion temperatures simultaneously with full light ion composition, and Fig. 5 in particular represents a novel result for the observatory.  atively slowly with altitude again, and the topside plasma density gradient scale height varies more gradually with altitude than at any other time of day. The electron-ion temperature ratio is as high as 2 below the F peak where EUV absorption is strongest and remains greater than unity up to about 800 km. In future analyses, it should be possible to infer the neutral temperature from the temperature ratio over this altitude range (Aponte et al., 2001). Most remarkably, a thin helium ion layer with a peak density of about 10% is evident near 750 km altitude. Layers such as these appear intermittently early in the day but at no other times (see below).
The features of the individual profiles recur in the summary plot shown in in Fig. 7 for the 24-h run. Data between 03:00-06:00 LT when the galaxy was inside the antenna beam have been excised, as have intervals when residual echoes from space debris interfered with the analysis. While there is significant scatter in some of the parameters, particularly in the postsunrise temperatures, this is a consequence of low signal levels and consistent with our error analysis.
During the day, the electrons are at least somewhat warmer than the ions throughout most of the thermosphere. Two situations lead to drastic temperature differences. In the bottom-side F region where photoelectron production is most rapid, the electron temperature is elevated in a narrow band of altitudes, but the ions are cooled through strong coupling to the dense neutral atmosphere, and a large T e /T i temperature ratio results. The ratio is greatest between 13:00-14:00 LT. At sunrise, meanwhile, a large temperature ratio also arises over a wider altitude span in the lower topside. Here, the plasma density is very low, the electrons and ions are weakly coupled thermally, and the ion temperature again falls well short of the electron temperature. Low F layer densities at sunrise may also permit photoelectrons generated lower on the flux tube to penetrate to the apex points and heat the electrons in the topside. The role of photoelectrons in the thermal balance of the equatorial ionosphere will be examined in subsequent analyses. The temperature ratio diminishes by 09:00 LT with the increase in topside plasma density.
Throughout the day, the electron and ion temperatures in the thermosphere exhibit a positive density gradient. This is associated with the parallel thermal conductivity and plasma heating up to the apex point from the warmest strata in the off-equatorial bottomside. The plasma temperature appears to approach the temperature of the peak below the bottomside.  Figure 7 clearly shows a helium ion later between about 500-800 km altitude in agreement with Fig. 6. The layer is patchy, and the ion fraction seldom exceeds 10%. The layer occurs near but below the protonospheric transition height and is visible only early in the day. Figure 7 also shows that the protonospheric transition height only varies by about 250 km throughout the day. Helium ion layers have not been observed at Jicamarca for decades and never before as the output of an automated analysis. The diurnal behavior of the H + /O + transition height has been observed using alternating codes in the past, but only during high solar flux conditions when the excursions were much larger. In April of 2006, the 10.7 cm solar flux index was 70, and the transition height stayed remarkably low throughout the 24-h period.
The diurnal behavior of the O + /H + transition height is similar to that reported by González et al. (1992) based on in situ observations of the AE-E satellite at low latitudes during solar minimum. That study also found that the He + concentration had a daytime maximum during equinox, although the peak concentration reported was an order of magnitude smaller than what we find here. At middle latitudes, the He + ion fraction has been shown to reach 10-20% during solar minimum, but at night rather than during the day (González and Sulzer, 1996). The mechanics of helium ion layer generation has been modeled by Heelis et al. (1990). At middle latitudes, helium ions produced by photoionization move upward along flux tubes under the influence of the ambipolar electric field where they tend to concentrate at or above the O + /H + transition height. A reservoir of helium ions thus forms at altitudes where dissociative recombination is nearly absent. During solar maximum, helium ions become the dominant species at some altitudes (Heelis et al., 1990;González et al., 2004). At the magnetic equator, however, helium ions are bound at the apex point of the flux tubes. The diurnal behavior of the E×B drift forces the flux tubes to low altitudes after sunset when helium ions are removed through dissociative recombination with molecular nitrogen. Consequently, helium ion fraction peaks during the day and is limited by the amount of helium that can be produced in a single day.
In summary, the salient features of Fig. 7    to model predictions for equinox solar minimum conditions over Jicamarca made by Huba et al. (2000) except that the measured plasma temperatures are somewhat smaller at sunrise. These features have also been suggested by gated analyses in the past and so do not appear to be creations of the fullprofile analysis. However, further validation studies using in situ satellite temperature and composition measurements and numerical modeling are warranted and will be undertaken.

SAMI2 model comparison
The Naval Research Laboratory has developed a two dimensional low-to mid-latitude ionosphere model: SAMI2 (Huba et al., 2000). SAMI2 model results have been compared with ARGOS satellite data, and good agreement was found (Huba et al., 2002). SAMI2 has also been used to model the ionospheric disturbances observed at Millstone Hill driven by solar wind fluctuations (Huba et al., 2003) and a penetration electric field (Swisdak et al., 2006). SAMI2 models the plasma and chemical evolution of seven ion species (H + , He + , N + , O + , N + 2 , NO + 2 and O + 2 ) in the altitude range 85 km-20 000 km. The complete ion temperature equation is solved for three ion species (H + , He + and O + ) as well as the electron temperature equation. Ion inertia is included in the ion momentum equation for motion along the geomagnetic field. An IGRF-like magnetic field model is used. The plasma is modeled from hemisphere to hemisphere in the low-to mid-latitude ionosphere (±60 • magnetic latitude) and to ∼2000 km in the high latitude ionosphere. The neutral species are specified using the empirical NRLMSISE00 model (Picone et al., 2002) that is based on MSIS86 (Hedin, 1991) and a new horizontal wind model HWM07 (D. Drob, private communication, 2007).
The simulation results presented in this paper use a grid (n z , n f )=(201,86) where n z is the number of points along a dipole field line and n f is the number of points in altitude along the magnetic apex (i.e., number of magnetic field lines). The altitude range of the simulation study is 90-2500 km. The geophysical parameters used are the following: day 110, F10.7=80, F10.7A=80, and A p =11. The simulations were started at 00:00 UT and run for 48 h; results are presented from the final 24 h.
The comparison between the Jicamarca data and SAMI2 simulation results for the electron density, electron temperature, ion temperature, and hydrogen ion fraction is very good. Jicamarca observations during a period of low www.ann-geophys.net/26/59/2008/ Ann. Geophys., 26, 59-75, 2008 Fig. 7. Summary Jicamarca data for the full profile analysis. From top to bottom, the panels represent electron density, electron temperature, ion temperature, hydrogen ion fraction, and helium ion fraction (18% full scale). A data gap exists when the galaxy was directly over the radar. Fig. 7. Summary Jicamarca data for the full profile analysis. From top to bottom, the panels represent electron density, electron temperature, ion temperature, hydrogen ion fraction, and helium ion fraction (18% full scale). A data gap exists when the galaxy was directly over the radar. solar flux activity found helium ions comprised ∼ 15% of the electron density in the daytime F region in the altitude range 500-900 km. The helium ion layer formed shortly after sunrise at ∼500 km and rose to higher altitudes during the day. SAMI2 model results showed similar behavior for the helium ion layer just after sunset. For nominal ionospheric conditions, helium ions were ∼12% of the electron density shortly after sunrise. The helium ion layer rose during the day to altitudes above 1000 km. Simulation results using different meridional neutral wind conditions indicate that the neutral wind can substantially modify the fraction of helium ions in the equatorial F region. We intend to investigate this issue in a future paper.

Summary
Gated analysis of long-pulse incoherent scatter data is impractical at Jicamarca, where the pulses are necessarily longer than typical ionospheric density and temperature vertical length scales. We have therefore used a full-profile analysis for data inversion and parameter estimation. The method involves solving simultaneously for ionospheric profiles consistent with available data. The problem is mixed determined, and a broad family of profiles exist which give statistically indistinguishable model prediction errors. Fullprofile analysis does not create the problem but instead addresses it transparently.
Ann. Geophys., 26, 59-75, 2008 www.ann-geophys.net/26/59/2008/  Certain profiles are ruled out of the problem on the basis of prior physical knowledge (we expect T e ≥T i ) or mathematical necessity (the total ion fraction is unity). Of the remaining candidate profiles, the best candidates are the most featureless or regular ones. Regularization is consistent with Occam's Razor: to the extent the data can be reproduced without a feature in the model, that feature is not supported by the data and should not be considered.
For each regularization cost function, there is a weight or regularization parameter to be set. Their determination is arguably subjective. Here, the weights on the temperature and composition constraints are quite large, as there is little point in exploring parameter regimes that are physically or mathematically impossible. The weights on the global sec-ond derivative terms meanwhile are adjusted so as to prevent oscillatory structure in the profiles which is not observed in alternating code and double pulse experiments. In both cases, neither doubling nor halving the weights changes the results qualitatively. Validation studies may help us to refine the weighting in the future.
The full-profile analysis can be improved in a number of ways. The full error covariance matrix can be included throughout the inversion with little added computational burden using sparse-matrix routines. Noise estimation and debris removal can be further improved. Most importantly, the results would benefit from a more comprehensive statistical inversion, with both double-pulse and Faraday angle data included in the overall forward model. The knots of the Bwww.ann-geophys.net/26/59/2008/ Ann. Geophys., 26, 59-75, 2008 splines should in that case be redistributed to allow for the fine structure often evident in the profiles below 450 km. Eventually, it should be possible to compute model parameters not just for all altitudes but for all times simultaneously, exploiting temporal as well as spatial continuity in the inversion.

Appendix A
Here, we fill in the steps leading to the stated result in (2). We first express an estimate of the real part of a lag product as: where angle brackets represent the expectation and where is the statistical error associated with approximating and ensemble average with a finite sum of samples. Then RE(ρ 12 − ρ 12 )= A 12 12 , and RE(ρ 12 − ρ 12 )RE(ρ 34 − ρ 34 ) = Â 12 A 34 12 34 (A2) where the last term in the angle brackets can be determined from the relationship A 12 A 34 = A 12 A 34 (1+ 12 34 ), which holds for zero-mean Gaussian random variables. Next, define the estimator A 12 by: Consequently, we may write With the help of a theorem for even moments of Gaussian random variables (Reed, 1962) we then obtain: A 12 A 34 = A 12 A 34 + 1 2k RE ρ 13 ρ * 24 + ρ 23 ρ *

(A4)
Combining this expression with (A2) and the relationship immediately below it yields the desired result stated in (2).