Franc ¸ois Bailly

Justin Carpentier

Mehdi Benallegue

Bruno Watier

Philippe Souères

Estimating the Center of Mass and the Angular Momentum Derivative for Legged Locomotion -A recursive approach

Keywords:

Estimating the center of mass position and the angular momentum derivative of legged systems is essential for both controlling legged robots and analyzing human motion. In this paper, a novel recursive approach to concurrently and accurately estimate these two quantities together is introduced. The proposed method employs kinetic and kinematic measurements from classic sensors available in robotics and biomechanics, to effectively exploits the accuracy of each measurement in the spectral domain. The soundness of the proposed approach is first validated on a simulated humanoid robot, where ground truth data is available, against an Extend Kalman Filter. The results demonstrate that the proposed method reduces the estimation error on the center of mass position with regard to kinematic estimation alone, whereas at the same time, it provides an accurate estimation of the derivative of angular momentum. Finally, the effectiveness of the proposed method is illustrated on real measurements, obtained from walking experiments with the HRP-2 humanoid robot.

I. INTRODUCTION

T He precise estimation of both the center of mass (CoM) position and the angular momentum derivative is an important issue for legged robots control [START_REF] Orin | Centroidal dynamics of a humanoid robot[END_REF], [START_REF] Hofmann | Exploiting angular momentum to enhance bipedal center-of-mass control[END_REF]. It also has various applications in humans for motion analysis in biomechanics or for disorders diagnosis in medicine. The main challenge of this estimation problem lies in the fact that these two quantities are not directly measurable, and depend at the same time on fixed parameters (e.g. mass distribution of the system, dimensions of the limbs, etc.) as well as on varying quantities (e.g. joints configurations).

A. State of the art

In order to estimate the CoM position, several methods exploiting different sources of information are commonly applied [START_REF] Lafond | Comparison of three methods to estimate the center of mass during balance assessment[END_REF], [START_REF] Stephens | State estimation for force-controlled humanoid balance using simple models in the presence of modeling error[END_REF], [START_REF] Bae | Novel state estimation framework for humanoid robot[END_REF]. They can be grouped into two broad categories, namely the kinematic methods which only exploit kinematic measurements and the kinetic methods which estimate the CoM position from kinetic signals. In the following, estimation methods from both robotics and Fig. 1: Illustration of the measurement apparatus. The several physical quantities involved in the estimation framework are displayed, as well as a simplified sketch of the estimation algorithm.

biomechanics are reviewed, showing that the goals and means of the two disciplines are very similar on this topic.

Kinematic methods are certainly the most popular ones. They assume that the system is composed of rigid links and articulations, which is reasonable for non-soft robots, but more questionable in the context of humans [START_REF] De Rosario | Kinematic description of soft tissue artifacts: quantifying rigid versus deformation components and their relation with bone motion[END_REF]. The masses and joint angles required for the kinematic estimation of the CoM can be either obtained from CAD (Computer-Aided Design) data and encoders measurements in robotics, whereas they can be obtained from anthropometric tables coupled with motion capture (MoCap) systems in biomechanics [START_REF] Leva | Adjustments to zatsiorsky-seluyanov's segment inertia parameters[END_REF]. Then, the CoM position is function of a weighted sum accounting for both the mass distribution and the position of each segment's CoM. It is worth noticing that kinematic methods are subject to noises of various types, mostly due to uncertainties about the mass distribution of the system (discrepancy between the CAD model and the real robot, statistical anthropometric tables in humans), and to soft tissue artifacts [START_REF] De Rosario | Kinematic description of soft tissue artifacts: quantifying rigid versus deformation components and their relation with bone motion[END_REF].

As for kinetic methods, the estimation of the CoM position is achieved by only measuring either the so-called Zero Moment Point (ZMP) [START_REF] Sardain | Forces acting on a biped robot. center of pressure-zero moment point[END_REF] or directly from the contact forces measurements [START_REF] Lafond | Comparison of three methods to estimate the center of mass during balance assessment[END_REF], [START_REF] Rotella | Humanoid momentum estimation using sensed contact wrenches[END_REF]. For instance, in humanoid robotics, the robot is often modeled as a linear inverted pendulum. Under this assumption, the CoM position can be directly estimated from the mathematical relationship linking the ZMP to the CoM [START_REF] Kajita | The 3d linear inverted pendulum mode: A simple modeling for a biped walking pattern generation[END_REF], [START_REF] Morasso | Computing the com from the cop in postural sway movements[END_REF]. These approaches then assume that the CoM trajectory can be assimilated to the oscillations of a linear inverted pendulum, and the acceptable range of motion of such methods is mostly limited to a neighborhood of the balance configuration. In [START_REF] Caron | Estimating the centre of gravity of the body on the basis of the centre of pressure in standing posture[END_REF], it was proposed to filter the ZMP trajectory with a low-pass filter to estimate the horizontal components of the CoM. The accuracy of this approach was limited to slow movements. Finally, by assuming a constantheight CoM, performing a double integration of external forces together with an estimation of the angular momentum derivative expressed at the CoM, [START_REF] Shimba | An estimation of center of gravity from force platform data[END_REF] and [START_REF] Levin | An iterative model for estimation of the trajectory of center of gravity from bilateral reactive force measurements in standing sway[END_REF] have provided another way to estimate the horizontal components of the CoM. Some recent works have proposed to remove simplifying assumptions and to fuse kinetic and kinematic data through Extended Kalman Filtering (EKF) [START_REF] Rotella | Humanoid momentum estimation using sensed contact wrenches[END_REF], [START_REF] Piperakis | Nonlinear state estimation for humanoid robot walking[END_REF]. Nethertheless, an important difficulty when estimating the CoM position from forces/moments data, comes from the fact that this quantity is only partially observable [START_REF] Carpentier | Centerof-mass estimation for a polyarticulated system in contact-A spectral approach[END_REF]. This means that at each time instant not all the components of the CoM can be observed, and the full observability can only be recovered by exploiting the correlation between the past sensor measurements, which requires to save these data in memory.

Concerning the estimation of the derivative of angular momentum, to the best of our knowledge, only kinematic estimation has been proposed in the literature [START_REF] Dapena | A method to determine the angular momentum of a human body about three orthogonal axes passing through its center of gravity[END_REF], [START_REF] Bruijn | Coordination of leg swing, thorax rotations, and pelvis rotations during gait: the organisation of total body angular momentum[END_REF]. Such approaches are affected by the same limitations as the ones encountered in estimating the CoM position, regarding mass distribution errors, with additional uncertainties coming from the numerical derivation of segments poses. Kinetic data can also be used as suggested in [START_REF] Herr | Angular momentum in human walking[END_REF], by exploiting the relation between the angular momentum and the external torques applied to the system. It should be noticed, however, that in most cases, the derivative of angular momentum is neglected for the sake of simplicity. For instance, this hypothesis was made in a previous work [START_REF] Carpentier | Centerof-mass estimation for a polyarticulated system in contact-A spectral approach[END_REF], where the authors relied on a complementary filter in the spectral domain, to merge both the kinematic and kinetic data in order to exclusively estimate the CoM position, assuming that the angular momentum derivative is negligible. Nethertheless, as it has been highlighted in [START_REF] Wieber | Holonomy and nonholonomy in the dynamics of articulated motion[END_REF], these time variations of the angular momentum are present and necessary in biped locomotion, making this assumption very restrictive.

B. Contributions

In this paper, we introduce a generic framework to simultaneously and accurately estimate both the CoM position and the derivative of angular momentum by exploiting the mathematical properties that link these two quantities (see Fig. 1 for a schematic overview of the proposed method). Compared to the previous approaches of the state of the art, the contributions of this paper are threefold:

(i) We introduce a recursive algorithm in order to improve, over the iterations, the estimated CoM position from the estimated angular momentum derivative and vice-versa; (ii) Kinematic and kinetic data are complementary fused together for estimating the CoM (as in [START_REF] Carpentier | Centerof-mass estimation for a polyarticulated system in contact-A spectral approach[END_REF]) and to the best of our knowledge, for the first time the derivative of the angular momentum; (iii) The proposed method does not rely on any simplifying assumptions such as neglecting the derivative of the angular momentum (we provide an estimation for this Fig. 2: Schematic representation of the spectral accuracy of the different input signals for the estimation of c and Lc . quantity), or assuming that the center of mass lies on a horizontal plane, therefore, our formulation is generic and allows one to work with any number of contacts, in any multi-contact scenario which are now common in robotics. The consistency of the proposed approach is first validated on a simulated humanoid, by demonstrating the accuracy of our estimated quantities against ground truth data and against a state-of-the-art algorithm (EKF). Then, in order to illustrate its efficiency in experimental conditions, our approach is applied to a real dataset by analyzing a walking motion performed by the HRP-2 robot.

