{'fraction_non_alphanumeric': 0.03746978819558657, 'fraction_numerical': 0.010339945368476752, 'mean_word_length': 4.571189546649669, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Distributed Information-based Source Seeking

In this paper, we design an information-based multi-robot source seeking algorithm where a group of mobile sensors localizes and moves close to a single source using only local range-based measurements. In the algorithm, the mobile sensors perform source identification/localization to estimate the source location; meanwhile, they move to new locations to maximize the Fisher information about the source contained in the sensor measurements. In doing so, they improve the source location estimate and move closer to the source. Our algorithm is superior in convergence speed compared with traditional field climbing algorithms, is flexible in the measurement model and the choice of information metric, and is robust to measurement model errors. Moreover, we provide a fully distributed version of our algorithm, where each sensor decides its own actions and only shares information with its neighbors through a sparse communication network. We perform intensive simulation experiments to test our algorithms on large-scale systems and physical experiments on small ground vehicles with light sensors, demonstrating success in seeking a light source.


I. INTRODUCTION
Multi-agent source seeking is a robotics task that uses autonomous vehicles with sensors to locate a source of interest whose position is unknown. The source of interest can be a light source [1], a radio signal transmitter [2], or a chemical leakage point [3]. The source seeking vehicles, or mobile sensors, can measure the source's influence on the environment and use this information to locate the source.
A large body of source seeking research investigates field climbing methods [4]- [9]. Assuming the source signal gets stronger as the sensor-source distance shortens, the mobile sensors can "climb" the source signal field to physically approach the source. These methods do not require explicit knowledge of the measurement model, making them easy to implement for different applications. However, field climbing methods are not necessarily the most effective source seeking methods. Firstly, they only exploit local information of the source field, with the typical requirement that sensors must maintain a tight formation to make a reasonable ascent direction estimate, as is the case in [6], [7]. Furthermore, the sensors cannot move too fast as a group for the measurement value to achieve a stable increase.
An alternative approach is to perform source identification/localization using various estimation methods, such as (Extended) Kalman Filter (EKF) [10], [11], Particle Filter (PF) [12], [13], and so on, to keep a constantly updated estimate of the source location over time. These methods enable the fusion of measurements from multiple sensors to The   identify a global view of the measurement field. A central theme of these works is to improve the estimation through sensor movements, and studies on optimal sensor placement propose to optimize specific information metrics, in particular, variants of Fisher information measures [2], [3], [11], [12], [14]- [17], which relate closely to the famous Cramér-Rao lower bound [18], [19]. However, these studies typically focus on deriving closed-form solutions of optimal sensor placement for particular types of measurement models [2], [12], [16]. Therefore, results often do not generalize, not to mention that closed-form solutions might not exist for many general measurement models. Besides, the methods are often not robust to the modeling error.
This paper draws advantages from the two methods mentioned above (field climbing and source localization) to develop multi-robot source-seeking algorithms. In particular, we assume the measurement model is available and propose an algorithm using range-based measurements. Each iteration of the algorithm consists of three steps: 1) Collect range-based measurements.
3) Move along the negative gradient direction of information-based loss function L, defined as the trace of the inverse of Fisher information. In particular, step 3) seeks to increase the information about the source location in the measurement, rather than looking for stronger source signals as in field climbing. Therefore, we name our method as information-based source seeking.
Our contributions: We first introduce our information-based method in the centralized setting (Section III), where we formally define the loss function L and the three-step outline of our method. We also provide theoretical justifications about the choice of loss function L, as minimizing L improves the estimation and gets the sensors close to the source location. Then we extend our method to the distributed setting (Section IV), which makes our method scalable in the number of sensors. We show how distributed estimation and gradient calculation can be effectively done using simple consensus schemes. In Section V we conduct extensive numerical experiments to study the performance of our method in both centralized and distributed settings. The results show that our method outperforms field climbing methods and is flexible and robust in multiple aspects. Lastly, we implement our algorithm on small ground vehicles carrying light sensors to seek a light source in a dark room (Section V-E), as shown in Fig. 1. The hardware implementation further demonstrates the effectiveness of our algorithm.
The advantages of our methods can be summarized as • Compared with field climbing algorithms, numerical studies quantitatively demonstrate that our algorithm converges much faster to the source and performs more consistently over repeated trials. See V-A2. The algorithm takes advantage of the sensing capacity of multi-sensors in the sense that the performance improves as the number of sensors increases. • Our algorithm, especially its gradient-guided movement, provides flexibility in handling various measurement models and picking different information metrics as loss functions. Section V-A3 confirms at least three applicable metrics. • The algorithm is more robust to modeling error than the source localization with stationary sensors, see V-A4. • Our distributed algorithm achieves comparable performance as the centralized algorithm. See V-B. Moreover, results show that the distributed algorithm is more robust than the centralized implementation to the error in initial guesses and communication delay. See V-C, V-D.
We note that this paper is the journal extension to our previous conference work [20]. Compared to [20], besides including more technical details for the performance guarantee of the methods, this paper extends the centralized algorithm in [20] to the distributed setting, where a central controller is absent, and the sensors decide their actions individually, as described in Section IV. Numerical results in Section V-B∼V-D provide detailed comparisons between the centralized and distributed algorithm. Lastly, the hardware implementation migrates from the centralized robotics platform in the previous work to a distributed robotics platform in this paper, with the programming environment switching from ROS 1 to ROS 2.

