Performance evaluation of the particle swarm optimization algorithm to unambiguously estimate plasma parameters from incoherent scatter radar signals

Simultaneously estimating plasma parameters of the ionosphere presents a problem for the incoherent scatter radar (ISR) technique at altitudes between ~ 130 and ~ 300 km. Different mixtures of ion concentrations and temperatures generate almost identical backscattered signals, hindering the discrimination between plasma parameters. This temperature–ion composition ambiguity problem is commonly solved either by using models of ionospheric parameters or by the addition of parameters determined from the plasma line of the radar. Some studies demonstrated that it is also possible to unambiguously estimate ISR signals with very low signal fluctuation using the most frequently used non-linear least squares (NLLS) fitting algorithm. In a previous study, the unambiguous estimation performance of the particle swarm optimization (PSO) algorithm was evaluated, outperforming the standard NLLS algorithm fitting signals with very small fluctuations. Nevertheless, this study considered a confined search range of plasma parameters assuming a priori knowledge of the behavior of the ion composition and signals with very large SNR obtained at the Arecibo Observatory, which are not commonly feasible at other ISR facilities worldwide. In the present study, we conduct Monte Carlo simulations of PSO fittings to evaluate the performance of this algorithm at different signal fluctuation levels. We also determine the effect of adding different combinations of parameters known from the plasma line, different search ranges, and internal configurations of PSO parameters. Results suggest that similar performances are obtained by PSO and NLLS algorithms, but PSO has much larger computational requirements. The PSO algorithm obtains much lower convergences when no a priori information is provided. The a priori knowledge of Ne and Te/Ti\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{e}/{T}_{i}$$\end{document} parameters shows better convergences and “correct” estimations. Also, results demonstrate that the addition of Ne\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${N}_{e}$$\end{document} and Te\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${T}_{e}$$\end{document} parameters provides the most information to solve the ambiguity problem using both optimization algorithms.


Introduction
Incoherent scatter radar (ISR) is one of the most powerful ground-based sounding techniques to estimate plasma parameters of the ionosphere at altitudes between 80 and ~ 1000 km (Evans 1969;Beynon and Williams 1978). Multiple parameters are estimated simultaneously through the fitting of the backscattered signal autocorrelation function (ACF) or its Fourier transform, the incoherent scatter spectra (ISS), at each altitude range defined by the transmitted pulse width of the radar. Nevertheless, the estimation process presents problems due to the existence of a mixture of different ion species (Oliver 1979). Different combinations of ion temperatures as well as molecular and atomic ion concentrations generate Open Access *Correspondence: mimartinezl@udec.cl 1 Physics Department, Universidad de Santiago de Chile, Av. Ecuador 3493, Santiago, Chile Full list of author information is available at the end of the article Page 2 of 25 Martínez-Ledesma and Jaramillo Montoya Earth, Planets and Space (2020) 72:172 almost identical signals, producing an ambiguous estimation of plasma parameters (Aponte et al. 2007). Indeed, molecular ions NO + and O + 2 are predominant at low altitudes (Evans and Oliver 1972;Hoffman et al. 1969), but are not distinguishable to the radar because of the similarity of their molecular masses (Oliver 1979). Therefore, these molecular ions are typically assumed to behave as a single mixed molecular species ( M + ) with a mass of 30.5 amu (i.e., a mixture of 25% O + 2 and 75% NO + ; Cabrit and Kofman 1996;Lathuillere et al. 1983). On the other hand, the topside ionosphere is dominated by atomic oxygen ions ( O + ), while at higher altitudes the protonosphere is mainly governed by lighter atomic ions ( H + and He + ; Schunk and Nagy 2009;Kelley 1989). Therefore, in the transition region between molecular and atomic oxygen ions, found at altitudes between ~ 130 to ~ 300 km, two possible solutions of plasma parameters can be obtained (Lathuillere et al. 1983;Wu et al. 2015). Different electron and ion temperatures are acquired depending on the relative abundance of molecular ions value ( p = n(M + )/N e ). This effect is known as the "temperature-ion composition ambiguity" (TICA) problem of ISR (Martínez-Ledesma and Díaz 2019).
Moreover, the multiple couplings of the magnetosphere-ionosphere-thermosphere (MIT) system (Sarris 2019; Borovsky and Valdivia 2018) and internal ionospheric processes provide a large variability to the ion composition at the ionospheric valley region. A study conducted by Lathuillere and Pibaret (1992) showed a large diurnal and seasonal variability of the transition altitude between molecular and atomic oxygen ions ( z 50 ). Additionally, this transition altitude revealed a clear dependence on the Kp geomagnetic activity index. Similar ion composition variabilities were obtained by Shibata et al. (2000) in the auroral zone. At high latitudes, magnetospheric activity generates ion outflows, carrying fluxes of atomic oxygen ions from the ionosphere to the magnetospheric plasma sheet (Borovsky and Valdivia 2018). In these outflow events, molecular ions are upwelled to high ionospheric altitudes during magnetic storms. These molecular ions may be transported to mid-latitudes during strong storm conditions. At mid-low latitudes, the ion downwelling might carry atomic ion species to lower ionospheric altitudes (Richmond and Lu 2000). Satellite measurements in the auroral zone have detected atomic oxygen ion depletions and increments of molecular ions related to ion outflow events (Brinton et al. 1978). Even during low magnetic activity and solar minimum conditions, simple ionospheric models show a large variability of ion composition at high latitudes (Sojka et al. 1981). These models agree with the ion densities obtained by the AE-C satellite located at ~ 300 km altitude (Brinton et al. 1978). Furthermore, recent satellite measurements have observed molecular ions in the magnetospheric ring current during small magnetic storms, indicating that the ion outflow is in fact a highly frequent phenomenon (Seki et al. 2019). Moreover, molecular ions were not found in the magnetosphere during geomagnetic quiet periods, suggesting that magnetic storms are a key driver of the low-altitude ionospheric ion loss through ion outflow processes. To date, it is still unknown what mechanisms control the ion outflow, as well as its spatial and temporal scales, and the amount of flux transported through this process (Borovsky and Valdivia 2018). The correct determination of ion composition and temperatures by the ISR method could provide an extremely relevant insight into these questions, helping disentangle several of the most outstanding current unknowns about the Earth's MIT system (Sarris 2019).
Different methods have been studied in the literature to solve the TICA problem. The most common solution is the use of theoretical models of plasma parameters to restrict the solutions of the ambiguous estimation. Initial methods considered equal ion and neutral temperature profiles or smooth variations of plasma parameters with altitude (Waldteufel 1971;Evans and Oliver 1972). Parametric models of ion composition and temperature were implemented by Oliver (1979) based on early rocket and satellite measurements. Models of ion composition were also employed to simultaneously estimate parameters of the entire ionospheric column using the full profile inversion method (Cabrit and Kofman 1996;Litvine et al. 1998). Furthermore, in the auroral region, different models were developed to analyze Joule heating events (Kelly and Wickwar 1981) and strong electric fields (Blelly et al. 2010;Zettergren et al. 2011). Alternatively, several methods have been developed to provide additional information from the plasma line of the radar, constraining the feasible solutions to solve the ambiguity (Wand 1970;Bjørnå and Kirkwood 1988;Fredriksen et al. 1989;Fredriksen 1990;Nicolls et al. 2006;Aponte et al. 2007). As the ion acoustic band of the radar presents the most pronounced features of the spectrum, plasma parameters have been commonly estimated from this frequency band (Akbari et al. 2017). Alternatively, the plasma line is a frequency peak in the megahertz range, which depends on electron density, electron temperature, and the magnetic aspect angle (Yngvesson and Perkins 1968). The plasma line is not routinely detected by most ISRs, but it is occasionally enhanced by nonthermal effects such as local or magnetic conjugate photoelectrons, secondary electrons, and electron beams (Akbari et al. 2017). Plasma parameters that can be estimated from the plasma line resonance frequency are the electron density ( N e ) and the electronto-ion temperature ratio ( T e /T i ;Wand 1970;Waldteufel 1971;Aponte et al. 2007), the electron temperature ( T e ;  Kofman et al. 1981;Bjørnå and Kirkwood 1988;Nicolls et al. 2006), and the electron bulk velocity ( V e , Behnke and Ganguly 1986;Akbari et al. 2017). Moreover, Oliver (1979) proposed a different method to solve the ambiguity based on the use of ISR signals with very small fluctuations. This study determined the two possible parameters solutions and posteriorly selected the solution that provided the smoother ion composition profile at different altitudes. Later, this method was implemented in Lathuillere et al. (1983) and Lathuillere and Pibaret (1992) by post-integrating signals with long integration times of 5 min.
The most recent study of the TICA problem (Martínez-Ledesma and Díaz 2019) presented a theoretical framework to quantify the probability of solving the TICA problem and determined the signal fluctuation thresholds for which the problem is unambiguously solved. This study evaluated the unambiguous estimation performance of the most commonly used non-linear least squares (NLLS) optimization algorithm in ISR analyses, for both range-gate and full profile estimation methods (Erickson 1998;Swoboda et al. 2017;Hysell et al. 2008;Nikoukar et al. 2008): the Levenberg-Marquardt (L-M) algorithm (Levenberg 1944;Marquardt 1963). In their study, the authors assessed the convergence of the estimated solutions by using the statistical distribution of the reduced Chi-square ( χ 2 r ) minimization cost function. By selecting only convergent solutions, the TICA problem was solved at a very small signal fluctuation level. The addition of a priori information from the plasma line increased the fluctuation level at which the signal could be unambiguously estimated. Furthermore, the study of Martínez-Ledesma and Díaz (2019) demonstrated that the optimization search was not always successful, even when the initial parameters of the search were near the correct solution. However, it was determined that the success rate depends on both the distance of the initial search parameters to the true solution and the signal fluctuation level. To overcome the initial assumptions constraint, the authors suggested a new estimation method based on Monte Carlo simulations and L-M optimization: the Monte Carlo Levenberg-Marquardt (MCLM) method. In this new method, the analyzed fittings used different initial parameters uniformly selected from the search space of the parameters. Hence, the most frequently obtained set of parameters is most likely to be the solution of the estimation problem. In this study, the probabilities of convergence and unambiguous estimation of the MCLM method were calculated for different combinations of parameters assumed to be known a priori with different levels of uncertainty.
A different fitting algorithm was evaluated by Wu et al. (2015) to unambiguously estimate plasma parameters of ISR signals obtained at the Arecibo Radio Observatory: the particle swarm optimization (PSO) algorithm . This study analyzed the performance of the PSO algorithm by fitting ISR signals with different combinations of parameters assumed to be known a priori. Their results suggested that the PSO algorithm outperformed the standard NLLS optimization algorithm commonly used in ISR fittings. Nevertheless, three elements from this work are important to note and might be interesting to explore: (i) the signals analyzed were obtained with signal-to-noise ratios (SNR) between 15 and 20, owing to the simultaneous frequency transmission (Sulzer 1986a) and coded long pulse (Sulzer 1986b) techniques used at Arecibo. This low signal fluctuation characteristic is difficult to achieve by other ISR facilities other than the Arecibo Radio Observatory. (ii) The PSO algorithm presents the problems of stagnation and early convergence to local minimums, which were not evaluated or analyzed in the study. And (iii) the PSO implementation constrained the search range of plasma parameters using theoretical assumptions of the ion composition that may be inappropriate in certain cases.
The work of Wu et al. (2015) proposed an interesting and powerful tool for parameter estimation of ISR signals. The PSO is a widely used optimization algorithm capable of determining optimal solutions in multi-dimensional complex search spaces (Freitas et al. 2020;Zhang et al. 2015). However, no previous study has determined the capacities of the PSO algorithm to solve the TICA problem in different noise regimes. The objective of the present study is to analyze the unambiguous estimation performance of the standard PSO algorithm at different signal fluctuation levels. Besides, our study provides new techniques to detect stagnation and determine incorrect estimation cases of PSO fittings. Furthermore, a comparison study of the PSO algorithm is done with both limited and full ion composition search ranges. Additionally, the performance of the standard PSO algorithm is compared to the MCLM method to determine which optimization technique provides the best estimation characteristics to solve the TICA problem.
This article is organized in four sections. "Methods" section describes the PSO algorithm implemented, the simulation methods, and the statistical analysis of simulation results. "Results and discussion" section shows the results and discussion, and "Conclusions" section indicates our conclusions.