II. ESTIMATION ALGORITHM

In this section, we first recall how the position of the CoM (denoted by c) and the angular momentum derivative expressed at the CoM (denoted by Lc ) are linked together through the so-called Euler's equation of motion. Then, we identify the different sources of information which contribute to the estimation of Lc and c, and we explain how to fuse them by exploiting their level of accuracy in the spectral domain.

A. Link between the CoM and the angular momentum

Let m be the mass of the system, and c, ċ and c the CoM position, velocity and acceleration respectively. Let f and τ 0 be the resulting forces and moment exerted on the system by the environment and expressed at the origin 0 of the world. They constitute the external contact wrench that can be translated at any other point of the space using the Varignon formula. According to the Euler's equations of motion, the angular momentum derivative at the CoM is equal to the moment of the external forces expressed at the CoM, yielding:

Lc = τ 0 + f × c (1)
Eq. ( 1) reveals the physical coupling between the two quantities we want to estimate, namely the CoM position c and the angular momentum time derivative Lc .

B. Input measurements

Our estimation framework assumes that both the total mass m of the system and the gravitational acceleration g = (0, 0, -g) are known, with g the standard gravity value. These are rather reasonable assumptions since the mass is easily measurable from force sensors in stationary position and the gravity is precisely known. Using joints encoders or a MoCap Recursive estimation

+ + L COMPLEMENTARY FILTER c COMPLEMENTARY FILTER
Angular momentum calculation, Eq. ( 7) Eq. ( 8) Eq. (6b)

Measurements

Fig. 3: Flow chart of the recursive complementary estimation framework, with s the Laplace variable. Dotted lines represent the update step of the algorithm. Complementary filtered components are summed in order to output Lc est and c est , the estimates of Lc and c respectively. system and 6-D force measurement units, one can measure the kinematic and kinetic data related to the physical quantities of interest:

• the resulting wrench of contact (f , τ 0 ) expressed at the origin, via the force sensors (embedded sensors, force plates in biomechanics, etc.); • the position of the CoM, denoted by c kinematic , directly deduced from the kinematics or the MoCap system associated with the mass distribution of the system; • the angular momentum at the CoM, denoted by L c kinematic , which can be obtained from the angular velocities and the inertias of the different links.

C. Multi-source fusion for estimating the CoM position

This section explains how to compute and process the three different input signals that carry information about the position of the CoM. The spectral accuracy of each signal is illustrated in Fig. 2.

1) From kinetic data: The second Newton law states that:

c = f m + g. (2)
By integrating twice the right hand side of this equation, one can obtain c force , the estimation of the center of mass from the contact forces, thanks to the information provided by the force sensors. However without the knowledge of integration constants (initial position and initial velocity), this method leads to a quadratic drift visible in the low and medium frequencies [START_REF] Carpentier | Centerof-mass estimation for a polyarticulated system in contact-A spectral approach[END_REF].

2) From kinematic data: Kinematic computations using either a model of the robot or reflective markers combined with anthropometric tables are usually used to estimate the CoM position. This source of data (c kinematic ) suffers mainly from low-frequency biases due to modeling errors of the mass distribution. It can also be altered by the high frequency sensor noise due to the MoCap technology [START_REF] Winter | Biomechanics and motor control of human movement[END_REF]. Then, the error between this signal and the real position of the CoM mostly lies in the low and high frequency domains.

3) From the central axis of the contact wrench: Considering the contact moment field, there exists one unique axis called the central axis of the contact wrench and denoted by ∆ [START_REF] Dimentberg | The screw calculus and its applications in mechanics[END_REF], [START_REF] Bailly | A mechanical descriptor of human locomotion and its application to multicontact walking in humanoids[END_REF], where the moment of contact forces is collinear to f . The projection c ∆ of the CoM c onto ∆ is given by:

c ∆ = f × τ 0 + (c • f )f f 2 , (3)
where • and × denote the dot and cross products, respectively. This projection can be computed from the force and moment signals provided by the force/torque sensors and the current estimate of the CoM, denoted by c est as follows:

c ∆ = f × τ 0 f 2 + (c est • n)n, (4)
where n = f f is the normalized direction of f . Knowing an estimate of Lc , denoted by Lc est , and f , and using the properties of the vector triple product, the difference between c and c ∆ can be computed as:

c -c ∆ = Lc est × f f 2 . ( 5
)
This leads to a third estimate of the CoM position that will be denoted by c axis and expressed as:

c axis = c ∆ + Lc est × f f 2 , ( 6a
)
c axis = f × τ 0 f 2 + (c est • n)n + Lc est × f f 2 . (6b)
4) Complementary filter of the CoM: The theory of complementary filtering provides a simple and effective way to fuse different sources of information in the presence of noise. The idea is to filter each signal in the spectral domain in order to preserve it in the frequency bandwidth where the signal has the best accuracy. To this end, complementary filters have to be designed, requiring that the sum of their transfer functions is equal to one, maintaining the energy of the signal with zero phase modification [START_REF] Higgins | A comparison of complementary and kalman filtering[END_REF]. In theoretical conditions, where the spectral distribution of the noise is perfectly known, this simple CoM

ω 2 lc s 2 + 2ζω lc s + ω 2 lc s 2 s 2 + 2ζω hc s + ω 2 hc 2ζω hc s 3 + [ω 2 hc -ω 2 lc + 4ζ 2 ω hc ω lc ]s 2 + 2ζω lc ω 2 hc s s 4 + 2ζ(ω lc + ω hc )s 3 + [ω 2 lc + ω 2 hc + 4ζ 2 ω lc ω hc ]s 2 + 2ζω lc ω hc (ω hc + ω lc )s + ω 2 lc ω 2 hc Lc ω 2 l L s 2 + 2ζω l L s + ω 2 l L s 2 + 2ζω l L s s 2 + 2ζω l L s + ω 2 l L - TABLE I:
Transfer functions of the complementary filters associated to the center of mass (first row) or the time derivative of the angular momentum (second row). s is the Laplace variable, ω lc and ω hc denote the low-cut and high-cut frequencies of the CoM filter, ω l L denotes the cutoff frequency of the Lc filter. ζ is the damping ratio of the filters.

k ← 0 cest[k] ← c kinematic Lc est [k] ← Lc kinematic while (||∆cest|| > and ||∆ Lc est || > ) do Lc force [k + 1], c axis [k + 1] ←update(cest[k], Lc est [k]) cest[k + 1] ←comFilter(c force , c axis [k + 1], c kinematic ) Lc est [k + 1] ←amFilter( Lc kinematic , Lc force [k + 1]) ∆ cest ← cest[k + 1] -cest[k] ∆ Lc est ← Lc est [k + 1] -Lc est [k] k ← k + 1
Fig. 4: Recursive estimation algorithm. The update() function implements the computation of Eqs. (6b) and ( 8). The comFilter() and amFilter() functions computes the complementary filtering of c est and Lc est respectively. is a tolerance threshold to be chosen according to the desired convergence accuracy. strategy is the optimal way to combine the signals without loss of information.

Similarly to [START_REF] Carpentier | Centerof-mass estimation for a polyarticulated system in contact-A spectral approach[END_REF], we rely on a complementary filtering approach to fuse the three signals, c force , c kinematic and c axis , in order to retrieve an accurate estimate of the CoM position. Thus, a complementary filters consisting of a highpass, a band-pass and a low-pass, as evident in Fig. 3 and Tab. I, is designed. First, we design the low-pass and high-pass transfer functions, with ω lc and ω hc their cutoff frequencies. The band-pass transfer function is then analytically deduced so that three transfer functions sum to 1.

D. Multi-source estimation of the angular momentum derivative

In this section, we explain how to compute and process the two input signals for estimating Lc from kinematic and kinetic signals.

1) From kinematic data: When the angular and linear velocities of the segments and their inertia matrices are known, then the angular momentum of the system can be readily computed as:

L c kinematic = i I i ω i + -→ cc i × m i v i , (7)
where I i , ω i , c i , m i and v i denote the inertia matrix, angular velocity, CoM position, mass and relative velocity of the ith segment respectively. The time derivative of this quantity provides an estimation of Lc from this first source of information at the price of amplifying high frequency noise (see Fig. 2). Indeed, the gain of the derivative operator increases linearly with frequency, and consequently, high frequencies (where the signal to noise ratio is usually bad, due to sensor noise) are greatly amplified. Therefore, the output of the derivative operator needs to be low-pass filtered as sketched in Fig. 3.

2) From the moment field properties: The second source of information for the estimation of Lc , lies in Eq. ( 1):