A. Related Work
Field climbing methods. Field climbing originates from scientific studies of animal behavior in exploring nutrition or chemical concentration fields [21], [22]. The studies inspire source seeking methods that use field value measurements to estimate the field gradient and apply formation control to climb along the gradient [5], [6]. Gradient-free field climbing methods are also studied in the literature: [7] maintains the sensors in a circular formation and uses measured field values as directional weights to guide the overall movement. Later works such as [8], [9] extend the algorithms above to the distributed setting by applying consensus in the calculation of ascent direction. Single-agent field climbing is also possible if a gradient ascent control law is combined with proper zerothorder gradient estimation algorithms, for example, using the control framework in [23].
The main differences between field climbing and our algorithm are that i) field climbing maximizes the source field value directly, while our algorithm exploits the Fisher information, an indicator of both estimation accuracy and source-sensor distance, ii) field climbing does not assume a given measurement function, but ours does, iii) field climbing algorithms often require a tight sensor formation for a stable field ascent, which only exploits local information; while our algorithm uses the sensors to collect global information for source localization. We will show in Section V that the sensors under our algorithm tend to spread out so that measurements contain more diverse information about the source location.
Source localization and optimal sensor placement. If a measurement model is available, the source location can be estimated using methods including EKF [10], [11], [24], [25], PF [12], [13], and so on. It is even possible to reconstruct the entire source field [3], [26]. Such estimation of source location is also known as source localization. An important problem in source localization is how to improve estimation via sensor movements. This problem is often investigated under the optimal sensor placement framework, in which sensors move to optimize various information metrics, including covariance [27], [28], mutual information [13], [29], and Fisher information measures [2], [3], [11], [12], [14]- [17], [24], [25] such as the determinant of Fisher information, the largest eigenvalue of its inverse, and the trace of its inverse (D-, E-, and Aoptimality criterion respectively).
Our method also employs Fisher information but is different from the above as follows: i) Most of the previous works focus on deriving closed-form solutions of optimal sensor placement for a particular type of measurement model, for example, RSS [2], pollutant diffusion [12], and gamma camera [16]. Our focus is not on providing closed-form solutions. Instead, we provide a gradient-based method applicable to a large class of range-based measurement models. This method also provides flexibility in choosing different information metrics as loss functions. ii) Many previous works are about finding the optimal angular placement of the sensors at a fixed distance to the source or in a restricted area [2], [11], [14], [17]. In contrast, we allow the sensors to move freely and eventually reach the source. iii) Some studies relax the restrictions on the sensor movement [15], [16]. However, these methods produce a spiraling sensor movement which is inefficient for source seeking purposes. Our method does not produce such movement. iv) The line of work by [24], [25] also allows the sensors to move freely, but the control framework is very different from our work. The proposed ergodic control, which involves improving a control objective in integral form, is conceptually much more complex than our method-a gradientbased control. Ultimately, the ergodic control is designed for a very general set of mobile sensor applications, while our work is catered for the source seeking problem.
Bayesian inference and optimization. The recent advances in Bayesian learning inspire many source seeking studies to adopt the Bayesian methods [30]- [32]. These studies view the environment as a field characterized by an (unknown) density function related to measurement and use Gaussian Process or other likelihood models as a surrogate to guide the sensor movements for new measurement collections. The computation (for running the posterior update and Bayesian optimization) and memory (for storing historical measurements as in non-parametric Bayesian methods) demands of these methods are usually much higher than our method. Overall, Bayesian-based methods and our method are designed from different principles. A detailed comparison is left for future work.
Lastly, it is worth mentioning that our work is significantly inspired by [11], which studies the optimal sensor placement problem on a surveillance boundary. We leverage their ideas to introduce the Fisher information in our objective, but we change the loss function from the determinant of Fisher information to the trace of its inverse. We also generalize the measurement model to fit into our experiments and other real-world measurement settings. To handle such more general measurement models, the sensors only compute the gradient of the loss function rather than solving for its optimum at each time step.

II. THE PROBLEM STATEMENT
Consider the problem of using a team of mobile sensors to find a source whose position is unknown. The source can be any object of interest, such as a lamp in a dark room. Mobile sensors are robotic vehicles that can measure the influence of the source on the environment, such as small ground vehicles carrying light sensors.
Specifically, we use q ∈ R k to denote the source position, and use p 1 , p 2 , · · · , p m ∈ R k to denote the sensor positions, where m is the number of sensors and k is the spatial dimension. The measurement y i made by the ith sensor is modeled by where h i is a known continuously differentiable function and ν i is the measurement noise. The value of h i depends on the position of sensor i but not on other sensors, and different sensors may have different h i . We let y = [y 1 , y 2 , . . . , y m ] denote the vector of all measurement values, let ν = [ν 1 , ν 2 , . . . , ν m ] be the noise vector, and use p = p 1 p 2 · · · p m to denote the joint location of all mobile sensors. We let H : R mk × R k → R m be the mapping that describes the joint measurement made by all the sensors, i.e., . . .
We define the global measurement model as This model relates all measurements to the source location and all sensor locations.
The source seeking objective is to have at least one of the sensors get within 0 distance to the source. In other words, we want to achieve where 0 is some small positive number.

III. INFORMATION-BASED SOURCE SEEKING
This section presents the key components of our information-based source seeking, including the three-step algorithmic flow and the loss function design. We assume a centralized setting for the conciseness of description, meaning that there is a central server that gathers all sensor locations p and measurements y, and decides new locations for all the sensors. This assumption will be removed in Section IV where we introduce the distributed setting.

A. The Algorithm
Our information-based source seeking consists of three consecutive steps in one iteration: measurement, source location estimation, and sensor movement. The algorithm is illustrated in Figure 2 and detailed in Algorithm 1.  Estimate the location of the source byq ← E(y, p).
In our implementations, the Extended Kalman Filter (EKF) is employed as the estimation algorithm. We define z := [q , v ] as the source state, where v is the velocity of the source. We do not assume any prior knowledge about the motion of the source, except that it evolves according to the second-order dynamics defined by: Equation (5) models the source moving at some unknown but constant velocity v. The applicability of this motion model is backed by extensive practical evidence in GPS problems [33], [34] and moving target tracking [35]. In general, it is sensible to make a minimal assumption as above on the motion model in source localization. Another example is to define the state z := q so that the transition model is z + = f (z) := z. See the third paragraph of section IV.B in [24]. However, we found the estimation result from this model is not as stable as using the model in (5). At each time step, the EKF takes in the sensor measurements and sensor locations and returns an estimateẑ of the true source state z. Please see [36,Definition 3.1] for details about the EKF update.
We will explain the motivation of this loss function later in Section III-B, by relating it to the Fisher information matrix and the Cramér-Rao lower bound [18], [19]. Sensor movement includes two parts: waypoint planning and waypoint tracking. a) Waypoint planning: In Lines 5-8, with the source location estimateq obtained in Line 4, we generate a set of waypointsp(0), . . . ,p(T ) by applying gradient descent on L with respect to p, with q =q fixed. In Line 7, α t > 0 is the step size and M t 0 is a directional regularization matrix. b) Waypoint tracking: In Lines 9-10, after the joint waypointsp(0), . . . ,p(T ) have been calculated by the gradient descent on L, we extract the waypointsp i (0), . . . ,p i (T ) for each mobile sensor i, and use the motion planner M P to calculate a sequence of control inputs u i (1), . . . , u i (T ) and apply the control at the first instance u i (1) to the corresponding mobile sensor. 1 Remark 1. Note the algorithm aims to minimize L(p,q) rather than L(p, q), which raises the question of whether improving the former leads to a decrease in the value of the latter. Intuitively, it is expected that ifq and q are sufficiently close to each other, i.e., the estimation error is small, then improving L(p,q) will decrease L(p, q). These two algorithm components rely on each other to function as a whole.