The particle swarm optimization algorithm
The PSO algorithm was originally proposed and developed in 1995 by R.C. Eberhart and J. Kennedy in two different papers Eberhart Page 4 of 25 Martínez-Ledesma and Jaramillo Montoya Earth, Planets and Space (2020) 72:172 and . The PSO algorithm is a method for optimizing nonlinear functions based on social behaviors found in nature, such as bird flocks seeking for food. It is a stochastic algorithm that resembles the collective behavior of a decentralized and self-organized system (i.e., swarm intelligence; Zhang et al. 2015). Furthermore, this algorithm has been proven to work efficiently and stably with noise (Wang et al. 2018). This method is based on the premise that the knowledge is spread between individuals of a population, but also relies on the social distribution of information among different generations of this population (Freitas et al. 2020). Multiple variants, enhancements, and extensions that improve the original PSO algorithm have been proposed in the literature to address different kinds of problems (for reviews of this technique, see Freitas et al. 2020;Wang et al. 2018;Sengupta et al. 2019;Zhang et al. 2015).
In the PSO algorithm multiple individual agents, called particles, are potential solutions of the optimization problem (Wang et al. 2018). The positions and velocities of the n particles that comprise the swarm are initialized randomly using a uniform distribution in the search space (which prevents the PSO algorithm being sensitive to initial conditions). Then, supposing a P-dimensional space of parameters of the objective function, the particles move iteratively through the solutions space to find the global minimum of the problem by adjusting their trajectories and velocities. At each iteration, the algorithm stores the best fit of each particle ( pbest i ) as well as the current global best fit of all particles ( gbest ). The optimization ends when the number of iterations of the algorithm reaches a maximum value or when the global best fit value reaches a specific tolerance. In this case, the optimization result corresponds to the global best location of all particles.
The way in which the PSO algorithm has been designed allows the particles to explore and exploit the search space with the aim of finding the global optimum (Trelea 2003). At each iteration, the algorithm calculates the positions and velocities of the particles by sharing information between them. Each new position is a combination of the individual knowledge (the previous velocity and position of the particle and its best historical position; pbest i ) and the collective knowledge of the swarm ( gbest ; Yang 2014). In addition, the movement of each particle contains a stochastic component that contributes to maintaining an adequate level of diversity in the swarm population (Freitas et al. 2020). Therefore, the velocity ( v ij ) and position ( x ij ) of each particle are updated according to the following equations (Shi and Eberhart 1998): where v ij is the velocity vector of particle i into the j dimension, x ij is the position of particle i into the j dimension, ϕ 1 and ϕ 2 are uniformly distributed random numbers ( ϕ 1 , ϕ 2 ∈ [0, 1] ), ω is the inertia weight, and c 1 and c 2 are the cognitive and social components, respectively.
Nevertheless, the PSO algorithm frequently selects local optimum solutions in mid-and large-size problems (Lin et al. 2014). This selection of local optimums depends on the complexity of the objective function, but it can also be related to the fast disappearance of the diversity of the particles in the search space (Wang et al. 2018). Furthermore, the PSO algorithm presents the problem of premature convergence (Wang et al. 2018). Previous studies have verified the convergence uniqueness of the PSO algorithm, demonstrating that the particles converge in mean square to the best position found by the swarm (Wang and Shen 2012;Jian et al. 2007). Nevertheless, this convergence position is not obliged to be a global or local optimum of the objective function (Jian et al. 2007;Yuan and Yin 2015). This early convergence effect is known as a "stagnation problem" and it occurs when particles stop changing their positions, generating an early convergence to a position not guaranteed to be the global or local minimum (Freitas et al. 2020). Additionally, it has been theoretically and experimentally demonstrated that the PSO algorithm would not converge to the optimum solution if the number of used particles is very small for the problem (Raß et al. 2015).