Lc force = τ 0 + f × c, (8)
and is provided by force sensors measurements with an estimation of c.

3) Complementary filter for the angular momentum variation quantity: In the same way as before, complementary filtering is used to add both contributions to the estimation of Lc . In the Laplace domain, two complementary filters are designed (high-pass and low-pass, see Fig. 3 and Tab. I for details) such that their transfer functions sum to 1. To our knowledge, it is the first time that the estimation of Lc is performed using both kinematic and kinetic data.

E. Recursive estimation of the CoM position and of the angular momentum derivative

The proposed recursive approach is motivated by the fact that the estimation of Lc force depends on an estimation of c (see Eq. ( 8)) and that the estimation of c axis depends on an estimation of Lc (see Eq. (6b)). Hence, the estimation accuracy of c can be increased by refining the estimation of Lc , and vice versa. Fig. 4 depicts the recursive algorithm that lies at the core of the method. Finally, Fig. 3 illustrates a detailed diagram of the proposed estimation framework.

III. EXPERIMENTAL VALIDATION OF THE ESTIMATION

FRAMEWORK

In this section, the efficiency of the estimation framework is first demonstrated on simulated data. In addition, an EKF-based estimator was implemented, as in [START_REF] Rotella | Humanoid momentum estimation using sensed contact wrenches[END_REF], in order to compare both approaches. Then experimental data are processed to illustrate the performance of the proposed method on a real application scenario. In both cases, the estimation is performed recursively as sketched in Fig. 3. A. Generation of ground truth and noisy measurements in simulation Initially, a dynamically consistent walking motion has been simulated for a humanoid robot. Then, by adding artificial noise to the several quantities involved in the estimation framework (with noise values representative of classic robotics setups), we compare our center of mass and angular momentum derivative estimation to ground truth data (available in simulation).

The simulated avatar consisted of a humanoid model endowed with 36 degrees of freedom. The simulated motion corresponds to a walking motion on horizontal ground at 0.5 m/s using the dynamically consistent approach introduced in [START_REF] Carpentier | A versatile and efficient pattern generator for generalized legged locomotion[END_REF]. Ground truth kinetic and kinematics data were computed at 200 Hz using the modeling software Pinocchio, which includes rigid body dynamics algorithms for polyarticulated systems in contact [START_REF] Carpentier | The Pinocchio C++ library -A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives[END_REF]. Band-limited centered white noise was added to each signal, in accordance with our analysis of the noise spectral density. The low, medium and high frequency domains were delimited from 0 to 4 Hz, 16 to 24 Hz and 92 to 100 Hz respectively, in order to be able to discriminate the effect of each noise and to limit the slopes and the order of the filters. If needed, the final user would only need to adjust the cutting frequencies inside the complementary filters, which depend on the sampling frequency as well as on the nature of the recorded motion. The standard deviations of the noise were 1 N and 1 Nm in the low and medium frequency domains for the force and moment measurements respectively, 0.1 m and 0.1 Nm in the high frequency domain for c kinematic and Lc kinematic respectively, in order to model the noise of the 6-axis sensors. Inertial data (mass, CoM position and inertia matrix) of each segment were randomly biased in order to simulate the error due to mass distribution uncertainties (multiplicative random Fig. 7: Simulation. Evolution of the error 2-norm integrated over the whole trajectory for c kinematic , c axis and c est . This plot is the result of the algorithm in Fig. 4. The integrated errors are normalized with regard to the error of c kinematic , which is the initialization of the algorithm and is not updated throughout the iterations (and therefore constantly 1.) value following a 1-centered normal distribution with standard deviation of 1). The cut-off frequencies were 1 Hz and 25 Hz for the low-pass and high-pass second order filters of c est respectively. The band-pass filter was deduced by complementarity (see Tab. I for detailed transfer functions). The cut-off frequency for the low pass filter of Lc est was 2 Hz, the high-pass filter was similarly deduced by complementarity. For both estimates, the damping ratio of the transfer functions was set to its critically damped value (ζ = 1).

The estimation procedure was initialized with the value of c kinematic before filtering, because it is the best guess one can make before adding further data. The recursion was stopped when the last estimates of Lc and c were close (10 -3 ) to the current ones (see Fig. 4 for a detailed presentation of the recursive algorithm). Fig. 7 depicts the error of c kinematic , c axis and c est integrated over the whole trajectory and normalized with regard to the error of c kinematic . This result shows how the complementary filter improves the initial guess c kinematic by reducing the norm of the error with regard to the ground-truth values (which are exactly known in simulation). In Fig. 5, the time residuals of the error of c kinematic , c axis and c est are displayed. Both the 1-norm and the mean of these residuals are displayed in order to demonstrate how our estimation method improves the kinematic estimation of the CoM, and that, most of the time, it is more accurate than the EKF approach. It should be noted that the reason why c axis is that good is because it is updated at each iteration, and it therefore improves in the same way as c est .

In more details, Fig. 6 shows the X, Y and Z components of c and Lc for the different estimates. Ground truth data computed without noise during the simulation are also displayed. On the top row of Fig. 6, c kinematic is displayed, to exhibit the improvement provided by our method. A low frequency bias between c kinematic and the ground truth value is noticeable, in particular on the Y and Z components because of the scale of the plots. This bias comes from the error on the inertial parameters that we purposely introduced in the simulated data. Fig. 6 also shows that the final estimation Fig. 8: Experimentation. Snapshots of the HRP-2 robot, walking on horizontal ground in the experimental conditions. c est is better than the initial guess, as corroborated by the error analysis in Fig. 7. For instance, on the Z axis, the average absolute error of the kinematic estimate is of 3.6 cm, whereas it is reduced to 1.1 cm for the final estimate. On the Y axis the average absolute error is reduced from 2.1 cm to 2.9 mm by the presented framework. On the X axis, which is the main direction of locomotion, no significant improvement is measured: the average absolute error remains close to 6 mm. Fig. 6 bottom displays the estimation of Lc , and the result follows the dynamics of the ground truth trajectory. Nethertheless, some overshooting can be noticed in all axes. As a preliminary conclusion, these simulated results demonstrate the efficiency of the proposed method on purposefully noised data. Indeed, they show that recursively taking Lc into account not only improves the estimation of c but also provides a good estimation for this quantity which is usually neglected for the sake of simplicity. However, one should notice that neglecting Lc in the context of robot locomotion seems to be a rather harsh assumption, given the values that it takes, as also emphasised in [START_REF] Wieber | Holonomy and nonholonomy in the dynamics of articulated motion[END_REF].

B. Application to real robot experimental data

To complete the simulation work presented above, the proposed method is employed on a real dataset to demonstrate its effectiveness. To this end, the estimation framework was applied on real kinematic and kinetic data obtained during robot walking experiments, in order to retrieve the estimates of c and Lc .

Data from the HRP-2 robot were used for this application (see Fig. 8). The dataset was taken from [START_REF] Stasse | Benchmarking the hrp-2 humanoid robot during locomotion[END_REF]. In this experience, the robot was programmed to walk on a flat terrain without feedback on its pose in the room. The values of the encoders and of the 6-axis forces and torques sensors located in the ankles of the robot were collected at 200 Hz. The kinematic estimation of the CoM was performed by implementing an odometry algorithm which assumed that, when a foot is in contact, it can not slip. The contact detection of each foot was performed using the collected force data, and by setting a threshold of 20 N on the vertical component of the force. Classic rigid body algorithms implemented in Pinocchio [START_REF] Carpentier | The Pinocchio C++ library -A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives[END_REF] were used to compute the position of the CoM and the value of Lc . Fig. 9 displays a variety of curves related to the estimation of c and Lc for this experimental walking on horizontal ground. As there are no ground truth values available to assess the quality of the recursive estimation with regard to kinematic estimation alone, only a comparative and quantitative analysis is relevant. The order of magnitude of the differences between c est and c kinematic are consistent with the improvement towards the ground-truth values which were observed in simulation. This result suggests two remarks: the levels of noise introduced in simulation were realistic and the differences observed on real data after recursion correspond to an improvement in the estimate. Interestingly, one can notice that both the kinematic, the EKF and the recursive estimates corroborate that the robot tended to drift towards the left of its heading direction (≈ 20cm after 3m of walking), a fact which was observed during the experiments. However, the recursive estimate and the EKF one show wider oscillations in the frontal plane (Y axis) than kinematic estimation. This can be due to the fact that the non perfectly planar contacts of the soles with the floor cannot be detected by kinematic estimation, whereas they can be captured by kinetic data. Another interesting results lies in the correction brought by the recursive estimation of the CoM, on the vertical axis. Both the kinematic and EKF estimation lead to an altitude of the CoM which drifts up to ∼ 5 cm, while it is reduced to 2 cm by the recursive estimation. Fig. 9 bottom illustrates the estimation of Lc , which is consistent with the dynamic of the motion and the weight of the robot. The angular momentum derivative profiles obtain with the reccursive method appear to be more regular and less noisy than the one computed with the EKF. The values tend to be a bit smaller than the ones obtained in simulation, this is due to the fact that the motion was performed slower during the experiments.