B. Information-based Loss Function for Sensor Movement
One naíve approach to determine sensor movement after Step 2) could be heading directly towards the estimated location. However, we observed that such movement often makes sensors cluster tightly, which reduces the diversity in measurements and is undesirable. As we later show in Fig. 3c and 4c, the resulting estimation and seeking are unstable.
Creating a good condition for estimation through sensor movements is crucial for successful source seeking. Therefore, instead of following the haphazard approach above, our algorithm takes the gradient steps in Step 3). The idea is to let the sensors maximize some information metrics related to the source location. In specific, we consider the Fisher information about the source location, which is defined as follows: let y, p be the joint-measurement and joint-position of all the sensors, and H be the corresponding global measurement function, as in (3). We treat y as a random vector, with joint probability density f (y; q) parameterized by the true source location q. The Fisher information(FIM) about q is considered to measure the amount of information about the unknown q contained in measurement y. It is defined as Fisher information is closely connected with estimation quality, as it is a lower bound on the covariance of unbiased estimators. This result is formally known as the Cramér-Rao Lower Bound.
Theorem 1 (Cramér-Rao Lower Bound (CRLB) [18], [19]). For any unbiased estimatorq of q, the following matrix inequality holds As a special case, if we assume y follows multi-variate Normal distribution described by then the Fisher information takes the form of Under the further assumption that Motivated by the form of FIM in (11), we use as the loss function. Note that by minimizing (12), the CRLB is driven closer to the zero matrix. Furthermore, if H satisfies Assumption 1, then minimizing L also results in approaching the source as stated in Propositions 1 and 2.
Assumption 1. We make the following assumptions on the measurement functions h i : 1) Isotropic measurement: The measurement values depend only on source-sensor distance, i.e., for some function g i : (0, +∞) → R, where r i := p i − q . 2) Monotonicity: The absolute value of the derivative of each function, |g i (r)|, is monotonically decreasing in r. 3) Non-degeneracy: Letr i denote the unit direction vector from the source to the i-th mobile sensor, i.e.,r i = ( Proposition 1 (Reaching the Source I). Under Assumption 1, which means that if L(p, q) decreases, then max i |g i (r i )| tends to increase and consequently min i r i tends to decrease. Please see Appendix B for the proof.
In Proposition 2, we assume the p(t) sequence converges for technical convenience, but it is also reasonable to assume the sensors will reach a stationary formation when they have minimized L. Also, the infimum of L is taken over the set{p : p i = q, ∀i} because the FIM could be undefined if p i = q for some i.
Remark 2. Assumption 1.3 is a regularity condition that ensures L is well-defined, which is essential for the analysis of the relationship between L and the source-sensor distance. We enforce this assumption in the implementations by adding a positive definition matrix δI, with δ > 0 being a small constant, to the estimated FIM when calculating L and its gradient. This technique proves to eliminate most of the pathological numerical behaviors while maintaining the soundness of source-seeking performance. We have added this clarification to Remark 2 of the revised paper.
Remark 3. The full analysis of Algorithm 1 requires an effort that is beyond the scope of the current paper and must be left for a future paper, but we will outline the high-level approach to deriving a convergence guarantee below. In a nutshell, the analysis accounts for optimizing a parameterized objective function, where the parameter of the objective function needs to be estimated online. In the analysis, we consider the Gradient Descent and EKF components as a whole and analyze the convergence of this combined system. Suppose {δ t := ||q t − q||} t=0,1,2,... are the estimation errors at each time step t. We show that (1) For any > 0, there isδ > 0, such that if the δ t ≤δ for all t, the gradient descent on L with appropriate step sizes can converge to the point that is −close to the source location.
For the above to hold, there must be a bounded neighborhood of the source q where there are no other stationary points, and the gradient of L is Lipschitz. The initial state of the agents must be within this neighborhood, and the step size on gradient descent should be chosen according to the Lipschitz constant as well as the error boundδ. The specifics of this item are mainly based on standard results in non-linear optimization [38], [39], and also the online convex optimization(OCO) literature studying the properties of gradient descent under imperfect knowledge about the loss function, see [40], [41] and references therein. (2) In item (1), the assumption that δ t ≤δ for all t may seem too limiting, but if the estimation error of EKF decreases fast enough so that δ t ≤δ for most of the time steps, the conclusion in (1) still holds. To prove the above, the linearized model for EKF must have bounded system norms, the measurement function must be Lipschitz within a bounded neighborhood of the source location, and the initial estimation error must be small enough. The specifics of this item are given in [36]. (3) Note that the convergence of error in item (2) does not consider the fact that item (1) contributes to a lower CRLB, which is beneficial to estimation. So far, we do not have a definite answer on how to incorporate this factor to potentially sharpen the convergence rate in (2). However, it will be an interesting direction to look into.

IV. DISTRIBUTED INFORMATION-BASED SOURCE SEEKING
In this section, we present our distributed source seeking algorithm. We no longer assume a central controller in the distributed setting. Instead, the mobile sensors decide on their actions individually. Meanwhile, the sensors can communicate through a network to exchange information with their neighbors about measurements, sensor positions, and estimations. We model the communication network by a directed graph, denoted by G = (N , E), where N is the set of nodes(sensors) and E is the set of edges. If an edge (i, j) ∈ E, then sensor j can directly receive information from sensor i. Let Here m j denotes the number of its incoming neighbors except itself. For each (i, j) ∈ E, there is a consensus weight w ji > 0 satisfying i∈Ni w ji = 1 for all j (see [42, Assumption 2.1-2.3] for conditions on the consensus weights).
We assume the information available to sensor j includes those from its measurements and direct communication. In particular, mobile sensor j knows the following information.
1) The measurement functions in its neighborhood, or equivalently the local joint measurement function H j 2) The joint position and joint measurement of its neighborhood(in bold font).
We assume the local joint measurement model for j is where ν j = [ν j , ν j2 , ..., ν jw ] is a random vector following multi-variate Normal distribution N (0, R j ).
We assume the sensors act cooperatively, that is their objective is still to have some members in the team to reach the source, as defined in (4).

A. The Distributed Algorithm
Algorithm 2 Distributed Information-based Source Seeking(for sensor j) Take measurements and communicate with neighbors to obtain p j , y j , {ẑ k ,F k } k∈Nj .

11:
Execute control action u j (1). 12: until ||p j −q|| ≤ 0 In the distributed algorithm, each sensor j maintains an estimate for the source location and velocity, denoted bŷ z j := [q j ,v j ] . The sensor also maintains an estimate for the global Fisher information, denoted byF j . Our source seeking algorithm requires sensor j to repeat three consecutive steps: information gathering, source location estimation, and movement. Algorithm 2 contains the detailed operations.
1) Information Gathering (Line 3): Sensor j takes measurements and communicates with its neighbors to obtain the sensor positions p j and measurements y j . Sensors also communicate their source location estimatesẑ j and their Fisher Information estimatesF j , as required in the subsequent steps.
2) Source Location Estimation (Line 4): The gathered information is used to form an updated source location estimateq j , through consensus-based estimation algorithms, elaborated in Section IV-B. For ease of exposition, we assume the source location estimator to be the consensus EKF (20), though other estimators could also be adapted to the distributed setting following a similar procedure.
3) Sensor Movement (Line 6-11): Sensor j calculates the partial gradient of the loss function L (6) through Distributed Information Gradient Update (Line 6-9), elaborated in Section IV-C . Then it calls a motion planner MP to compute the corresponding control actions for the gradient steps and executes the first action (Line 10-11).
The key differences between the distributed Algorithm 2 and Algorithm 1 are: • Algorithm 1 runs on the central controller. The distributed algorithm runs on individual sensors. • Comparing Step 1), the Algorithm 1 knows the measurements and positions of all sensors. The distributed algorithm only knows those from the neighborhood. As a result, the Algorithm 2 modifies Steps 2) and 3) to accommodate the distributed setting. Next, we elaborate on how agents perform distributed source location estimation and distributed information gradient update in Algorithm 2.