Relevance of the PSO parameters selection
The selection of the parameters c 1 , c 2 , and ω of the PSO algorithm has a large impact on the performance of the optimization (Zhang et al. 2015). In multimodal problems, where multiple areas of the search space may contain local solutions, the determination of those parameters is even more relevant (Freitas et al. 2020), and it is difficult to improve the accuracy of the algorithm by their adjustment (Lin et al. 2014). The cognitive component ( c 1 ) represents the individual knowledge of each particle and encourages the particle to move towards its own best-known position. Alternatively, the social or cooperation component ( c 2 ) represents the collaborative knowledge of the swarm and moves the particle toward the global solution. Both parameters are positive constants ( c 1 > 0 and c 2 > 0 ), also called "acceleration constants" or "learning factors" (Wang et al. 2018;Zhang et al. 2015). Commonly used cognitive and social components are those originally suggested by Kennedy and Eberhart (1995), i.e., c 1 = c 2 = 2 (Wang et al. 2018; Page 5 of 25 Martínez-Ledesma and Jaramillo Montoya Earth, Planets and Space (2020) 72:172 Yang 2014). Nevertheless, Van den Bergh and Engelbrecht (2006) demonstrated that the critical condition c 1 + c 2 < 2(ω + 1) is required to guarantee convergence. Furthermore,  suggested that the sum of these constants should not exceed 3 (i.e., c 1 + c 2 ≤ 3 ), and also recommended a larger social component ( c 2 > c 1 ) to provide better performance. Alternatively, the inertia weight ( ω ) parameter was first introduced by Shi and Eberhart (1998) to balance the search of the particles between the local and global environments (i.e., exploration and exploitation of the search space, respectively; Zhang et al. 2015). The use of this parameter highly improved the optimization performance of the original PSO algorithm (Wang et al. 2018). Therefore, the use of ω became extensive by the community, being commonly referred to as the standard PSO (SPSO; Freitas et al. 2020). A high value of ω provides higher relevance to the individual knowledge of the particle. Alternatively, a small value of ω increases the use of the collective swarm knowledge, providing the necessary momentum for particles to travel across the search space (Zhang et al. 2015), and preventing the algorithm from converging to local minimums (Freitas et al. 2020). As this parameter has the greatest influence on the algorithm performance, multiple studies have been focused on evaluating different methods to configure the inertia weight value (e.g., see Wang et al. 2018). Due to the difficulty of dynamically adapting the value of ω to every particular search space, Shi and Eberhart (1998) proposed the idea of a time-varying inertia. In a later study, Shi and Eberhart (1999) demonstrated that a linearly decreasing ω that starts at 0.9 ( ω max ) and ends at 0.4 ( ω min ) with c 1 = c 2 = 2 provided fast convergences to multiple experimental tests, independently of the number of particles used. More recently,  suggested that ω should decrease linearly from 0.8 to 0.5 to guarantee the critical condition obtained by Van den Bergh and Engelbrecht (2006).
The study of Wu et al. (2015) implemented a SPSO algorithm to estimate plasma parameters from ISR spectrum signals. Their parameter configuration was initially suggested by Shi and Eberhart (1999). In a previous study, Chen et al. (2013) demonstrated that this parameter configuration provides better results than the configuration proposed by Eberhart and Shi (2000), which uses a fixed inertia weight ω = 0.729 and c 1 = c 2 = 1.49445.
Alternatively, the selection of the number of particles of the swarm population ( n ) depends on the objective function to optimize. Although the algorithm performance is not very sensitive to the size of the swarm (Wang et al. 2018), an insufficient number of particles may not converge to the optimum solution (Raß et al. 2015). In the study of Wu et al. (2015), a total of 100 particles were used to estimate more than one plasma parameter, but only 10 particles were used when fitting one parameter. Furthermore, complex problems require not only increasing the number of particles, but also increasing the maximum number of iterations. A maximum of 5000 iterations was used by Wu et al. (2015) to solve the problem of fitting one plasma parameter, but much larger maximum values were used for fitting two and three parameters (20,000 and 30,000 iterations, respectively). Nevertheless, it is assumed that the PSO algorithm implementation of Wu et al. (2015) increases the iteration counter each time a particle is evaluated. Therefore, configurations of a maximum number of 500, 200, and 300 iterations of all particles were used in their study when evaluating 1, 2, and 3 plasma parameters, respectively.

PSO parameters choice
To allow a direct comparison with the results obtained by Wu et al. (2015), in our work, we used the parameter configuration proposed by Shi and Eberhart (1999), i.e., c 1 = c 2 = 2 and a linearly decreasing ω from 0.9 to 0.4. From herein, this configuration is known as "Param. 1". Besides, to evaluate the performance of PSO using a different configuration, a second parameter configuration was implemented ("Param 2"): c 1 = 1.2 , c 2 = 1.8 , as well as the linearly decreasing ω from 0.8 to 0.5 suggested by . It is important to note that the later parameters c 1 and c 2 verify the conditions c 2 > c 1 and c 1 + c 2 ≤ 3 recommended by .
In our work, more than one plasma parameter was estimated simultaneously. Therefore, we selected n = 100 , the number of particles used by Wu et al. (2015) when a similar number of plasma parameters were estimated. Furthermore, to verify the sensitivity of the optimization performance to this configuration, several fitting simulations were done with different numbers of particles. Results of these simulations are explained in Additional file 1: Text S1 and Figures S1, S2, and S3.
To properly determine the maximum number of iterations required to estimate plasma parameters from ISR signals (i.e., max_iterations ) and reduce the overall computation time of all Monte Carlo simulations, different fitting simulations were done with max_iterations values of 1000, 500, and 300. Results of these simulations are shown in Additional file 1: Text S2 and Figures S4, S5, S6, and S7. These simulation results suggest that simulations with max_iterations = 500 provide the best tradeoff between computing time and estimation performance. Note that, in our PSO implementation, the iteration counter increases each time all particles are evaluated.