C. Portability of the method to human motion analysis

The authors would like to emphasize the fact that the presented method can be directly applied to human motion analysis. Indeed, in this context, the kinematic estimates of c and Lc are obtained using 3D point clouds of reflective markers for computing the pose of the segments, along with anthropometric tables that contain regression equations to predict segments lengths and mass distributions. For such cases of application, the kinetic data would typically come from 6-axis forces sensors such as force plates. For processing human experimental data, the only parameters that need to be tuned are the cutoff frequencies of the complementary filters. Though ground truth on kinematic and inertial parameters can hardly be provided when it comes to human body, this method could be used for benchmarking existing anthropometric tables and other estimation methods. One promising outcome of using such an estimation framework in biomechanics would be to correct the error brought by the regression equations (computed in average over a population) thanks to the exploitation of the kinetic data (specific to the subject under study). Preliminary work has been done in this direction on a horizontal walking motion dataset [START_REF] Bailly | Recursive estimation of the human body's center of mass and angular momentum derivative[END_REF], demonstrating the effectiveness of the method for non-rigid locomotion systems, and further work will be achieved on more dynamical motions, involving multi-contact with the upper limbs.

IV. CONCLUSION

In this paper, a novel recursive filter was introduced, to simultaneously estimate the CoM position and the derivative of angular momentum expressed at the CoM. By complementary filtering kinematic and kinetic data, this method is able to get the best of both worlds. Our estimation framework does not rely on any simplifying hypothesis, such as neglecting the angular momentum derivative [START_REF] Carpentier | Centerof-mass estimation for a polyarticulated system in contact-A spectral approach[END_REF], or assuming that contact surfaces are parallel [START_REF] Levin | An iterative model for estimation of the trajectory of center of gravity from bilateral reactive force measurements in standing sway[END_REF], [START_REF] Morasso | Computing the com from the cop in postural sway movements[END_REF]. Thus, the method is generic to any in-contact locomotion scenario, and the estimation is improved with regard to former approaches, because the mathematical equations that are used are exact and exploit the coupling between c and Lc .

On a simulated walking robot, the results showed that the proposed framework enabled us to reduce the error between the estimation of the CoM position and its real value, when compared to kinematic or EKF estimation. This improvement was a direct consequence of the proposed recursive algorithm exploiting the spectral properties of each measurement in a complementary filter. More precisely, this improvement was mainly due to a low frequency bias correction brought by a better estimation of the angular momentum derivative, exploited in the geometrical link between the CoM position and the central axis of the contact wrench. Simulation results demonstrated that the estimation of the derivative of the angular momentum was not perfect, and this is one explanation to the small bias that remained on the estimation of the CoM position after convergence of the algorithm, in addition to its non-observability on the Z axis. The overall estimation could be enhanced by improving the estimation of the angular momentum derivative which is left as a future work.

In order to demonstrate the applicability of the method in experimental conditions, the estimation algorithm was also applied to real data acquired from the HRP-2 robot during a walking task on horizontal ground. The results suggest that our method is more accurate than kinematic or EKF estimation.

We would like to emphasize that estimating both the center of mass and derivative of angular momentum is a good step towards centroidal state feedback for online motion generation algorithm in legged robotics [START_REF] Carpentier | A versatile and efficient pattern generator for generalized legged locomotion[END_REF]. Taking into account the derivative of angular momentum in these motion generation and control pipelines is a requirement for achieving stateof-the-art dynamical motion. To this end, it is essential to contribute with estimation schemes that accurately and efficiently estimate the latter quantities. In this perspective, as a future work, we plan to integrate this filter in the motion control pipeline of the Gepetto team [START_REF] Carpentier | Multicontact locomotion of legged robots in complex environments-the loco3d project[END_REF], in order to get an online feedback of the centroidal state of our humanoid robot.

Filters

-0. 10 -Fig. 5 :Fig. 6 :

 1056 Fig. 5: Simulation. Time residuals of c kinematic , c ekf , c axis and c est after convergence of the recursive algorithm, for a horizontal walking motion of the HRP-2 robot in simulation. The 1-norm and the mean of the residuals are displayed.X axis Y axis Z axis

Fig. 9 :

 9 Fig. 9: Experimentation. Estimation of c and Lc for a horizontal walking motion of the HRP-2 robot. Top row c est , c kinematic and c ekf . Bottom row : Lc est and Lc ekf .


Nassim Athmouni

Hatem Baloudi

AND Mondher Damak

Marwa Ennaceur

THE DISCRETE LAPLACIAN ASSOCIATED TO EDGES WITH APPLICATION TO MAGNETIC ADJACENCY MATRIX

Keywords: 2010 Mathematics Subject Classification. 05C12, 05C63, 47B25, 81Q35 spectral graph theory, infinite graphs, line graph, adjacency matrix

We consider the notion of χ-completness of the locally finite graph and we approach this geometric hypothesis for the weighted magnetic graph. Moreover, we establish a link between the magnetic adjacency matrix on line graph and the magnetic discrete Laplacian on 1-forms.

Introduction

Spectral graph theory represents an active area of research. In the past few year, the question of essential self-adjointness of the discrete magnetic Laplacian on magnetic graph has attracted a lot of interests. For the reference to the literature on this topic, see [START_REF] Golénia | Hardy inequality and asymptotic eigenvalue distribution for discrete Laplacians[END_REF][START_REF] Milatovic | Essential self-adjointness of magnetic Schrödinger operators on locally finite graphs[END_REF][START_REF] Schmidt | On the existence and uniqueness of self-adjoint realizations of discrete (magnetic) Schroödinger operators[END_REF][START_REF] Milatovic | A Sears-type self-adjointness result for discrete magnetic Schrödinger operators[END_REF] and reference therein. In [START_REF] Anné | The Gauss-Bonnet operator of an infinite graph[END_REF], in the case of non-magnetic graph, the authors have addressed this question using the notion of the χ-completeness. This notion is related to [START_REF] De Verdière | Essential self-adjointness for combinatorial Schrödinger operators II: geometrically non complete graphs[END_REF][START_REF] De Verdière | Essential self-adjointness for combinatorial Schrödinger operators III-Magnetic fields[END_REF] and to the one of intrinsic metric, e.g.g, [START_REF] Anné | The Gauss-Bonnet operator of an infinite graph[END_REF][START_REF] Huang | A note on selfadjoint extensions of the Laplacian on weighted graphs[END_REF]. In particular, they show that the Hodge Laplacian is essentially self-adjoint, when the graph is χ-complete. This method was more recently adapted in the case of weighted triangulation in [START_REF] Chebbi | The Discrete Laplacian of a 2-Simplicial Complex[END_REF]. In the present study, the graphs are considered as one-dimensional simplicial complexes. If the graph contains cycles of length 3, we can consider it as two dimensional simplicial complex, see [START_REF] Danijela | Spectra of combinatorial Laplace operators on simplicial complexes[END_REF][START_REF] Duval | Critical groups of simplicial complexes[END_REF].

The current study has two major aims. The first aim is to consider the notion of χ-completeness and we approach this geometric hypothesis for the weighted magnetic graph. More precisely, we introduce the notion of χ θ -completeness where θ stands for the notion for the magnetic potential acting on edges. We start with showing that this notion covers many situations, we turn to see Theorem 4.4. In particular, in the case of a simple increasing arithmetic graph, we refer to Section 4.3 for a precise definition, we prove that G is χ-complete. In contrast, if we define the magnetic potential to be θ x,y := 1 d(x)d(y)

, where d(.) is the degree of vertices we refer to Section 2 for the definition, we have G is not χ θ -complete. Then, we study the question of essential self-adjointness of the discrete magnetic Laplacian via the notion of χ θ -completeness. More precisely, we prove that the discrete magnetic Laplacian is essentially self-adjoint when the magnetic graph is χ θ -complete. The technique of the proof is inspired from [START_REF] Anné | The Gauss-Bonnet operator of an infinite graph[END_REF]. At the end of this part, we deal with the question of essentially self-adjointness via Nelson' Lemma and quasi-analytic vectors approach. These two techniques are explained in [START_REF] Baloudi | The adjacency matrix and the discrete Laplacian acting on forms[END_REF] in the case of non magnetic graphs. Our remark is that these two criteria remain applicable for a magnetic graphs and are independent of the magnetic factor on the graph.