B. Distributed Source Location Estimation
In Step 2), each sensor must locally form a source location estimate based on the information available to it, and one simple way is to run local EKF based on its local information. For each sensor j, the local EKF maintains two running variables,ẑ j and P j , wherê is sensor j's estimate of the position-velocity vector of the source, and P j is a matrix quantifying the uncertainty in the estimate. These two variables are updated according tô where Q j , R j are pre-defined constant positive matrices, f is given in (5), and H j is given in (15). Q j and R j represent our prior belief of the covariances of process and measurement noises, therefore are usually the same among all j. The quantities A j , C j , and K j are defined by The local EKF may be one of the simplest estimation algorithms based on sensors' local information from neighboring sensors. However, its estimation is not great in our source seeking experiments, as we later show in V-B. One of the main reasons is that the communication network is typically sparse, and each sensor is connected to only a few neighbors. Therefore the measurements available to each sensor are few, and the resulting estimation is poor.
The consensus EKF, an extension to the Kalman Consensus Filter II algorithm [43], addresses the limitations of local EKF by introducing a consensus procedure among the estimates. The estimationẑ j is updated with the following formulâ with w ji ≥ 0, i∈Nj w ji = 1. The update of P j remains the same as the local EKF.
The update in (20) can be viewed as adding a consensus term w jiẑi to the difference between local EKF update and currentẑ j . Meanwhile, the local EKF can be viewed as a special case of (20) with w jj = 1 and w ji = 0 for i = j. We highlight here the advantage of the consensus EKF over the local EKF: • Since the neighbors' estimate uses the measurements from possibly outside N j , sensor j is indirectly exposed to more data than its immediate neighborhood. Therefore the overall estimation of the consensus EKF is much better than the local EKF. • With properly chosen consensus weights as specified in [42, Assumption 2.1-2.3], the sensors can reach consensus on estimates much quicker than the rate of change in local EKF update. Effectively, the limiting estimate at consensus is an estimator that uses the measurement from all the sensors in the network. It serves as a target location for the sensor to coordinate their movement, and we observe much better convergence to the source using the consensus EKF. Our consensus EKF algorithm is an extension to the Kalman-Consensus Filter II algorithm [43,Section IV,eq. (20)] to the non-linear measurement setting. A number of extended Kalman-like filters that also employ consensus share some similarities with our algorithm [44]- [46], but they take consensus over both the mean estimate and covariance estimate, while our algorithm only takes an average over the mean estimate to lower the communication load, which improves the responsiveness of source seeking in our real-time environment. This scheme is also backed by the theoretical result in [43,Theorem 2], that performing consensus on the mean estimate but not on the covariance estimate can still lead to convergence in estimation. It is also worth pointing out that although our consensus EKF shares similarities with the existing literature, our consensus-based distributed gradient update, described in Section IV-C, is a novel contribution.

C. Distributed Information Gradient Update
Similar to the centralized setting, our distributed algorithm requires the sensors to move to new locations along the gradient of the Fisher-Information-based loss function L (6). However, since L depends on the global measurement function H and positions p, the challenge is for individual sensors to compute the partial gradients in a distributed fashion. In specific, the partial gradient of L can be shown to be where A j = ∇ q h j (p j , q). Note that both A j and ∇ pj A j can be computed using the local information of sensor j already. Therefore, the remaining problem is to estimate the value of FIM that depends on global information H, p. We note that, meaning that the Fisher information is proportional to the equally weighted sum of the rank-one matrices which we name as partial FIM. Note each F j is again computable using the local information of sensor j. Therefore, the sensors can recover the global Fisher information through consensus. In specific, let each sensor maintain a running estimateF j of the global FIM, and updateF j according tô Here {w ji } i∈Nj are consensus weights, and F − j denotes the partial FIM value from the previous time step. To ensureF j converges to an equally weighted sum of F j 's, the consensus weights can be pre-defined to form a static doubly-stochastic consensus matrix. Alternatively, the update in (23) can be implemented to follow the parallel two-pass algorithm [42].
In summary, sensor j calculates an estimated gradient descent direction by calculating A j and ∇ pj A j locally, while substituting the FIM in (21) withF j , as shown in Line 9 of Algorithm 2. Remark 4. The analysis of Algorithm 2 has to be left for a future paper, but overall it is an extension to the analysis of Algorithm 1, where the information required in gradient computation and estimation is not fully available for each agent but needs to be estimated through consensus. Let ξ t be the maximum error in gradient estimation through consensus among agents at time step t. For each agent i, let {δ i,t := ||q i,t − q||} t=0,1,2,... be its estimation errors at time steps t ≥ 0. We show that (a) For any > 0, there areδ > 0,ξ > 0, such that if the δ i,t ≤δ, ξ t ≤ξ for all (i, t), the gradient descent on L with properly chosen step sizes can converge to the point that is −close to the source location. The additional error resulting from ξ t can be modeled as an adversarial disturbance to the gradient update, and the analysis can apply the techniques in adaptive control with adversarial disturbances. The work by [41] provides an inspirational framework to analyze this item. Regarding the choice of step size, we refer to the techniques in distributed optimization as reviewed in [47]. (b) Similarly, although the assumption that ξ t ≤ξ for all t may seem too limiting, if ξ t decreases fast enough through consensus, the conclusion in (a) still holds. This item requires the consensus weights on gradient consensus to be properly chosen according to [42]. (c) Also, for the estimation errors to satisfy δ i,t ≤δ in the distributed setting, the conditions on EKF 's convergence must be met for each agent in addition to the effect of consensus in the loop. The specifics of this item require the tools in [43] and [36].