Optimization algorithm implementation
Despite the fact that there are multiple possible variants of the PSO algorithm that could be considered to improve different aspects of the standard algorithm, our study is mainly focused on evaluating the performance of the most common PSO implementation to solve the TICA problem, as in Wu et al. (2015). In this way, to allow a direct comparison with the study of Wu et al. (2015), we have implemented a SPSO algorithm to fit simulated spectra (ISS) with known ionospheric parameters. The flowchart of the implemented algorithm is illustrated in Fig. 1. Different modifications of the PSO algorithm and hybridizations with other optimization algorithms (such as those indicated in Freitas et al. 2020;Wang et al. 2018;Zhang et al. 2015) should be considered in future studies to evaluate possible improvements to this technique, but the discussion of those other algorithms is beyond the scope of this paper. The main difference between our implementation and the SPSO algorithm employed by Wu et al. (2015) is the use of a counter to identify stagnation (indicated by a dashed box in the flowchart in Fig. 1). The stagnation problem is generated when the best historical position of particles and the global best position remain constant during some iterations. In these cases, particles would only be able to leave their current positions if their previous velocities were non-zero (Freitas et al. 2020). In our implementation, a counter named count_const increases when the global best fit value remains equal during successive iterations; otherwise it is set to zero. If this counter reaches the maximum value ( max_count_const ) the search is stopped. This stop condition reduces the computing time when stagnation occurs, but does not affect the estimation performance of the algorithm. It is relevant to highlight that, when all particles were initially located far from the local or the global optimum, the global best fit would remain constant during the initial part of the search. Therefore, erroneous stagnation detections may occur during the initial exploration of the search space. To avoid possible stagnation misinterpretations during the initial iterations of the algorithm, the stagnation counter is allowed to increase only when the minimum number of iterations is reached ( count_const_init).
To detect the stagnation of the optimization search, the maximum number of consecutive repetitions with identical best fit ( max_count_const ) was set to 200. This value was selected to ensure that the number of repetitions should be less than half the number of maximum iterations of the fitting ( max_iterations ). To avoid erroneous initial stagnation detections, the count_const_init parameter was configured to 100. Furthermore, simulations were performed using the following values: count_const_init of 100 and 200, and their results are shown in Additional file 1: Text S2 and Figures S4, S5, S6, and S7.
By the addition of a penalty function, the PSO algorithm can be applied to problems with restricted search spaces (i.e., constrained optimization problems). To implement this feature, particles with positions out of the allowed space of parameters ( x ij < x j,min or x ij > x j,max ) should have their objective functions increased by a large penalty factor (Freitas et al. 2020). As a result of this addition, the minimization process of the new objective function bounds the search space by always selecting solutions without penalization. In our work, a penalty function has been used to constrain the search to the allowed plasma parameters ranges, as in the work of Wu et al. (2015). The ranges of plasma parameters considered in our study are indicated in Sect. 2.5. Additionally, to improve the optimization performance, the velocity of the particles can be restricted ( v j,min < v ij < v j,max ). This restriction helps avoid an overfly of potential solutions or an insufficient exploration of the search space (Freitas et al. 2020;Clerc and Kennedy 2002;Eberhart and Shi 2000). In our study, no velocity limitation has been implemented, as in Wu et al. (2015).
At each iteration, the PSO algorithm requires to simultaneously evaluate the objective function of multiple particles. To reduce the computing time required, the ISR theoretical spectrum of Martínez-Ledesma and Díaz (2019) was implemented using vectorization. Despite the fact that the PSO algorithm is known to be a parallelizable algorithm (Wang et al. 2018), the ISR spectrum model employed is not easily parallelizable due to the lack of parallelization of the Faddeeva Dawson function (Johnson 2012) it uses. Using this vectorizing computing scheme, multiple sets of plasma parameters were simultaneously evaluated. This process finally obtained a matrix of theoretical spectrums that can be evaluated to minimize the cost function of multiple sets of plasma parameters. Computing times obtained with this vectorized scheme were linearly increased with the input vector length, as it can be verified in Additional file 1: Figure  S3. Additionally, to improve the overall computing time of Monte Carlo simulations, the fittings of different spectra were parallelized by executing each optimization in a different CPU core at the nodes of a supercomputing infrastructure (of the Chilean National Laboratory for High Performance Computing, NLHPC), as in Martínez-Ledesma and Díaz (2019).

Monte Carlo data fitting simulations
Simulation methods were identical to those implemented in Martínez-Ledesma and Díaz (2019). The simulations were done considering a radar frequency of 450 MHz  . 1 Flowchart of the standard PSO algorithm implemented. n : swarm population size; P : number of plasma parameters; max_iterations : maximum number of generations; t : generations counter (from 1 to max_iterations ); i : particle's id counter (from 1 to n ); j : parameter's dimension from 1 to P ; v ij (t) and x ij (t) : i th particle's j th dimension's t th generation's velocity and position value, respectively; count_const: counter of stagnation cases (i.e., identical fit values of gbest ; from 1 to max_count_const ); count_const_init : minimum generation to start counting for stagnation cases and without magnetic aspect angle effects . Synthetic radar signals ( m ) were created using a theoretical ISS model ( f (x) ) based on Kudeki and Milla (2011) plus an additional noise contribution ( ε ; Vallinkoski 1988). These radar signals were assumed to be obtained with weak SNR conditions, obtaining a signal-independent noise related to thermal and sky noises (Huuskonen and Lehtinen 1996;Sulzer 1986a). Signals from multiple radar pulses are required to be integrated to provide an accurate statistical estimate (Farley 1969). According to the central limit theorem, these multiple integrations provide Gaussian characteristics to the random noise (Vallinkoski 1988). Therefore, in our simulations, an additive white Gaussian noise (AWGN) was used with zero mean and a standard deviation ( σ ) defined by the selected signal fluctuation level ( δ ) multiplied by the maximum absolute value of the theoretical spectrum signal.
As a consequence of the noise normalization with respect to the maximum of the power spectrum, the same signal fluctuation level may correspond to different SNR values. As the shape of the ISR signal spectrum depends on the plasma parameters, the maximum of the power spectrum is dependent on those parameters. Nevertheless, the average signal fluctuation level should be correctly related to the SNR level of ISR signals, as explained in Additional file 1: Text S3 and Figure S8.
Plasma parameters were estimated from ISR signals using the PSO optimization algorithm by the minimization of the reduced Chi-square cost function ( χ 2 r ). This objective function is a normalization of the Chi-square least squares estimator ( χ 2 ) commonly used in ISR fittings (Erickson 1998;Vallinkoski 1988) by the number of degrees of freedom ( DoF ). Therefore, this objective function is given by (Bevington and Robinson 2003;Taylor 1997): where m i , σ 2 i , and f i (x) are the M vector components of the measured signal, measured signal variance, and theoretical model function, respectively. In our work, the ionacoustic ISR spectra are obtained at frequencies between ± 10 kHz and with a vector length M=50.
The reduced Chi-square criterion is the maximum likelihood estimator for signals with AWGN stochastic noise characteristics and a diagonal measurement variance-covariance matrix (Swoboda et al. 2017;Vallinkoski 1988). The DoF of the estimation problem was calculated as the vector length of the ISR signal minus the number of parameters to estimate (i.e., DoF = M − P ). The χ 2 r criteria allows the 'goodness of the fit' to be statistically (Bevington and Robinson 2003), providing information about the solutions that converged to a minimum of the cost function.
Estimated plasma parameters were different than those selected in the PSO analysis of Wu et al. (2015). In our work, plasma parameters were electron density ( N e ), electron temperature ( T e ), ion temperature ( T i ), and molecular ion concentration ( p ). The ion drift radial velocity ( V i ) was assumed to be known a priori, as it does not affect the determination of the other parameters (Martínez-Ledesma and Díaz 2019; Wu et al. 2015). On the other hand, the study of Wu et al. (2015) assumed that N e and V i parameters were obtained independently. Also, instead of estimating the T e parameter, these authors estimated the ratio of electron to ion temperature ( T e /T i ). Therefore, they considered the following list of fitting parameters: T i , T e /T i , and p . Furthermore, in our work, different combinations of plasma parameters were assumed to be known a priori from the plasma line (given N e , N e and T e /T i , and N e and T e parameters), based on previous studies of the plasma line of the radar (Aponte et al. 2007;Nicolls et al. 2006;Bjørnå and Kirkwood 1988;Waldteufel 1971). Also, this selection of plasma parameters allows a direct comparison with the estimation results of the MCLM method (Martínez-Ledesma and Díaz 2019). Alternatively, Wu et al. (2015) studied different scenarios in which the ion composition or the ion temperature were assumed to be known a priori, although these parameters are currently not possible to be obtained independently from the plasma line. Also, note that in Wu et al. (2015) all a priori known plasma parameters were considered without uncertainty. Although the addition of a priori parameters without uncertainty is unrealistic, it provides the best-case estimate that can be obtained by the optimization algorithm, as demonstrated by Martínez-Ledesma and Díaz (2019).
To determine the unambiguous estimation performance of the PSO algorithm, Monte Carlo simulations of the plasma parameters estimation process were conducted. Simulations were done at different signal fluctuation levels to determine the impact of increasing the noise level of ISR signals on the estimation performance. To allow a direct comparison with the results of the MCLM method, these fittings were done for the same 2000 true input parameters ( x true ) used in the Monte Carlo analysis done by Martínez-Ledesma and Díaz (2019). Furthermore, to provide enough statistical representation, the fitting process for each set of input parameters was repeated 100 times adding a different random noise. The same set of true input plasma parameters was used in all simulations, allowing a direct comparison of results for different noise fluctuation levels. These true input parameters were Page 9 of 25 Martínez-Ledesma and Jaramillo Montoya Earth, Planets and Space (2020) 72:172 uniformly selected from the space of parameters defined by the following ranges (Martínez-Ledesma and Díaz 2019): electron density, 10 9 ≤ N e ≤ 10 12 m −3 , electron temperature, 300 ≤ T e ≤ 5000 °K, ion temperature, 300 ≤ T i ≤ 3000 K, electron-to-ion temperature ratio, 0.1 ≤ T e /T i ≤ 5 , and ion composition, 0 ≤ p ≤ 1.
On the other hand, in our work, the plasma parameter search intervals of the optimization algorithm were identical to the previously indicated true input parameter ranges. However, electron and ion temperatures were configured more broadly than the input parameter ranges, with values ranging from 200 to 6000 °K, as in Martínez-Ledesma and Díaz (2019). Alternatively, in the study of Wu et al. (2015) the plasma parameter search ranges of the PSO algorithm were confined with a priori information of the true input plasma parameters. The search interval for the molecular ion fraction was limited by Wu et al. (2015) to [ max(p true − 0.5, 0) , min(p true + 0.5, 1) ], where p true is the true input ion composition parameter of the simulation. This limited ion composition search range is shown in Fig. 2. Moreover, the search ranges for T i and T e /T i were confined by Wu et al. (2015) to [100, 2T i,true ] and [1, 3], respectively.