The second aim of this paper is to study the magnetic adjacency matrix acting on graph. First, we define the difference magnetic operators d θ acting on 0-forms. We denote by δ θ it formal adjoint. This permits to define two different types of discrete operators associated to a locally finite magnetic graph: ∆ 0,θ = δ θ d θ and ∆ 1,θ = d θ δ θ . In the case of non-magnetic graphs (i.e., θ = 0), the expression of ∆ 1,0 is different from that obtained in [START_REF] Baloudi | The adjacency matrix and the discrete Laplacian acting on forms[END_REF]. Indeed, the discrete operator obtained in [START_REF] Baloudi | The adjacency matrix and the discrete Laplacian acting on forms[END_REF] was built by choosing a skew-symmetric statistic, namely fermionic one. Second, we study the link between the magnetic adjacency matrix A θ acting on line graph and the operator ∆ skew 1,θ = ∆ 1,θ | D(∆ 1,θ )∩C skew (E) where D(∆ 1,θ ) denotes the domain of ∆ 1,θ and C skew (E) denotes the set of 1-cochains corresponds to fermionic statistics. In the case of bi-partite graph, by analysing the structure of ∆ skew 1,θ , we prove the existence of a multiplication operator Q(R) such that A θ + Q(R) and ∆ skew 1,θ are unitarily equivalent.

The paper is structured as follows: The next Section is devoted to some definitions and notations for graph theory. In the third section, we define two discrete magnetic Laplacians operators ∆ 0,θ associated to vertices and ∆ 1,θ associated to edges. In the fourth section, we discuss the notion of χ θ -completeness and we give some examples. In the last section, we establish a link between the magnetic adjacency matrix acting on vertices on line graph and the discrete magnetic Laplacian action on edges. Acknowledgment. The authors thank Echi Nadhem for the useful discussions and comments on the paper.

Generalities about graphs

We start with some definitions and fix our notations for graphs. We refer to [START_REF] Chung | Spectral graph theory[END_REF][START_REF] Mohar | A survey on spectra of infinite graphs[END_REF][START_REF] Jeribi | Spectral theory and applications of linear operators and block operator matrices[END_REF] for surveys on the matter. A weighted graph G is a triple (V, E, m) where V is a set at most countable (vertices), E : V × V -→ [0, +∞) symmetric (edges) and m : V -→ (0, +∞). We say that two vertices x, y ∈ V are neighbors if E(x, y) = 0. In this case, we write x ∼ y. We say that there is a loop in x ∈ V, if x ∼ x. The set of neighbors of x ∈ V is denoted by

N G (x) := {y ∈ V : x ∼ y}.
A graph is locally finite if ♯N G (x) for all x ∈ V. The weighted degree of vertices is given by

d V (x) := 1 m(x) y∈V E(x, y).
A graph G is connected, if for all x, y ∈ V, there exists an xy-path, i.e., there is a finite set (x 0 , ..., x n ) ∈ V n such that x 0 = x, x n = y and E(x i , x i+1 ) > 0 for all i ∈ {0, ..., n -1}. The minimal possible n is denoted by ρ V (x, y) and called the distance between x and y. If no edges appear more than once in (x 1 , ..., x n ), the path is called a simple path. The path is called a cycle or closed when the origin and the end are identical, i.e., x 1 = x n . An n-cycle is a cycle with n vertices. Finally, as we are dealing with magnetic fields, we fix the magnetic potential:

θ : V × V -→ [-π, π], such that θ(x, y) = -θ(x, y).
We denote θ(x, y) by θ x,y . We say that G := (V, E, m, θ) is an unoriented weighted magnetic graph. A magnetic graph G is simple if m = 1, E ∈ {0, 1} and θ = 0. When G is simple, we have d V (x) = ♯N (x) for all x ∈ V.

In the sequel, we assume that all graphs are locally finite, connected and have no loop.

The discrete magnetic Laplacian

In this section, we will always consider a weighted magnetic graph G = (V, E, m, θ). We define the set of 0-cochains (or 0-forms)

C(V) := f : V -→ C .
We denote by C c (V) the 0-cochains with finite support in V. Let us define the Hilbert space ℓ 2 (V) as the sets of 0-cochains with finite norm, i.e.,

ℓ 2 (V) := f ∈ C(V) : f 2 := x∈V m(x)|f (x)| 2 < ∞ .
The associated scalar product is given by

f, g := x∈V m(x)f (x)g(x)
for all f, g ∈ ℓ 2 (V). Concerning the set of edges, we define:

C(E) := f : V × V -→ C . Set C sym (E) := f : V × V -→ C : f (x, y) = f (y, x) and C skew (E) := f : V × V -→ C : f (x, y) = -f (y, x) .
where C skew (E) corresponds to fermionic statistics and C sym (E) corresponds to bosanic statistics. The set of functions in C(E) with finite support is denoted by C c (E). We turn to the Hilbert structure on the set of edges:

ℓ 2 (E) := f ∈ C(E) :
1 2

x,y∈V

E(x, y)|f (x, y)| 2 < ∞ .
The associated scalar product is given by f, g := 1 2

x,y∈V E(x, y)f (x, y)g(x, y) for all f, g ∈ ℓ 2 (E). The difference magnetic operator is the operator

d θ : C c (V) -→ C c (E)
given by

d θ f (x, y) := f (y) -e iθx,y f (x).
The magnetic coboundary operator is the formal adjoint of d θ . It is the operator

δ θ : C c (E) -→ C c (V)
given by

δ θ g(x) := 1 2m(x) y E(x, y)(g(y, x) -e iθy,x g(x, y))
for all g ∈ C c (E) and x ∈ V. Thus it satisfies

d θ f, g = f, δ θ g
for all f ∈ C c (V) and g ∈ C c (E). The operators δ θ and d δ are closable. Indeed,

δ θ : ℓ 2 (E) -→ ℓ 2 (V) (resp. d θ : ℓ 2 (V) -→ ℓ 2 (E)) is with dense domain then δ θ ( resp. d) is closable.
The magnetic discrete Laplacian operator acting on 0-form is given by

∆ 0,θ (f )(x) := δ θ d θ (f )(x) = 1 m(x) y E(x, y)(f (x) -e iθy,x f (y)),
for all f ∈ C c (V) and x ∈ V. The case ∆ 0,0 is so-called physical Laplacian and ∆ 0,π is the t he singless Laplacian, see [START_REF] Cvetković | Towards a spectral theory of graphs based on the signless Laplacian, II[END_REF]. If G is bi-partite, i.e, is a graph whose vertices can be divided into disjoint sets

V 1 and V 2 such that V = V 1 ∪ V 2 and E(x, y) = 0 for all (x, y) ∈ V 2 1 ∪ V 2 2
, then ∆ 0,0 and ∆ 0,π are unitarily equivalent, see [START_REF] Baloudi | The adjacency matrix and the discrete Laplacian acting on forms[END_REF]. If G is simple, the discrete Laplacian ∆ 0,0 is essentially self-adjoint on C c (V), see [START_REF] Wojciechowski | Stochastic compactetness of graph[END_REF]. The magnetic discrete Laplacian operator acting on C c (E) is given by