A. The advantage of Information-based Source Seeking
This subsection compares Algorithm 1 with field climbing methods via numerical experiments, showing the advantage of information-based source seeking. All experiments are implemented in a centralized way. The algorithm performance under actual robot dynamics is studied in simulations in Section V-A1. In the subsequent studies, we remove the robot dynamics in simulations to efficiently conduct repetitive trials and assume the sensors follow the gradient steps exactly. We study the influence of the number of sensors in Section V-A2, the difference of various information metrics in Section V-A3, and the robustness to modeling error in Section V-A4.
1) Gazebo Numerical Experiments: The following numerical experiments are carried out using the Gazebo simulation toolbox [48], with virtual mobile sensors simulating the same dynamics as the actual robots. We generate simulated measurement values of the sensors by with ν i drawn independently from N (0, 0.01). The measurement function h i (r) = 1/r 2 i is given to the EKF for estimation. Stationary Source In the first set of simulations, we use three mobile sensors to seek a stationary source. The source is fixed at position (6.0, 6.0), while the mobile sensors are initially placed at (1.0, 2.0), (2.0, 2.0), (3.0, 2.0). The initial guess of source location given to the EKF is (3.0, 4.0). The terminal condition threshold 0 = 0.5. We compare the convergence to a stationary source among three algorithms: (a) our algorithm; (b) the field climbing algorithm introduced by [7] that only maximizes measured signal strength; (c) following straight lines to the estimated location. (c) is included to show the importance of exploiting Fisher information in obtaining accurate estimation. The results are displayed in Figure 3.
First, notice that the straight-line algorithm fails to converge to the source, as shown in Fig. 3c. We suspect the reason is sensors cluster together quickly as they move to the (same) estimated location and cannot provide sufficiently rich, diverse measurements for a reasonable estimation. Consequently, the estimate gradually deviates from the source location, as do the sensors. On the other hand, if taking a trajectory that improves the Fisher information, the sensors cover the space more thoroughly, resulting in a stable decrease in the estimation error and the final success of reaching the source, as shown in Figure 3a.
Comparing Fig. 3a and 3b, note that sensors using our algorithm first spread out to estimate the source location better and then converge to the source, whereas sensors doing field climbing maintain a tight formation while steadily approaching the source. Since we use constant rather than diminishing step sizes for 3b, the virtual robots do not stop completely near the source and perform a looping behavior. Although our algorithm and the field climbing algorithm [7] are both successful with a stationary source, our algorithm converges faster consistently over repetitive trials, as is shown later in Section V-A2.
Moving Source In this set of experiments, all parameters are kept the same as the stationary case, except the source moves in a circular motion with constant speed. See Figures 4a, 4b,  4c. Note the straight-line algorithm again leads to a sensor formation that causes the estimation to deviate from the actual source location. Both our algorithm and the field climbing algorithm [7] successfully get close to the source. However, the field climbing method exhibits unnecessary irregular motion when the sensors are near the source. We suspect that the field climbing direction becomes very sensitive to the source movement as the sensors get close to the source, which leads to such motion. In comparison, sensors following our algorithm trace much more stable paths. We later discuss more advantages of our algorithm.
2) Performance with Sensor Swarms: We now investigate the influence of the number of mobile sensors on the algorithm performance. We perform repetitive simulations to study the convergence rate of source-sensor distances and the evolution of estimation error. We use a stationary source fixed at position (6.0, 6.0) and randomize the initial sensor locations following uniform distribution in a 3.0-by-3.0 rectangle, as shown in Fig.  6a.
The first set of experiments studies the convergence rate of the sensors to the source. We look at mobile sensor team performance with sizes of 3, 10, 20, and 50. Our algorithm is compared with two field-climbing methods [6] [7]. The color bands in Fig. 5-6 represent the standard deviation across 100 repetitive experiments. The results show that our algorithm converges much faster than the others.
In addition, we see that increasing the number of sensors substantially reduces all algorithms' variance while not affecting the convergence time. Because the speed of sensors has an upper bound regardless of the number of sensors, the overall convergence rate of source-sensor distance is limited. Nevertheless, having more sensors provides more measurements, which reduces the variance and results in more consistent, stable trials. In particular, our algorithm benefits the most from having more sensors since more measurements contribute to better estimation, thus, faster convergence with smaller variance.
(a) Information-based seeking1 (b) Field climbing [7]. (c) Following straight lines to the estimated location.   To further test the above conjecture that sensors moving freely and in the direction of improving Fisher information brings richer measurement and therefore enhances the estimation, we conduct the following three experiments: 1) The sensors move freely guided by our proposed algorithm; 2) The sensors are restricted to staying outside a radius of 3.0 from the source, performing projected gradient descent of our proposed loss function at the boundary; 3) The sensors do not move but only perform location estimation. The results are shown in Fig. 6. We see that a larger team corresponds to a faster decline in the expected estimation error with a smaller variance. Moreover, sensors moving freely as in Figure 6b results in a quicker decrease in estimation error with smaller variance.
3) Comparison of Different Information Metrics: Although we specify the loss function as the A-optimality criterion in our algorithm, one can, in principle, replace it with other alternatives. In the following experiments, we test our algorithm's performance with four different loss functions: i) Tr(FIM −1 ); ii) λ max (FIM −1 ); iii) − log det(FIM), and iv) Tr P where P is the posterior covariance of the EKF defined by [27, Equation (15)], P = (∇ q HR −1 ∇ q H + P −1 0 ) −1 , with P 0 being the current estimation covariance and R being the measurement noise parameter of the EKF.
All loss functions are tested with a stationary source and three freely-moving sensors for 100 repetitive trials. Each trial is initialized randomly in the same way as in the previous subsection. Fig. 7 shows the gradient descent trajectories for the tested metrics. Figure 8 shows that all the tested metrics except − log det(FIM) yield relatively good performance when using our algorithm, with Tr(FIM −1 ) and λ max (FIM −1 ) achieving a better balance in convergence and estimation than the covariance metric Tr(P ). These results confirm the generality of our algorithms.
Remark 5 (Rationale behind the metrics). Assume the FIM is a k ×k matrix. Since the FIM is always positive semi-definite, it has non-negative eigenvalues λ 1 , λ 2 , ..., λ k . The goal of minimizing CRLB is to make the matrix FIM −1 as close to the zero matrix as possible, which means 1 λ1 , ..., 1 λ k should all be as close to zero as possible. Therefore, Tr (FIM −1 ) = k i=1  Fig. 6: The time evolution of source location estimation error. In Fig. 6b, the variance of the 50 sensor curve after 10 steps is small but not zero.
1 λi is very close to zero while all others remain large. The covariance metric is essentially Tr((FIM + P −1 0 ) −1 ), which can be thought of as a regularized version of the first metric.
Remark 6 (Geometric Properties). Although the trajectories in Fig. 7 may look very different, they share some important common geometric properties that make sense intuitively. Specifically, they all encourage a separation among the sensors at some point and create an overall trend of approaching the source. We need to look at the explicit form of FIM to understand these behaviors. In the proof in Appendix B, we show that where m is the number of sensors, g i is the measurement function for the i'th sensor, r i is the distance between the i'th sensor and the source, andr i is the unit vector pointing from the source to the i'th sensor. Having FIM −1 close to zero means FIM itself is 'very positive definite'. And to achieve that, either 1)|g i (r i )| are large or 2)r ir i are 'very linearlyindependent', or both. Under the assumption that |g i (r i )| becomes larger as r i decreases, item 1) contributes to the overall trend of approaching the sensor, which seeks a stronger signal-to-noise ratio. Further, item 2) encourages the unit vectorsr i to point towards different directions, making the sensors separate from one another and creating different lineof-sights. The covariance is perhaps the one that stands out the most among the metrics as it exhibits some follow-straight-line behavior at the beginning and separates only when close to the source. This phenomenon can be explained via the dominance between P −1 0 and ∇HR −1 ∇H , and we include a detailed discussion in Appendix D.