Statistical analysis of results
To determine the impact of the PSO algorithm on the estimation ambiguity in different noise scenarios, results from the Monte Carlo simulations were analyzed. Simulations performed with different signal fluctuation levels were evaluated with the statistical methods defined by Martínez-Ledesma and Díaz (2019), which are briefly explained below.
First, the convergence of the fitting was determined by selecting the solutions with χ 2 r value smaller or equal to the maximum cost function value ( χ 2 r,max ). This maximum cost function value was selected to have a 4σ probability criterion (i.e., P χ 2 r > χ 2 r,max = 0.00317% , provided by the Chi-square statistical distribution of the χ 2 r cost function), ensuring that a fit with χ 2 r > χ 2 r,max is obtained by an invalid convergence of the algorithm. As in Martínez-Ledesma and Díaz (2019), maximum and minimum values of χ 2 r,max are 2.073 and 2.0317 for 5 and 2 plasma parameters ( P ), respectively. Figures S9 and S10 of Additional file 1 show the histograms of the estimation error of all plasma parameters, providing a graphical verification that the fittings with χ 2 r > χ 2 r,max obtain large errors in all plasma parameters.
Secondly, convergent solutions were selected as "correct" or "incorrect" depending on the distance of the estimated solution ( x ) to the true input plasma parameters ( x true ). The ion composition parameter was selected as a discriminator of correctness because of the known existence of two different solutions in the limited range [0, 1] (Lathuillere et al. 1983;Oliver 1979;Wu et al. 2015). The statistical distributions of both "correct" and "incorrect" ion composition solutions are assumed Gaussian. Therefore, the distribution of solutions follows a bimodal Gaussian mixture model (GMM) given by the following: where α is the weight of the mixture distribution ( 0 ≤ α ≤ 1 ), and N x|µ 0 , σ 2 0 and N x|µ 1 , σ 2 1 are the probability density functions (PDF) of "correct" and "incorrect" results, respectively.
For each set of true input plasma parameters, the statistical distribution of the ion composition estimation error ( ε p = p true − p ) was analyzed using the expectation maximization (EM) algorithm (Dempster et al. 1977). The determination of the GMM distribution by the EM algorithm allows the clustering of both "correct" and "incorrect" solutions.
(3) Finally, the performance of the optimization algorithm was evaluated using different probabilities. The probability of valid convergence of the fit ( P fit valid ) was calculated as: where N fit valid is the number of convergent results of a simulation ( χ 2 r ( x) ≤ χ 2 r,max ), and N total is the total number of parameters simulated ( N total = N MC · N rep , and where N MC is the number of different input parameters of the Monte Carlo simulation, N MC =2000, and N rep is the number of fitting repetitions of each input parameter of the simulation, N rep =100).
The probability of 'correct' estimation ( P correct ) was calculated as: where N correct is the number of "correct" results of a simulation obtained by the EM clustering algorithm.
And, the probability of valid convergence and 'correct' estimation ( P fit valid & correct ) was calculated as:

Results and discussion
The following sections show the results of the Monte Carlo simulations of fittings using the PSO algorithm with different combinations of plasma parameters assumed to be known a priori from the plasma line. The following case studies were analyzed: a) 4 unknown parameters ( N e , T i , T e , and p ); b) 3 unknown parameters ( T i , T e , and p ) given a priori N e ; c) 2 unknown parameters ( T i and p ) given a priori N e and T e /T i ; and d) 2 unknown parameters ( T i and p ) given a priori N e and T e .