∆ 1,θ (f )(x, y) := d θ δ θ (f )(x, y) = 1 2m(y) z E(z, y)(f (z, y) -e iθz,y f (y, z)) - e iθxy 2m(x) z E(x, z)(f (z, x) -e iθz,x f (x, z)),
for all f ∈ C c (E) and x ∼ y. Both operators are symmetric and thus closable. We denote the closure by the same symbol ∆ 0,θ (resp. ∆ 1,θ , its domain by D(∆ 0,θ ) (resp. D(∆ 0,θ ), and its adjoint by ∆ * 0,θ (resp.

∆ * 1,θ ). In this paper, the graphs are considered as one-dimensional simplicial complexes. If the graph contains cycles of length 3, we can consider it as two dimensional simplicial complex. In this case, the expression of ∆ 1,0 is different because contains a part acting on the triangular faces, see [START_REF] Baloudi | The discrete laplacian acting on 2-forms and application[END_REF][START_REF] Chebbi | The Discrete Laplacian of a 2-Simplicial Complex[END_REF][START_REF] Lim | Hodge Laplacians on graphs, Geometry and Topology in Statistical Inference[END_REF]. The magnetic Gauβ-operator is defined on C c (V) ⊕ C c (E) by:

D θ := d θ + δ θ ∼ = 0 δ θ d θ 0 .
This operator is of Dirac type and is motived by the Hodge Laplacian:

∆ θ,HL := ∆ θ := (d θ + δ θ ) 2 ≃ ∆ 0,θ ⊕ ∆ 1,θ .
We refer to [START_REF] Anné | The Gauss-Bonnet operator of an infinite graph[END_REF] for the case of non magnetic field. Notation: We denote by N the non-negative integer. In particular, 0 ∈ N. For x ∈ V, we set d(x) := ♯N G (x) and by N * the one of positive integers.

Essential self-adjointness

In this section, we discuss the question of essential self-adjointness via the notion of χ-completeness which was introduced and investigated by Torki in [START_REF] Anné | The Gauss-Bonnet operator of an infinite graph[END_REF]. We approach this geometric notion and we discuss some examples.

Geometric Hypothesis.

We first recall the notion of χ-completeness in the case of non-magnetic graph obtained in [START_REF] Anné | The Gauss-Bonnet operator of an infinite graph[END_REF]:

Definition 4.1. The graph G := (V, E, m) is χ-complete if
there exists an increasing sequence of finite set S n such that V = ∪ n S n and there exist related functions χ n satisfying:

1) χ n ∈ C c (V), 0 ≤ χ n ≤ 1, 2) χ n (x) = 1 for all x ∈ S n . 3) ∃M > 0, ∀n ∈ N, x ∈ V, such that 1 m(x) y E(x, y) | dχ n (x, y) | 2 ≤ M.
In [START_REF] Anné | The Gauss-Bonnet operator of an infinite graph[END_REF], in the case of non-magnetic graph, the authors schow that the two discrete Laplacian operators are essentially self-adjoint, when the graph is χ-complete. In the case of magnetic graph, we do no know that if the χ-completeness hypothesis implis that ∆ 0,θ and ∆ 1,θ are essentially self-adjoint. We introduce the following definition: Definition 4.2. Let G := (V, E, m, θ) be a magnetic weighted graph and let (O n ) n be an increasing sequence of finite sets such that V = ∪ n O n and there exists

(η n ) n and φ n with i) η n ∈ C c (V), 0 ≤ η n ≤ 1 and η n (x) = 1 for all x ∈ O n , ii) φ n ∈ C(V) such that φ n converge to 0, Let χ n = η n e iφn .
We say that G is χ θ -complete if the exist is a non-negative c such that for all x ∈ V and n ∈ N we have

(1) 1 m(x) y E(x, y) | d θ χ n (x, y) | 2 < c. Remark 4.3. Let G := (V, E, m, θ) be a weighted magnetic graph. If G is χ θ -complete then (2) 1 m(x) y E(x, y) | sin θ xy 2 | 2 < c.
Now, we ask the following question: Is there a graph that must χ-complete if it is without magnetic potential and not χ θ -complete if with nonzero magnetic potential? We recall that the graph is assumed connected. For x 0 , x ∈ V and n ∈ N, we set

B n (x 0 ) := {y ∈ V : ρ V (x 0 , y) ≤ n}, |x| := ρ V (x 0 , x) and d Bn(x0) := d V | Bn(x0) .
Theorem 4.4. Let G := (V, E, m, θ) be a magnetic graph with the phase 4 ).

θ x,y =      1 d 2 (x) -d 2 (y) , if d(x) = d(y), 0, otherwise. Let x 0 ∈ V. Suppose that sup x∈V y E(x, y) m(x)d(x)d(y) < ∞ and d Bn(x0) = O(n 1
Then G is χ θ -complete.

Proof. Let n ∈ N * , x ∈ V and set

η n (x) :=                n |x| - n|x| -n 2 n 8 -n 5 + 2n 4 -n + 1 , if n ≤ |x| ≤ n 4 + 1, 1, if x ∈ B n (x 0 ), 0, otherwise and
φ n (x) := |x| n .
We note that n 4 + 1 satisfies the following equation

(3) -nX 2 + n 2 X + n 9 -n 6 + 2n 5 -n 2 + n = 0. Let J(x) = 1 m(x) y E(x, y) | d θ χ n,θ (x, y) | 2 .
The fist case: If x ∈ B n (x 0 ), we take into account the equation ( 3), then for some c, C > 0 independent from x and n, we obtain

J(x) ≤ 4 m(x) y∈Bn(x0) E(x, y) | sin φ n(y) -φ n (x) + θ x,y 2 | 2 + 2 m(x) |y|=n+1 E(x, y) n 8 -n 5 + 2n 4 + n 2 + 1 (n + 1)(n 8 -n 5 + 2n 4 -n + 1) 2 + 8 m(x) |y|=n+1 E(x, y) n 9 -n 6 + 2n 5 -2n 2 (n + 1)(n 8 -n 5 + 2n 4 -n + 1) sin( φ n (y) -φ n (x) + θ y,x 2 ) 2 ≤ c m(x) y∈Bn+1(x0) E(x, y) 1 n 2 + 1 d(x)d(y) ≤ C.
The second case: |x| = n 4 + 1, we take into account the equation (3) then, for some C > 0 independent from x and n, we obtain

J(x) = |y|=n 4 E(x, y) 2n 4 -n + 1 n 11 -n 8 + 2n 7 -n 4 + n 3 ≤ C m(x) y∈B n 4 (x) E(x, y) n 7 .
The thirst case: If x ∈ B n 4 \B n , then

J(x) ≤ c m(x) |y|=|x| E(x, y) | sin( φ n (y) -φ n (x) + θ y,x 2 ) | 2 + c m(x) |y|=|x|-1 E(x, y) | sin( φ n (y) -φ n (x) + θ y,x 2 ) | 2 + c m(x) |y|=|x|-1 E(x, y) | (n 8 -n 5 + 2n 4 -n + 1) + n|x|(|x| -1) |x|(|x| -1)(n 8 -n 5 + 2n 4 -n + 1) | 2 + c m(x) |y|=|x|+1 E(x, y) | sin( φ n (y) -φ n (x) + θ y,x 2 ) | 2 + c m(x) |y|=|x|+1 E(x, y) | (n 8 -n 5 + 2n 4 -n + 1) + n|x|(|x| + 1) |x|(|x| + 1)(n 8 -n 5 + 2n 4 -n + 1) | 2 ≤ C m(x) y∈B n 4 +1 E(x, y) 1 n 2 + 1 d(x)d(y)
.

the fourth case: If |x| > n 4 + 1, we take into account the equation (3), we obtain J(x) = 0.

Remark 4.5. Let G := (V, E, m, θ) be a magnetic graph with the phase

θ x,y =      1 d 2p (x) -d 2p (y) , if d(x) = d(y), 0, otherwise.
By repeating the proof of Theorem 4.4, we obtain that if

sup x∈V y E(x, y) m(x)d p (x)d p (y) < ∞ and d Bn(x0) = O(n 1 4 ), then G is χ θ -complete.
Example 4.6. Let α > 1. Let G := (V, E, m, θ) be a magnetic graph with the phase given in Theorem 4.4 and

E(x, y) := m(x)m(y)d(x)d(y) (m(x) + m(y))(|x| + |y|) α and d(x) ≤ n 1 4 for all x ∈ B n . Then, G is χ θ -complete.
We construct an example of graph that is χ-complete and not χ θ -complete.

Case of the class of anti-tree.

Let x 0 ∈ V and set S n (x 0 ) := {y ∈ V : ρ V (x 0 , y 0 ) ≤ n} the sphere of radius n ∈ N around x 0 . Definition 4.7. Let G := (V, E, m, θ) be a magnetic graph. We say that G is anti-tree if there is x 0 ∈ V such that

N G (x) = S n-1 (x 0 ) ∪ S n+1 (x 0 ) for all n ∈ N * and x ∈ S n (x 0 ). S 0 S 1 S 2
An anti-tree with sphere S 0 , S 1 , S 2 of cardinal 1, 17, 4.

Let G := (V, E, m, θ) be a magnetic anti-tree graph with m = 1, E ∈ {0, 1} and

θ x,y := |x| -|y| 2 π. If d(.) is unbounded, then G is not χ θ -complete. Indeed, y∼x sin 2 ( θ x,y 2 ) = d 2 (x) 2 ,
for all x ∈ V. Using Remark 4.3, we have G is not χ θ -complete. In contrast, In [START_REF] Baloudi | The adjacency matrix and the discrete Laplacian acting on forms[END_REF], the authors prove that if

n∈N 1 ♯S n + ♯S n+1 = +∞,
then G is χ-complete.

Case of increasing arithmetic graph.

We now strengthen the previous example and follow the ideas of [START_REF] Bonnefont | Essential spectrum and Weyl asymptotics for discrete Laplacians[END_REF].

Definition 4.8. A 1-dimensional decomposition of the graph G is a family of finite sets (S n ) n≥0 which forms a partition of V, that is V = ⊔ n≥0 S n , and such that for all x ∈ S n , y ∈ S m ,

x ∼ y =⇒ |n -m| ≤ 1.

We now introduce an example of 1-dimensional decomposition:

Definition 4.9. Let G := (V, E, m) be a weighted locally finite graph and let (S n ) n be a 1-dimensional decomposition of the graph G. We say that G is increasing arithmetic graph if 1) x, y ∈ S n =⇒ E(x, y) = 0 2) There are two integer a, k ∈ N such that d(x) = a + nk for all x ∈ S n .