4) Robustness to Measurement Modeling Error:
We next investigate whether our algorithm can function despite the error in the measurement model. We simulate source seeking with 10 mobile sensors and a stationary source, in which the measurement is generated by y = 1/r 2 +ν, ν i.i.d ∼ N (0, 0.01). We provide imperfect measurement models to the EKF in the form of h(r) = 1/rb = 1/r 2+∆b , with ∆b taking values in 0, ±0.1, ±0.5. We study the estimation error in two settings: 1) The sensors move freely using our algorithm; 2) The sensors are stationary. The results are shown in Fig. 9. The robustness of our algorithm is two-fold: In Fig. 9a, when compared to using the perfect measurement model (∆b = 0), our algorithm shows no significant degrade in estimation when ∆b is small. Therefore our algorithm maintains reasonable estimation quality even if the measurement model is imperfect. Besides, comparing Fig. 9a and Fig. 9b, our algorithm shows more robustness than stationary sensors, whose estimations tend to divergence when ∆b is at the value of 0.5.

B. The effectiveness of our distributed algorithm
In this subsection, we showcase the effectiveness of our algorithm when extended to the distributed setting. In Fig. 10-16, the x-axes are time steps, and the color bands show the standard deviation across 100 repetitive experiments. Specifically, our experiments compare the following groups: • Our distributed algorithm 2: consensus is used in estimation and gradient update. The sensors have a shared loss function. • Variation I: consensus is used in estimation but not in gradient update. The sensors have different local loss functions. • Variation II: each sensor gathers measurement information from its neighbors, but the estimation is done on its own. No consensus in either estimation or gradient update. Local EKF is used instead of consensus EKF. The sensors have different local loss functions. • Algorithm 1. The first three groups are implemented as distributed algorithms. The fourth is Algorithm 1 with a central controller, which serves as the control group. Detailed descriptions for the two variants I and II are provided in Appendix C. All sensors use the same measurement model as defined in (24). We simplify the sensor communication network so that the distributed algorithms use a static, undirected circulant network with each mobile sensor connected with two neighbors.
Although the neighbors of a sensor in physical locations could change over time, the communication neighborhood is fixed and does not vary in our experiments. We program the robots to communicate via a single-hop, static network, where each robot communicates with and only with a fixed set of neighbors on the communication network. This feature is also possible in our hardware implementation because our lab space is small enough that all sensors are within the maximal range of communication with one another. Using a fixed communication network preserves the effect of limited information in the distributed setting while allowing simpler, more interpretable results from the experiments. However, fixed communication is indeed a limiting assumption. Our algorithm can be revised to accommodate a more general, potentially time-varying network, such as the ad-hoc network where the robots can only communicate with teammates within a limited distance. Based on the previous studies on distributed consensus and optimization over time-varying communication graphs, e.g., [47], [49], [50], we expect our algorithm to achieve good performance as long as the network remains connected for a sufficient fraction of the time. This is left as our future work.
We repeat the experiments of seeking a static source. The sensor locations are initialized in the same way as in the previous subsection V-A. The source is still situated at (6,6). We use the fixed location (3,3) as the initial guess of the source location for the centralized control group. As for the first three distributed groups, each sensor uses an independent random location in [0, 3] × [0, 3] as the initial guess.
The main results are shown in Fig. 10 of Section V-B1 showing the time evolution of median sensor-source distance and estimation error. All variations converge to the source in around 100 ∼ 125 steps. Sections V-B2 and V-B3 further analyze the data from Section V-B1 to demonstrate the advantage of our distributed algorithm over Variations I and II in terms of full-team convergence and success/perfect finish rates.
1) Discussion on the main results: The curves in the first row of Figure 10 represent median i ||p i − q||, the median of source-sensor distances within the team. The second row shows the median estimation error median i ||q i − q||. All the estimation error plots in Figure 10 show the distributed variations have an overall worse estimation than the centralized algorithm. Also, note that when the number of sensors is small, as in Figure 10a, the centralized algorithm has a clear advantage over the distributed variations both in estimation and convergence to the source. These results exemplify one of the weaknesses of the distributed variations: each agent has access to much less measurement information than the centralized algorithm.
However, as the number of sensors increases, we can observe that distributed variations improve in both estimation and convergence to the source. As the number of sensors exceeds 10, Variation II and our distributed algorithm achieve comparable convergence to the source with the centralized algorithm. The centralized algorithm still reaches the source the fastest, in around 100 steps, while Variation II and our distributed algorithm take slightly longer than that. But interestingly, the distance to the source for the centralized algorithm seems to decrease slower than these two variations in the initial steps, and the reason is likely that the centralized algorithm tends to spread out more than the distributed variations due to the stronger coordination in sensor movements, as shown in the illustrative example in Figure 11.
It is also interesting that Variation I seems to be consistently worse than other distributed variations in convergence to the source. Recall that Variation I performs consensus in estima- The figures illustrate the general trend that the centralized algorithm tends to spread out more than the distributed variations in the initial steps. Also, Variation I spreads out more than Variation II and our distributed algorithm.
tion but not in gradient update, giving it more information than Variation II. We do observe from Figure 10 that in the early stages, its estimation error is comparable with our distributed algorithm and is lower than Variation II. However, as illustrated in Figure 11, Variation I tends to spread out more than other distributed variations, yet it does not coordinate the movement among the agents, meaning it is slow to converge to the source. These factors potentially explain why Variation I does not have an advantage over Variation II in terms of the median i ||p i − q|| metric.
The main results suggest that with a large number of sensors, even with the limited information for each sensor, as in Variation II, the distributed variations can still perform well. The sheer size of the distributed sensor team allows a high chance of some sensors getting close to the source and obtaining good measurement readings. Nevertheless, our distributed algorithm still outperforms Variations I and II due to the additional consensus in estimation and movement coordination.
Overall, our distributed algorithm converges to the source at a comparable rate as the centralized implementation. It also converges to the source slightly faster than Variations I/II and exhibits more improvement when using more sensors. Furthermore, our algorithm maintains a clear advantage in estimation accuracy over Variations I/II throughout the seeking process.
2) The convergence of the entire team: The advantage of our distributed algorithm over Variations I and II is more prominent in terms of the entire team's convergence to the source. Fig. 12 shows the performance of the variations in terms of max i ||p i −q|| using the same data as in Fig. 10, where ||p i − q|| is the distance between sensor i and the source. The convergence of max i ||p i −q|| to 0 thus characterizes the entire team's convergence to the source, a condition stronger than the one we proposed in Eq. (4). Figure 12 shows the max i ||p i −q|| of our distributed algorithm converges to the source almost as fast as the centralized algorithm when the number of sensors exceeds 10; meanwhile, as the number of sensors increases, the performance of Variation I gets slightly worse, and Variation II gets significantly worse. These results demonstrate that with consensus in estimation and motion planning, our distributed algorithm can ensure all sensors come to the source sufficiently fast. In contrast, without sufficient coordination via consensus, as in Variations I and II, some sensors could perform much worse than others.
3) The Success and Perfect Finish Rates: Fig. 13 shows the perfect finish/success rates of different variations of our algorithm. The figures are generated using the same data from Fig. 10. A robot is deemed to have reached the source if its distance from the source is below 0.2. A 'success' is defined in Eq. (4), corresponding to a trial where at least one sensor reaches the source at the final step. A 'perfect finish' is a trial where all sensors achieve that, which, like the max i ||p i − q|| metric, is a notion characterizing the entire team's convergence. The success and perfect finish rates are computed over 100 trials. The results show that all variations have close to 13 (a) 4 Sensors (b) 10 Sensors (c) 20 Sensors (d) 40 Sensors Fig. 12: Comparison of the entire team's convergence to the source for different variations, in terms of max i ||p i − q||. The figures are generated using the same data as in Fig. 10.
100% success rates, which demonstrates the effectiveness of our information-based source-seeking framework. Meanwhile, the perfect finish rates of the distributed algorithm are very close to one consistently over any number of sensors, but the perfect finish rates for Variation I and II drop significantly as the number of sensors increases. The perfect finish rate results again demonstrate the crucial role of consensus in estimation and motion planning in our distributed algorithm: It allows the entire team to reach the source, achieving a more challenging control goal than Eq. (4). 4) Discussion on non-monotonic estimation error: In Fig.  10, we observe that the estimation errors do not decrease monotonically. Indeed, classical results for Kalman filters show that monotonic convergence can be achieved, but that applies only to linear systems, where both the measurement and motion models are linear. In contrast, the estimation of EKF on a non-linear measurement model, the estimation error of which is generally not monotonic and could diverge given bad initialization. Although the criteria for the convergence of EKF is discussed in detail in [36], we are unaware of any literature that discusses the criterion for monotonicity of the estimation error for EKF.
We also observe in Fig. 10 the somewhat counter-intuitive behavior that source-sensor distances can decrease as the estimation errors increase. Fig. 14 unveils the reason for this curious phenomenon. It shows a common situation in the experiments, where the estimation from the EKF could move toward and overshoot the true source location but then comes back to the source as the sensors spread out and move closer to the source. During the overshoot, the estimation error gets larger, but the estimation still guides the sensors to move toward the direction of the source. That is why we observe that the estimation error increases while the sensor-source distance decreases during this period. We also note that in Fig. 14d for Variation II, the estimations of two sensors are in the direction of the source, while the estimations for the other two are not in the direction of the source. But because the two sensors with the 'right direction' come closer to the source, their measurements become better/more informative, which is shared with the other two sensors and 'saves' them from divergence. This phenomenon is common in experiments and explains why Variation II performs well in terms of the median distance metric in Fig. 10 even though it uses no consensus in estimation and motion planning.