Effects of the limited ion composition search space
In Wu et al. (2015), the limited ion composition search range (shown in Fig. 2) restricts the possible ambiguous solutions obtained by the optimization algorithm. This search interval uses a priori information of the true input ion composition, which is information not available in real radar measurements. Therefore, it is assumed that the restricted search range applied by Wu et al. (2015) artificially enhances the unambiguous estimation performance of the PSO algorithm.
To identify the number of "incorrect" solutions not considered in the limited ion composition search interval of Wu et al (2015), Monte Carlo simulations were done with the full ion composition search range (i.e., 0 ≤ p ≤ 1 ) at different signal fluctuation levels. Furthermore, these simulations were done with the PSO (4) P fit valid = N fit valid N total , parameters used by Wu et al (2015) (i.e., parameters configuration "Param. 1" indicated in "PSO parameters choice" section). Figure 3 shows the estimated ion composition parameters obtained from simulations with the full ion composition search range at different signal fluctuations and (as black dashed lines) the boundaries of the limited search interval of Wu et al (2015).
Simulation results without a priori information (case study a in Fig. 3) at mid and large signal fluctuations ( δ > 0.05% ) obtained several ion composition estimates outside of the limited search interval. The addition of information from the plasma line (case studies b, c, and d) reduced the number of estimates located in these intervals. Particularly, the addition of a priori knowledge of N e and T e information (case study d) had almost no "incorrect" estimates. Therefore, in this case study, almost all solutions were found inside of the limited search range of Wu et al. (2015). This combination of a priori known parameters solves the TICA problem even at high signal fluctuations, as does the MCLM algorithm (Martínez-Ledesma and Díaz 2019). The number of convergent solutions found outside of the search interval of Wu et al (2015) is shown in Fig. 4 at different signal fluctuations for these different combinations of parameters known a priori. At large signal fluctuations ( δ > 10% ), the number of "incorrect" solutions found in the limited search range greatly increases in all case studies. It is relevant to highlight that the area outside of the search space of Wu et al. (2015) corresponds to 25% of the total available space. Therefore, the number of cases higher than 25%, obtained at high signal fluctuations, indicates that "incorrect" ion composition values were more often obtained at highly noisy scenarios. This result suggests that ion composition solutions were not uniformly distributed in the search space at high signal fluctuations. Alternatively, the number of cases outside the search interval was abruptly reduced at small signal fluctuations in case studies a and b (with δ < 0.2% and δ < 0.5% , respectively). It is assumed that this reduction is related to the decrease of convergent estimates that occurs at small signal fluctuations, which is explained below.
Additionally, to evaluate the impact of the limited search range on the unambiguous estimation performance of the PSO algorithm, similar Monte Carlo simulations were done using the ion composition search interval employed by Wu et al (2015). Figure 5 shows the probabilities obtained with a limited ion composition search range and with a full search range (dashed and continuous lines, respectively). This figure shows the probability of convergence of the optimization algorithm ( P fit valid ), the probability of having a "correct" estimation ( P correct ), and their joint probability ( P fit valid & correct ).
Page 11 of 25 Martínez-Ledesma and Jaramillo Montoya Earth, Planets and Space (2020) 72:172 Simulations with a limited ion composition search range showed slightly higher values of P fit valid for case studies a and b (without a priori information and knowing a priori N e parameter, respectively). The increase in convergent solutions of case study b was larger than those of case study a (with values ~ 11% and ~ 4%, respectively). Also, the probability of convergence in these two case studies decayed for signal fluctuations smaller than 0.5%. Alternatively, simulations of case studies c and d obtained almost perfect convergences ( P fit valid ≈ 100% ) in all noise scenarios, independently of the search range configuration. These results suggest that the limited ion composition search range of Wu et al. (2015) does not generate a significant impact on the convergence of the PSO algorithm. Nevertheless, a relevant result was obtained in simulations without adding a priori information. In these simulations, independently of the ion composition search range, more than half of the total number of estimations were not convergent ( P fit valid < 50% ) for small signal fluctuations ( δ < 0.1% ). Moreover, at very small noise scenarios ( δ = 0.01% ) almost no estimations were convergent ( P fit valid ≈ 10% ), indicating that the optimization algorithm was unable to find a minimum of the cost function in the majority of the cases. This result Fig. 3 Scatter plot of estimated and input values of ion composition, obtained from the Monte Carlo analysis of different combinations of parameters known a priori from the plasma line using the PSO algorithm: a without a priori information, b given N e , c given N e and T e /T i , and d given N e and T e . Each color represents results obtained by simulations with a particular signal fluctuation percentage ( δ(%) ). Simulations were done with the full ion composition search range. Black dashed lines represent the maximum and minimum values of the confined ion composition search range of Wu et al. (2015) Page 12 of 25 Martínez-Ledesma and Jaramillo Montoya Earth, Planets and Space (2020) 72:172 suggests that the selection of the PSO algorithm parameters ( c 1 , c 2 , ω , and number of particles n ) made by Wu et al. (2015) is not useful to estimate plasma parameters of ISR signals in very small noise scenarios without the addition of a priori information. Note that this case study (i.e., fitting 4 plasma parameters without a priori information) was not evaluated in their study. On the other hand, the use of the limited search range produced an impact on the probability of unambiguous estimation. The increments of P correct when using the search range of Wu et al. (2015) were ~ 17% and ~ 11% for case studies a and b, respectively. Note that the effect of using the limited ion composition search range in simulations without adding any a priori information generates a P correct similar to the case of knowing a priori N e . Alternatively, case study c displayed a probability increment smaller than 5% at signal fluctuations smaller than δ ≤ 10% . On the other hand, almost identical and perfectly unambiguous estimations were obtained when knowing a priori N e and T e parameters (case study d), as in the case of using the MCLM algorithm (Martínez-Ledesma and Díaz 2019). The small increase obtained in case studies c and d implies that the knowledge of N e and T e /T i or N e and T e provided more information than the use of the confined ion composition search range of Wu et al. (2015). Also, of relevance is that P correct value without a priori information (case study a) did not reach 100% at very small signal fluctuations ( δ = 0.01% ). This result is assumed to be related to the small number of convergent solutions at this signal fluctuation level ( P fit valid ≈ 10% ) which impacted the average number of "correct" solutions.
In summary, these results verify that the use of the confined search range provides a priori information that artificially increases the unambiguous estimation performance of the algorithm in almost all case studies. Nevertheless, this segmented search range did not impact the convergence of the algorithm of case studies c and d. Therefore, the knowledge of N e and T e /T i or N e and T e provides more information than the use of the confined search range. Furthermore, almost no performance improvement was obtained in the case of knowing a priori or N e and T e parameters. This suggests that the addition of those parameters already solved the ambiguity problem. Even so, to avoid "incorrect" estimates in case these ion composition assumptions were invalid, it is suggested to fit plasma parameters always using the full search range of parameters. The following sections show the probability results that were obtained using the full search range.