Let G := (V, E, m, θ) be a magnetic anti-tree graph with m = 1, E ∈ {0, 1} and

θ x,y := 1 d(x) -d(y) . If d(.) is unbounded then G is not χ θ -complete. Indeed, y∼x sin 2 ( θ x,y 2 ) = d(x) sin 2 1 2k .
for all x ∈ V.

Remark 4.10. Let G be a simple increasing arithmetic graph. Then

∞ n=1 1 a + n + a - n = ∞ n=1 1 √ a + nk = ∞, when a + n = sup x∈Sn ♯N G(x)) ∩ S n+1 and a - n = sup x∈Sn ♯N G(x)) ∩ S n-1 . Using [4, Theorem 3.20],
we infer that G is χ-complete.

For f, h ∈ C c (V) and g ∈ C(E), we set:

(f + g)(x, y) = f (y)g(x, y), (h + f -)(x, y) = h(y)f (x), (f -g)(x, y) = f (x)g(x, y). and f (x, y) = f (x) + f (y) 2 .
Remark 4.11. 1) For g ∈ C c (E) and f ∈ C c (V), we have

δ θ ( f g)(x) = f (x)δ θ g(x) 2 + 1 4m(x) y E(x, y)f (y) g(y, x
)e iθy,x g(x, y) .

2) For f, g ∈ C c (V), we have

d θ (f g)(x, y) = f (y)d θ (g)(x, y) + g(x)d θ (f )(x, y) + (e iθy,x -1)f (y)g(x)
.

3)

If d θ = δ * θ and δ θ = d * θ , then D δ is essentially self-adjoint. Theorem 4.12. Let G := (V, E, m, θ) be a χ θ -complete graph, then D θ is essentially self-adjoint.
Proof. Let f ∈ D(δ * ). By using Remark 4.3, we have

(4) (1 -χ n )f 2 ≤ 4 f 2 , (1 -χ n ) + d θ f 2 ≤ 4 d θ f 2 , (5) f + d θ (1 -χ n ) 2 ≤ x∈V m(x)|f (x)| 2 1 m(x) y E(x, y)|d θ (1 -χ n )(x, y)| 2 ≤ 2c f 2 and (6) (e iθ.,. -1)(1 -χ n ) -f + 2 ≤ x,y E(x, y) sin 2 θ y,x 2 |1 -χ n (y)| 2 |f (x)| 2 ≤ c f 2 .
On the other hand, we have

(7) d θ (f -χ n f ) = f (y)d θ (1 -χ n )(x, y) + (1 -χ n )(x)d θ f (x, y) + (e iθx,y -1)f (y)(1 -χ n )(x).
Taking into account Eqs(4-7) and the Dominated convergence theorem, we obtain

lim n (f -χ n f ) + d θ (f -χ n f ) = 0.
This proves that

d θ = δ * θ . Now, let f ∈ D(d * ), then (1 -χ n )δ θ f ≤ 2 δ θ f 2 and 1 4m(.) y E(., y)(1 -χ n )(y)(f (y, .
)e iθ.,y y(., y)) 2 ≤ 1 4

x∈V y

E(x, y) | f (y, x) -e iθy,x f (x, y) | 2 ≤ δ θ f 2 .
By the following equality

δ θ ((1 -χ n )f )(x) = (1 -χ n )(x)δ θ f (x) 2 + 1 4m(x) y E(x, y)(1 -χ n )(y)(f (y, x) -e iθy,x f (x, y)) we deduce that lim n (f -χ n f ) + δ θ (f -χ n f ) = 0.
So, δ θ = d * θ . Using Remark 4.11, we infer that D θ is essentially self-adjoint

Inspired from [START_REF] Anné | The Gauss-Bonnet operator of an infinite graph[END_REF]Proposition13] and [30, Theorem VIII], we obtain the following consequence:

d G (x) -1 2n = ∞.
The technique of the proof is based on the quasi-analytic vectors approach. We can adapt this approach for ∆ 1,θ . Indeed, let

g ∈ C c (E), k ∈ N, then 1 B 2 k (0) g, ∆ 1,θ 1 B 2 k (0) g = x∈V 1 4m(x) | y E(x, y) (1 B 2 k (o) g)(y, x) -e iθy,x (1 B 2 k (o) g)(x, y) | 2 ≤ x∈V 1 4m(x) t E(x, t) y E(x, y) | (1 B k (o) 2 g)(y, x) -e iθy,x (1 B k (o) 2 g)(x, y) | 2 ≤ x∈V 1 2m(x) t E(x, t) y E(x, y) | (1 B 2 k (o) g)(y, x) | 2 + | (1 B 2 k (o) g)(x, y) | 2 ≤ sup x∈B 2 k (o) d G (x) 1 B 2 k (o) g 2 .
This proves that if the hypothesis (8) is satisfied, then ∆ 1,θ is essentially self-adjoint on C c (E). We turn to the Nelson criterion, we recall [4, Theorem 5.13]:

Suppose that sup

x,y∈V,x∼y z∈V

1 m(x) E(x, z)|M(x, y) -M(x, z)| 2 < ∞. Then, ∆ 1,0 | C c skew (E) is essentially self-adjoint on C c skew (E) and ∆ 1,0 | C c sym (E) is essentially self-adjoint on C c sym (E).
This criterion remains applicable for a magnetic graphs. Indeed, the terms f (x, y) in the proof of [4, Theorem 5.13] are replaced by f (y, x)e iθy,x g(x, y). Remark 4.15. the Stieltjes and Nelson criterium are independent of the magnetic factor on the graph. On the contrary, the χ θ depends on the choice of the phase on the edges. 4.5. Some remarks on the spectrum. Let X be a Banach space. We denote by C(X) the set of all closed densely defined linear operators on X. For A ∈ C(X), we write σ(A) (resp. ρ(A)) for the spectrum (resp. the resolvent set ) of A. We denote by K(X) the set of all compact operators on X to itself. We define the essential spectrum of the operator A by

σ ess (A) = K∈K(X) σ(A + K).
In [START_REF] Ayadi | Spectra of Laplacians on forms an infinite graph[END_REF], the author studies the relationship between σ(∆ 0,0 ) and σ(∆ 1,0 | C c (E) ). More precisely, the author proves that the nonzero spectrum of σ(∆ 0,0 ) and σ(∆ 1,0 | C c (E) ) is the same when d V (.) is bounded by using the Weyl's criterion. This result remains true under the influence of magnetic potential on the edges of G. Indeed, ∆ 0,θ and ∆ 1,θ are non-negative operator, then ρ(∆ 0,θ ) and ρ(∆ 1,θ ) are non empty. By using [20, Theorem 1.1], we have

σ(∆ 0,θ )\{0} = σ(∆ 1,θ )\{0}.
Going over the same techniques of the proof of [START_REF] Ayadi | Spectra of Laplacians on forms an infinite graph[END_REF], we obtain

σ ess (∆ 0,θ )\{0} = σ ess (∆ 1,θ )\{0} when d V (.) is bounded.

Application to the magnetic adjacency matrix

The spectral theory of adjacency matrices acting on graphs is useful for the study of some gelling polymers, of some electrical networks, and number theory, e.g., [START_REF] Doyle | Random walks and electric networks[END_REF][START_REF] Davidoff | Elementary number theory, group theory, and Ramanujan graphs[END_REF][START_REF] Mohar | The spectrum of infinite graphs with bounded vertex degrees, Graphs, hypergraphs and applications[END_REF]. In this section, we investigate the relationship between the magnetic adjacency matrix and the magnetic discrete laplacian acting on pseudo-cochains. Let G := (V, E, m, θ) be a weighted magnetic graph. The magnetic adjacency matrix is given by

A θ G (f )(x) := 1 m(x) y E(x, y)f (y)e iθy,x
with f ∈ C c (V) and x ∈ V. It is symmetric and thus closable. We denote it closure by the same symbol. When G is simple, A 0 G is unbounded if and only if it is unbounded from above and if and only if the degree is unbounded see [START_REF] Golénia | Unboundedness of adjacency matrices of locally finite graphs[END_REF]. To achieve this goal, we analyze the structure of ∆ 1,θ . Note that