C. Sensitivity to the Initial Guess
Our distributed algorithm is not only as effective in source seeking as the centralized implementation; it is also more robust in multiple ways. For example, the distributed algorithm is much less sensitive to the initial guesses of the source location. In the following experiments, we compare the sensitivity of the distributed and centralized implementation to initial guesses in terms of convergence and estimation error. Five sensors are used in all experiments. The locations of the sensors and the source are initialized in the same way as in the previous experiments. Random initial guessesq 0 are given to estimators in both the distributed and centralized implementations. It is defined bŷ In other words,q 0 is a random location in a box with side length D centered at q. Different estimators may have different q 0 , but the level of deviation D is kept the same among estimators for one experiment. Fig. 15 shows the distributed algorithm exhibit clear advantages over the centralized implementation in convergence and estimation when deviation D increases. The distributed algorithm prevails because using more estimators adds robustness. Although the centralized implementation has more data for estimation, it runs only one estimator. If the initial guess is poor, the estimator may fail. In contrast, sensors of the distributed algorithm form an ensemble of estimators. As long as the majority of the initial guesses are decent, individual estimations should improve after rounds of consensus.

D. Robustness to Communication Delay
Apart from being robust to errors in initial guesses, the distributed algorithm is more robust than the centralized implementation to communication delay. Note that communication delay typically grows with the incoming information rate [51, M/M/1 system, Section 4.1.2], which usually increases with the number of incoming connections. Therefore, the delay increases with the number of incoming connections. The consequence is that the delay in centralized implementation increases with the number of sensors used, while the delay in the distributed algorithm can remain low as the network expands as long as the incoming degree of each node stays unchanged.
We show that the distributed algorithm is more robust to communication delay in the following numerical experiments. We assume each incoming connection to a node brings a 0.5 time step delay in information passing. That is, if m sensors are used in the centralized algorithm, the centralized controller receives all information with a delay of m 2 time steps; if there are k neighbors connected to sensor j in the distributed algorithm, sensor j receives all information with a delay of k 2 time steps. The delay is modeled in this way since each sensor takes and sends out measurements to the neighbors at a roughly constant rate. Initialization is the same as in the previous experiments. In particular, each sensor in the distributed algorithm is connected to two neighbors regardless of the total number of sensors. Fig. 16 shows that when the delay is in effect, the advantage of the centralized implementation over the distributed algorithm dissipates as the number of sensors grows. In the centralized implementation, the benefit of more sensor data does not offset the negative effect of increasing delay, and the performance gradually degrades. On the other hand, distributed sensors benefit from the increasing global information shared in the network while the delay remains low. Therefore, the centralized implementation is eventually overtaken by the distributed algorithm.

E. Lab Implementations
We implement our distributed algorithm on Turtlebot3 robotic ground vehicles to seek a light source in a dark room. The light source is an LED lamp. Light sensors are installed on Turtlebots to measure the local strength of light. The robots communicate through an ad-hoc WiFi network, each being an individual node. We only allow each robot to communicate with two other pre-determined robots throughout the experiment as a reasonable approximation of a non-fullyconnected network. An indoor motion capture system captures the positions of the robots and the source.
The algorithm is implemented as a ROS 2 package in Python language [52]. The package is installed on the onboard computer of every robot. By running the package, every robot gathers its position from the motion capture system, communicates information required in estimation and waypoint planning, makes estimations, and decides its control actions. As a measurement model is needed to estimate the source location, we fit a function for each robot describing the relationship between their measurements and distances to the source before deploying them in source seeking. The function is defined by where y i , r i are the measured data from the ith mobile sensor, and k i > 0, b i > 0, C 0,i , C 1,i are the associated model parameters.
In the implementation, a robot is deemed to have reached the source if its distance from the source is below 0.5m. A separate program that runs independently from the source seeking algorithm continuously monitors whether any robot has reached the source using the motion capture data and notifies the robots if they make contact with the source. The accompanying video shows that the robots implementing our distributed algorithm consistently converge to the source. Remark 7. The hardware implementation is constrained by the lab resources, and it could be more comprehensive than the numerical tests. This is one limitation of the experiments which could be improved. We consider the hardware experiments to serve as a proof-of-concept showing the proposed algorithms are implementable on physical systems with satisfactory performance. We hope the extensive numerical tests and the revised hardware experiments together demonstrate how the algorithm works.

VI. CONCLUSION
This paper proposes a method for multi-robot source seeking which utilizes the Fisher Information associated with the estimated source location to direct the robots' movements. We show that improving the trace of the inverse of Fisher Information aids estimation while guiding the robots to converge to the source. The method is extended to the distributed setting where a central controller is absent, and the robots make individual decisions while communicating via a distributed network. We propose consensus schemes that make distributed estimation and gradient update effective. The algorithms are verified in numerical experiments and on physical robotic systems.
Our work leads to a few future directions. One is to combine our work with Bayesian optimization methods to learn the measurement model while the robots seek the source. Progress in this direction removes the requirement of a known measurement model, making the method more flexible. It is also imperative to develop the ability to navigate an environment with obstacles during source seeking, especially when deploying our method in complex environments, for example, during search-and-rescue missions.