Impact of PSO parameters choice
To evaluate the impact of the selection of the PSO internal parameters, Monte Carlo simulations were also done Page 13 of 25 Martínez-Ledesma and Jaramillo Montoya Earth, Planets and Space (2020) 72:172 at different signal fluctuation levels using the configuration of parameters "Param. 2" (defined by . Figures included in Additional file 2 show the histograms of χ 2 r and the error of the ion composition estimate of case studies a, b, c, and d using this configuration of parameters. Figure 6 shows P fit valid and P correct values obtained from the simulations of the previously indicated combinations of plasma parameters known  Probability of convergence ( P fit valid ) and probability of 'correct' estimation ( P correct ) (in percentage) obtained by simulations (black continuous lines) using the MCLM method, and using the PSO algorithm with the parameter configuration (dash-dotted lines) "Param. 1" used by Wu et al. (2015) and (dashed lines) "Param. 2". Simulations were done at different signal fluctuation percentages ( δ(%) ) fitting different combinations of known a priori plasma parameters from the plasma line: without a priori information, given N e , given N e and T e /T i , and given N e and T e parameters Page 15 of 25 Martínez-Ledesma and Jaramillo Montoya Earth, Planets and Space (2020) 72:172 Simulation results without a priori information (case study a) showed different probabilities of convergence using the configuration of parameters "Param. 2". In this case, this configuration obtained an increase of P fit valid between 14 and 22% at signal fluctuations smaller than δ < 0.5% . Nevertheless, the other case studies (i.e., b, c, and d) showed almost identical values of P fit valid for both configurations of PSO parameters. Moreover, P correct values were almost identical in all case studies. The unambiguous estimation probability differences obtained in case studies b, c, and d were smaller than 1%, but differences up to 5% were found in simulations without a priori information. These results indicate that the use of a different configuration of PSO parameters did not affect the unambiguous estimation performance of the algorithm, but slightly increased the convergence of solutions when no a priori information was provided. Figure 7 shows the average standard deviation of both "correct" and "incorrect" distributions of the GMM PDFs calculated by the EM algorithm ( σ correct and σ incorrect , respectively) for simulations using the configuration of PSO parameters "Param. 1" and "Param. 2" (left and right columns, respectively). In small noise regimes (signal fluctuations smaller than δ < 2%), standard deviations increased with the signal fluctuation percentage. In this signal fluctuation region, estimated regressions of standard deviation ( σ est ) were obtained. In higher noise regimes, an almost constant standard deviation ( σ sat ) was obtained with a value of ~ 0.16. The signal fluctuation level where the standard deviation reaches this saturation value is approximately equal to δ cross , which is the signal fluctuation level at which P correct = 50% in Fig. 6. In case studies a and b, the estimated σ est were slightly smaller when the configuration "Param. 2" was used. Furthermore, in case study a (without a priori information) at small signal fluctuations the values of "correct" standard deviations were much larger than the "incorrect" values, an effect not obtained by simulations with information provided from the plasma line. Alternatively, identical estimates of σ est and δ cross were obtained in case studies c and d using both configurations of parameters. These results suggest that the estimation accuracy (i.e., ion composition estimation uncertainty) of the PSO algorithm did not improve significantly using different configurations of parameters. Moreover, when no a priori information was provided, larger estimation uncertainties were obtained by the optimization algorithm at small signal fluctuations.

PSO algorithm computation time
The number of iterations and computing times of the PSO algorithm for simulations using both configurations of PSO parameters ("Param. 1" and "Param. 2") are shown in Fig. 8. The average number of iterations (Fig. 8,top) obtained by simulations using the configuration "Param. 1" were equal to the maximum number of iterations ( max_iterations= 500 plus 1 for the determination of the initial random particle locations) in all case studies. Alternatively, simulations using the configuration "Param. 2" presented a smaller number of iterations. The number of iterations obtained was gradually reduced by the addition of different a priori information in each case study. Simulations assuming known a priori N e and T e parameters obtained the minimum number of iterations (~ 472 iterations). The maximum reduction of iterations was found in the high signal fluctuation range in all case studies. This reduction of iterations is related to the increase of stagnation detections obtained at high signal fluctuations when using the configuration "Param. 2", shown in Fig. 9. Note that, at high signal fluctuations almost all solutions were convergent ( P fit valid ≈ 100% in Fig. 6), indicating that stagnated solutions detected were the minimums of the cost function. Therefore, this result suggests that the configuration "Param. 2" improves the optimum detection speed of the PSO algorithm. Nevertheless, the reduction of iterations obtained by the stagnation detections did not generate a large impact on the overall computing time, as can be seen in Fig. 8 (middle). This effect is related to the large number of minimum iterations of the algorithm required to determine stagnation ( max_count_const+count_const_init= 300). As this number is large relative to max_iterations , the detection of stagnation cases generated only a small reduction in the total number of iterations.
Nevertheless, computing times were smaller when using the configuration of parameters "Param. 2" in all case studies. Boxplots of computing times for each case study are shown in Fig. 8 (bottom). Average computing time reductions were of 8.37%, 5.76%, 6.74%, and 5.96% for case studies a, b, c, and d, respectively. Therefore, this result indicates that the configuration of PSO parameters "Param. 2" obtains a reduction in the overall computing time of the PSO algorithm.

Comparison with MCLM optimization algorithm
Probabilities obtained by the MCLM optimization algorithm (Martínez-Ledesma and Díaz 2019) shown in Fig. 6 (continuous black lines) allow a visual comparison of the estimation performance with the PSO algorithm. Results of PSO simulations without a priori information (case study a) obtained much lower convergences ( P fit valid ) than MCLM simulations at very small signal fluctuations (with probability differences up to ~ 44% or ~ 58% depending on the configuration of PSO parameters). On Page 18 of 25 Martínez-Ledesma and Jaramillo Montoya Earth, Planets and Space (2020) 72:172 the other hand, PSO simulations of case studies b and c obtained larger P fit valid values than MCLM simulations, with differences of ~ 8% and ~ 14%, respectively. Finally, the PSO simulations of case study d were almost perfectly convergent ( P fit valid ≈ 100% ), as obtained by the MCLM algorithm. Of relevance is that the simulations of case study c have an almost perfect convergence using the PSO algorithm, while the use of the MCLM provided P fit valid ≈ 85% at small signal fluctuations. Therefore, the PSO algorithm improves the convergence of the MCLM algorithm when adding a priori information, particularly when assuming known a priori N e and T e /T i parameters. Approximately similar unambiguous estimation performance ( P correct ) was obtained by the PSO algorithm and the MCLM method. The most remarkable result was obtained in the case of simulations knowing a priori N e and T e parameters (case study d). In this case, identical values of P correct were obtained using both optimization algorithms. These probability values were almost perfectly "correct" up to high signal fluctuations. This result suggests that the use of the a priori information of N e and T e parameters solves the ambiguity and provides the maximum unambiguous estimation probability that could be obtained by any optimization algorithm, as indicated by Nicolls et al. (2006). On the other hand, the PSO simulation results without adding a priori information (case study a) obtained probabilities similar to the MCLM. In this case study, probability differences between PSO and MCLM algorithms were between ± 6%, and between − 2.8% and 4% when using the configuration of PSO parameters "Param. 1" and "Param. 2", respectively. Also, PSO simulation results of case studies b and c were identical to those of the MCLM method for signal fluctuations δ ≤ 0.2% . In those study cases, probability differences obtained were below 7%, with maximum differences located at signal fluctuations δ=2% and 5%, respectively. The signal fluctuation thresholds for obtaining a P correct ≥ 95.45% (a 2σ criterion) of the PSO algorithm were identical to those obtained by the MCLM algorithm only in case studies b and d ( δ th(Ne) PSO = 0.14% and δ th(Ne and Te)PSO = 7.9%, respectively). However, PSO simulations without adding a priori information obtained a value of P correct ≈ 95% at signal fluctuations δ = 0.01%, which is smaller than the global threshold ( δ th(no a priori) MCLM = 0.05%) obtained by Martínez-Ledesma and Díaz (2019). It is assumed that this small threshold value is related to the lower convergence of the algorithm at this signal fluctuation level ( P fit valid (δ = 0.01%) ≈ 10% or 25%, depending on the configuration of PSO parameters used). Aside, the signal fluctuation threshold obtained in simulations of the PSO algorithm knowing a priori N e and T e /T i parameters (case study c) was approximately δ th(Ne and Te/Ti) PSO = 1.4%. In this case study, the threshold value when using the MCLM algorithm was δ th(Ne and Te/Ti)MCLM = 0.568%. Therefore, the PSO algorithm improved the unambiguous estimation threshold of the MCLM algorithm when adding a priori N e and T e /T i parameters.
On the other hand, in some case studies, the joint probability P fit valid & correct obtained a much larger improvement. Of relevance is the increase in probability of ~ 7.3% and ~ 14% in case studies b and c, respectively. Simulations of the PSO algorithm knowing a priori N e and T e /T i parameters showed almost perfect P fit valid & correct values up to this signal fluctuation threshold value, as in the case of knowing a priori N e and T e parameters. Therefore, in these two cases of adding a priori information, the convergence of the PSO algorithm is not required to be verified.
The estimated tendencies of standard deviation of "correct" distributions calculated from simulations of the MCLM optimization algorithm ( σ MCLM obtained from Fig. 6 of Martínez-Ledesma and Díaz 2019) are shown in Fig. 7. In case studies b, c, and d, the results of simulations with the PSO algorithm have similar or identical increase tendencies to the results of the MCLM algorithm. Alternatively, when no a priori information was provided (case study a) the average values of "correct" standard deviation of the PSO algorithm were much larger than the results of the MCLM algorithm at small signal fluctuations. Furthermore, "incorrect" standard deviations obtained by the PSO algorithm were similar to the estimated standard deviation values of the MCLM algorithm. This result suggests that the PSO algorithm without a priori information obtained "correct" estimations with larger uncertainties than the MCLM algorithm at the small signal fluctuation range.
Moreover, PSO simulations were done to analyze the impact of using a priori information with uncertainty. The resulting probabilities of convergence and unambiguous estimation are shown in Fig. 10. These Monte Carlo simulations were done with the configuration of PSO parameters "Param. 2" at different signal fluctuation levels and with different uncertainty levels on the a priori parameters. To allow a direct comparison, this figure also shows the probabilities obtained by the MCLM method (gray dotted lines) when considering the uncertainty of the a priori parameters (obtained from Fig. 9 of Martínez-Ledesma and Díaz, 2019). Very similar probability results were obtained by both the PSO algorithm and the MCLM method. As indicated in Martínez-Ledesma and Díaz (2019), these results verify that the uncertainty of a priori parameters should be smaller than the signal fluctuation level of the ISR signal ( δ ≥ |ǫ| ) to ensure the convergence of the algorithm. Moreover, the use of a priori parameters with uncertainties |ǫ| ≤ 1% provide P correct values almost identical to simulations without uncertainty. Finally, the results at low signal fluctuation verify the existence of a global signal fluctuation threshold δ th = 0.05% to completely solve the estimation ambiguity ( P correct = 100% ) independently of the type of a priori information provided and the level of uncertainty of the a priori parameters.
As in the case of simulations without uncertainty, identical probabilities of convergence and of unambiguous estimation were obtained by both optimization algorithms from simulations knowing a priori N e and T e parameters (study case d). Nevertheless, some differences were found at 0.1% ≤ δ ≤ 10% in simulations with large uncertainties ( ǫ = ±50% and ǫ = ±100% ) due to the small number of convergent solutions obtained. Also, the convergence rates ( P fit valid ) obtained by the Page 20 of 25 Martínez-Ledesma and Jaramillo Montoya Earth, Planets and Space (2020) 72:172 PSO algorithm were larger than those of the MCLM method when assuming a priori known N e and N e and T e /T i (study cases b and c). Moreover, in these study cases, slightly larger probabilities of unambiguous estimation ( P correct ) were obtained at δ > 0.2% by the PSO algorithm. The average iterations required by the L-M algorithm used by the MCLM method are shown in Supporting Information Figure S7 of Martínez-Ledesma and Díaz (2019). Computing times of each fitting of the L-M algorithm were always less than 0.45 s. Also, computing times were reduced when adding information from the plasma line and when estimating radar signals with lower noise characteristics. Simulations of case studies b, c, and d were approximately 1.2, 2.3, and 4 times faster than simulations without a priori information, respectively. On the other hand, the computing times obtained by our vectorized implementation of PSO algorithm were between 6.5 and 7.7 s for each fit, depending on the addition of a priori information and the internal configuration of the PSO algorithm parameters (Fig. 8, bottom). The lack of computing time information in Wu et al. (2015) prevents the direct comparison with our PSO implementation. Furthermore, no stagnation detection was implemented in the study of Wu et al. (2015). Figure 11 shows the ratio between computing times of fittings using the PSO algorithm relative to fittings of the L-M algorithm. These ratios were highly dependent on the a priori information provided and the signal fluctuation level. The L-M algorithm was between 15 and 250 times faster than the PSO algorithm at the high and low signal fluctuation ranges, respectively. For noise Fig. 10 (Continuous colored lines) probability of convergence ( P fit valid ) and probability of unambiguous estimation ( P correct ) (in percentage) obtained by simulations considering uncertainty on the a priori parameters known from the plasma line: (left) given N e , (middle) given N e and T e /T i , and (right) given N e and T e parameters and at different signal fluctuation percentages ( δ(%) ). Simulations were done using the PSO algorithm with the parameter configuration "Param. 2". Probabilities obtained without uncertainty (shown in Fig. 6) are shown with dashed black lines. Results from simulations using the MCLM method (dotted gray lines) are shown for comparison Page 21 of 25 Martínez-Ledesma and Jaramillo Montoya Earth, Planets and Space (2020) 72:172 levels smaller than the signal fluctuation thresholds of the MCLM method (shown as vertical bars in Fig. 11), the probability of unambiguous estimation is high enough ( P correct ≥ 95.45% ) to assume that all estimates are "correct". In these cases, the MCLM method only requires evaluating one L-M fitting. On the other hand, for noise levels larger than the signal fluctuation thresholds, multiple fittings of the L-M algorithm are required to be done with different initial parameters of the search to select the "correct" solution in the MCLM method. In the Monte Carlo simulations of Martínez-Ledesma and Díaz (2019), fittings were repeated 500 times. This large number of repetitions was required to determine with high accuracy the unambiguous estimation performance of the algorithm, obtaining high quality statistical distributions of "correct" and "incorrect" estimates that were clustered by the EM algorithm. Also, the executions of the multiple L-M fittings of Martínez-Ledesma and Díaz (2019) were parallelized in 8 different CPU nodes to reduce the global computing time of the MCLM simulations, as in our implementation of the Monte Carlo simulations of the PSO algorithm. Therefore, the theoretical computing time of the MCLM method can be roughly approximated as T MCLM = T L−M · R/P (where T L−M is the L-M algorithm computing time, R is the number of repetitions of the fitting, and P is the number of parallelized nodes used). It is relevant to indicate that Martínez-Ledesma and Díaz (2019) did not suggest the number of fitting repetitions required by the MCLM algorithm to operationally estimate real radar measurements. Hence, we consider using 100 fitting repetitions in the MCLM method, although a much smaller number would be required to provide good estimates of the statistical distributions of "correct" and "incorrect" solutions. This number of repetitions of L-M fittings in the MCLM method is identical to the number of particles used in our implementation of the PSO algorithm, allowing a direct comparison between the computation requirements of both optimization algorithms. Therefore, the sequential execution of the MCLM method would require 100 times the L-M algorithm computing times at noise levels higher than the signal fluctuation thresholds. This result indicates that the computing times of the sequential execution of the MCLM are larger than the computing times of the vectorized PSO algorithm (Fig. 11). Nevertheless, if we parallelize the execution of the multiple fittings in 8 CPU nodes, the MCLM algorithm would only require 12.5 times the L-M algorithm computing times, which is much lower than the minimum computing time ratio of the PSO algorithm in any case study. Although further work is required to verify the computing time improvement of the parallelization of the PSO algorithm, the execution of the parallelized MCLM method would require lower computing times than the parallelized PSO algorithm. This is because the L-M algorithm requires a much lower number of iterations than the PSO algorithm to obtain a solution (as shown in Figure S7 of Martínez-Ledesma and Díaz, 2019).

Fig. 11
Ratio between computing times of fittings using the PSO algorithm relative to fittings of the L-M algorithm at different signal fluctuation percentages ( δ(%) ). Vertical bars indicate the signal fluctuation thresholds of the MCLM method, for which the MCLM method should require only one L-M fit with random initial parameters