∆ 1,θ | C c skew (E) := Q(∆ skew 1,θ )
Where Q is the operator multiplication by the function (x, y) -→ e i θx,y

+ e i θx,y 2 m(x) z E(x, z)e i θz,x 2 cos θ zx 2 f (x, z).
Concerning the symmetric choice, we set:

d + θ : C c (V) ⊂ ℓ 2 (V) -→ C c (E) ⊂ ℓ 2 (E) d + θ (f )(x, y) := f (y) + e iθx,y f (x). The formal adjoint of d + θ is given by δ + θ : C c (E) ⊂ ℓ 2 (E) -→ C c (V) ⊂ ℓ 2 (V) δ + θ (f )(x) := 1 2m(x) y E(x, y)(f (y, x) + e iθy,x f (x, y)).
We start with new and direct remarks: In this case, we say that

G is χ + θ -complete. Remark 5.2. If G is χ + θ -complete, then there is c > 0 such that 1 m(x) y E(x, y) cos 2 θ x,y 2 < c
for all x ∈ V. Now, set

D + θ := 0 δ + θ d + θ 0 .
By going through the same technique of the proof of Theorem 4.12, we have D + θ is essentially self-adjoint when G is χ + θ -complete. Remark 5.3. The proof of Theorem 4.4 does not apply line by line in the χ + θ -completeness case. As above, d + θ and δ + θ are both closable. We keep the same symbol for the closure. A direct computation gives that:

(∆ + 1,θ f )(x, y) := d + θ δ + θ (f )(x, y) = 1 2m(y) z E(z, y)(f (z, y) + e iθz,y f (y, z)) + e iθxy 2m(x) z E(x, z)(f (z, x) + e iθz,x f (x, z))
We set:

∆ sym 1,θ (f )(x, y) := e i θy,x 2 ∆ 1,θ | C c sym (E) (f )(x, y) = e i θy,x 2 m(y) z E(z, y)e i θz,y 2 cos θ y,z 2 f (z, y) + e i θx,y 2 m(x) z E(x, z)e i θz,x 2 cos θ z,x 2 f (x, z).
We have that ∆ sym 1,θ and ∆ skew 1,θ have the same expression. However, they do not act on the same space. In the case of a bi-partite graph, we prove that the two operators are unitarily equivalent. Proof. We consider the bi-partite decomposition {V 1 , V 2 }. Let F : ℓ 2 skew (E) -→ ℓ 2 sym (E) be the unitarily map given by

F (f )(x, y) := S(x, y)f (x, y) where S(x, y) :=    1, if (x, y) ∈ V 1 × V 2 , -1, if (x, y) ∈ V 2 × V 1 .
Let H be the following mapping from ℓ 2 sym (E) into ℓ 2 skew (E) given by:

H(f )(x, y) := S(x, y)f (x, y).
Then, F f, g = f, Hg , F H(g) = g and HF (f ) = f for all f ∈ ℓ 2 sym (E) and g ∈ ℓ 2 skew (E). So, we have Corollary 5.8. Let G := (V, E, m, θ) be a weighted magnetic bi-partite graph. If G is χ θ -complete then A θ G is essentially self-adjoint.

H(f ) = F -1 (f ) = F * (f ) For all f ∈ ℓ 2 sym (E). Therefore, F ∆ skew 1,θ F -1 (f )(x,
In [START_REF] Golénia | The magnetic Laplacian acting on discrete cusps[END_REF], they introduce the definition of the flux of a magnetic potential:

Definition 5.9. The space of cycles of G, denoted by Z 1 (G), is the Z-module with a basis of geometric cycles γ := (x 0 , x 1 ) + (x 1 , x 2 ) + .... + (x N -1 , x N ) with i=0,...,N-1, x i ∼ x i+1 , and x N = x 0 . A holonomy map is the map

Hol θ : Z 1 (G) -→ R/2πZ
given by Hol θ (x 0 , x 1 ) + (x 1 , x 2 ) + .... + (x N -1 , x N ) := θ x0,x1 + ... + θ xN ,x0 .

We denote by ∆ G 0, θ the magnetic physical discrete Laplacian acting on line graph.

Proposition 5.10. Let G := (V, E, m, θ) be a magnetic tree, i.e., a magnetic graph with no closed path. Then, ∆ G 0, θ is unitarily equivalent to ∆ G 0,0 .

Proof. Let γ ∈ Z 1 (G). Since G is tree, then there is γ 1 ∈ Z 1 (G) on the form

Corollary 4 . 13 . 4 . 4 .

 41344 Let G := (V, E, m, θ) be a weighted magnetic graph. If G is χ θ -complete then ∆ 0,θ is essentially self-adjoint on C c (V) and ∆ 1,θ is essentially self-adjoint on C c (E). Some remarks on the Stieltjes and Nelson criterium. Let G := (V, E, m) be a weighted graph and o ∈ V. In[START_REF] Baloudi | The adjacency matrix and the discrete Laplacian acting on forms[END_REF], they prove that ∆ 1,0 | C c skew (E) is essentially self-adjoint on C c skew (

Theorem 4 . 14 .

 414 Let G := (E, V, m) be a locally finite graph. Set M(x, y) := 1 + deg G (x) + deg G (y).

2 cos θ zy 2 f

 22 (f )(x, y) := e i θy,x 2 m(y) z E(z, y)e i θz,y (z, y)

Remark 5 . 1 .

 51 The geometric hypothesis given in Section 4.1 can be adapted in the context of d + . Indeed, we replace (ii) in the Definition 4.2 by φ n converge to π and we replace Eq.(2) by 1 m(x) y E(x, y) | χ n (y) + e iθx,y χ n (x) | 2 < c.

Lemma 5 . 4 .

 54 Let G := (V, E, m, θ) be a weighted magnetic bi-partite graph. Then, ∆ sym 1,θ and ∆ skew 1,θ are unitarily equivalent.

2 m 2 cos θ z,y 2 f 2 m 2 cos θ z,x 2 fDefinition 5 . 5 . 1 y 2 fCorollary 5 . 7 .

 222222551257 y) = S(x, y) e i θy,x (y) z E(z, y)S(z, y)e i θz,y (z, y)+ S(x, y) e i θx,y (x) z E(x, z)S(x, z)e i θz,x (x, z)Moreover, noting that S(x, y)S(z, y) = S(x, y)S(x, t) = 1 when E(z, y) > 0 and E(x, t) > 0. We get F ∆ skew 1,θ F -1 = ∆ sym 1,θ . Let G := (V, E, m, θ) be a weighted magnetic graph with the phase θ :E -→]π, π[. Set V = E/ ∼, where (x, y) ∼ (y, x). The magnetic line graph of G is the magnetic graph G := ( V, E, m, θ) where m = 1, E((x, y), (s, t)) := E(x, y)E(s, t) cos θ y=t )1 (x,y) =(s,t)and θ((x, y), (s, t)) := 1 2 (θ xyθ st )(1 x=s + 1 x=t ) + (θ tsθ yx )(1 y=t + 1 y=s ) A graph G A line graph GTheorem 5.6. Let G := (V, E, m, θ) be a weighted magnetic graph. Then, ∆ sym 1,θ is unitarily equivalent toA θ G + Q(R) where R(x, y) := E(x, y) cos θ and Q(R)is the operator of multiplication by R.Proof. Set F : ℓ 2 sym (E) -→ ℓ 2 ( V ) as the function given byF (f )(x, y) := E(x, y) cos θ x,y 2 f (x, y).It is clear thatF -1 (f )(x, y) = F * (f )(x, y) = 1 E(x, y) cos θx,(x, y) for all f ∈ ℓ 2 ( V). Notice now that on C c ( V): F ∆ 1,θ F -1 (f )(x, y) = z E(x,y)E(z, y) cos Let G := (V, E, m, θ) be a weighted magnetic bi-partite graph. If R is bounded, then A θ G ≥sup (x,y)∈E R(x, y).

γ 1 :- 1 i=0(Corollary 5 . 11 .

 11511 = ((x, x 0 ), (x, x 1 )) + ((x, x 1 ), (x,x 2 )) + ... + ((x, x n-1 ), (x, x n )) such that Hol θ (γ) : = Hol θ (γ 1 ) = nθ x,xiθ x,xi+1 ) = 0.Applying [19, Proposition 2.1], the result follows. Let G := (V, E, m, θ) be a magnetic tree. Then, ∆ skew 1,θ is unitarily equivalent to ∆ skew 1,0Proof. Combine Proposition 5.10, Theorem 5.6 and Lemma 5.4.