APPENDIX A MOTION PLANNING
The motion planner M P can be viewed as a device that transforms the planned waypoints into low-level actuation of a mobile sensor: by applying (u i (1), . . . , u i (T )) = M P (p i (0), . . . ,p i (T )) to sensor i, the sensor will follow the trajectory of the waypointsp i (0), . . . ,p i (T ). The motion planner M P typically requires the knowledge of the sensors' motion dynamics to compute the control inputs, and any method that can fulfill this task can be a motion planner in Algorithm 1. In our implementation, the motion planner combines spline-based motion generation described in [53] and the Linear Quadratic Regulator (LQR). Regarding the planning horizon T choice, in the Gazebo simulations and hardware implementation, we set T = 20 to ensure stability in sensor movement and robustness to disturbances. Meanwhile, in numerical studies where the robot dynamics are not simulated, it suffices to set T = 1.

APPENDIX B PROOF OF PROPOSITIONS 1 AND 2
Recall thatr i is the directional unit vector pointing from q (the source) to p i (sensor i's position), i.e., We first derive expressions for FIM and the loss function L.
Lemma 1. Suppose individual measurements are isotropic as in Assumption 1. Then Moreover, where λ i (FIM) is the ith eigenvalue of FIM.
Proof. By noticing that r i = p i − q and using the chain rule of calculus, one can show that . Then from the definition of H in (2), we have and so Then (28)  With the formulas of L and FIM, and by further imposing the monotonicity property of h i in Assumption 1.2, we can show that minimizing L leads to reaching the source as stated in Proposition 1.
So by Assumption 1, |g i (r i )| monotonically increases as r i decreases, so we have min i r i tends to decrease as L(p, q) decreases, which completes the proof.
We prove the theorem by contradiction. Assume min i ||p i (t) − q|| does not converge to 0. Then there must be constant > 0 such that Together with the continuity of L(p, q) in p on {p : p i = q, ∀i}, the assumption leads to the conclusion that the infimum of L in the theorem should be replaced with the minimum, and is attainable atp. We now construct a contradiction to the assumption in Eq. 30 showingp cannot be the minimum of L. Consider the location vector b = [b 1 , b 2 , ..., b m ] defined as So that ||b i − q|| = 2 .
Note that b satisfies b i = q for all i. We also note that Note that α > 1 sincer i > /2 and |g i (r i )| is assumed strictly monotone in r i in the theorem. Therefore, <L(p, q) Therefore, the location vector b, which is closer to q than p, achieves a lower L value thanp. We have reached a contradiction, and Eq. (30) cannot hold. Sor i = 0 for some i and lim t→∞ min i ||p i (t) − q|| = 0.

APPENDIX C DISTRIBUTED ALGORITHM VARIANTS I AND II
The two variations are formally defined as follows: 1) Sensor j in variation I and II no longer share a loss function with all sensors. Instead, it uses the following local loss that may be different for other sensors.
L j (p j , q) = Tr ∇ q H j (p j , q) · ∇ q H j (p j , q) −1 (31) This local loss reflects the Fisher information about q contained in the neighborhood measurements as opposed to the entirety of sensor measurements. 2) Sensor j in variation II performs the local Kalman Filter update defined in (18) using neighborhood information p j , y j . No consensus is required. 3) Sensor j in variation I and II skips line 5-7 in Algorithm 2, and changes line 8-9 to be: •p j (0) ← p j • Calculate joint-waypointsp j (1), ...,p j (T ) p j (t + 1) ←p j (t) − α t M t ∇ pj L j (p j =p j (t),q j ) (32) • Extract sensor j's waypointsp j (1), ...,p j (T ) from joint-waypointsp j (1), ...,p j (T ) Note that no consensus is required to compute the gradient of L j . The remaining parts of the algorithm are identical to our algorithm and the variations. In particular, sensor j in Variation I still performs consensus EKF update. Figure 7(d) shows that the sensors under the covariance metric approach the source in straight lines and then spread out a little when close to the source. We can understand the initial straight-line behavior by analyzing the direction of the gradient when sensors are far away from the source. Recall that the covariance metric defined in Section V.A 3) is P = (∇HR −1 ∇H + P −1 0 ) −1 . We set R = I in Figure 7, so the covariance metric is effectively Tr(P ) = Tr((∇H∇H + P −1 0 ) −1 ). Through matrix calculus, we can show that the covariance metric's gradient is given by the following.

TRAJECTORY OF THE COVARIANCE METRIC
G j :=∇ pj Tr((∇H∇H + P −1 0 ) −1 ) = − 2∇ pj A j (∇H∇H + P −1 0 ) −2 A j . Where A j := ∇ q g j (||p j −q||), and A j , ∇ pj A j can be shown to be A j = −g j (r j )r j ∇ pj A j = −g j (r j )r jr j − g j (r j ) r jt jt j Wherer j = pj −q ||pj −q|| is the unit vector pointing from q to p j , andt j is the unit vector orthogonal tor j .
Since we assume the individual measurement function |g j (r j )| gets smaller as r j increases, P −1 0 dominates ∇H∇H when r j are large(sensors are far away from the source). We note that we set P 0 = I in Figure 7(d). Therefore, when sensors are far away from the source, the gradient for sensor j is approximately Also, since the g j we consider are in the form of g j (r j ) = k j (r j − C 1,j ) −bj + C 0,j , its first and second order derivatives satisfies g j (r j ) ≤ 0 and g j (r j ) ≥ 0. Therefore, the descent direction −G j is approximately aligned with the direction of −r j when the sensors are far away from the source, causing the straight-line behavior.
When sensors are close to the source, analyzing the gradient direction is less trivial, but we do note that the covariance metric is roughly equivalent to Tr(F IM −1 ) in this condition since the ∇H∇H term dominates P −1 0 . We argued in Remark 6, that minimizing the Tr(F IM −1 ) amounts to approaching the source and creating separation between the sensors simultaneously, which explains why the sensors spread out a little at the end. Fig. 17 provides further numerical evidence showing the dominance between P −1 0 and ∇H∇H terms influence the sensor movements. It shows that 1) When P −1 0 is large, it dominates the covariance metric and gives rise to straight-line behavior; 2) On the other hand, when P −1 0 becomes smaller, the ∇H∇H term starts to dominate and the "bump" in the trajectory becomes increasingly larger. 3) When P −1 0 becomes extremely small, the trajectory shows the behavior of encouraging separation among the sensors from the very start, which matches the behavior of the Tr(F IM −1 ) metric, which is unsurprising since F IM = ∇H∇H . on the size of the "bump" of the gradient-descent trajectory of the covariance loss function. The matrix R is fixed to be 1.0I for all figures above.