Antonio Paolillo

Paolo Robuffo Giordano

Matteo Saveriano

Dynamical System-based Imitation Learning for Visual Servoing using the Large Projection Formulation

Nowadays ubiquitous robots must be adaptive and easy to use. To this end, dynamical system-based imitation learning plays an important role. In fact, it allows to realize stable and complex robotic tasks without explicitly coding them, thus facilitating the robot use. However, the adaptation capabilities of dynamical systems have not been fully exploited due to the lack of closed-loop implementations making use of visual feedback. In this regard, the integration of visual information allows higher flexibility to cope with environmental changes. This work presents a dynamical system-based imitation learning for visual servoing, based on the large projection task priority formulation. The proposed scheme enables complex and stable visual tasks, as demonstrated by a simulation analysis and experiments with a robotic manipulator.

I. INTRODUCTION

The level of ubiquity reached by robots imposes strict requirements of adaptability, interactivity, and practicality of use. Today robots have to be easily accessible, also to unskilled users. For this reason, the recent research effort is focused on providing robots with friendly control frameworks. In industries, for instance, the communication between technicians and machines can be eased with gesturebased interactions [START_REF] Paolillo | Towards the integration of a pointing-based human-machine interface in an industrial control system compliant with the IEC 61499 standard[END_REF] or other sorts of intuitive interfaces [START_REF] Lee | Gesture, posture, facial interfaces[END_REF].

In general, to allow non-expert users to conveniently use robots, coding duties have to be lightened or even avoided. This objective has been considered by Imitation Learning (IL) methods that undertake to learn a task from data, rather than explicitly coding it. Such a form of machine learning, otherwise known as programming-by-demonstration [START_REF] Billard | Robot programming by demonstration[END_REF] or learning-from-demonstration [START_REF] Argall | A survey of robot learning from demonstration[END_REF], avoids to explicitly program the robotic task, thus facilitating the large deployment of robots. Among others, Dynamical Systems (DS)-based IL allows to implement the data imitation strategy preserving the stability properties of the original controllers [START_REF] Khansari-Zadeh | Learning stable non-linear dynamical systems with gaussian mixture models[END_REF].

For instance, consider a common insertion task with a robotic manipulator as shown in Fig. 1. The task can be kinesthetically taught to the robot through manual demonstrations. DSs can be used to learn the robot joints motion to realize the task in a stable way without coding it, i.e. imitating the stored demonstrations. An obvious improvement of this kind of methods is the inclusion of exteroception, such as vision, which adds adaptability to the learning framework.

This work was supported by the European Union's Horizon 2020 (1-SWARM, Grant No. 871743), the Hasler foundation (EViRCo), and the Italian Ministry of the University (iNEST, Grant No. ECS 00000043). 1 Dalle Molle Institute for Artificial Intelligence (IDSIA), USI-SUPSI, Lugano, Switzerland antonio.paolillo@idsia.ch 2 CNRS, Univ Rennes, Inria, IRISA, Rennes, France prg@irisa.fr 3 Department of Industrial Engineering (DII), University of Trento, Trento, Italy matteo.saveriano@unitn.it In the considered scenario, in fact, a visual feedback would inform the control about possible changes of the insertion location. In this context, Visual Servoing (VS) [START_REF] Chaumette | Visual servo control. Part I: Basic approaches[END_REF], [START_REF]Visual servo control. Part II: Advanced approaches[END_REF] represents a very suitable strategy to realize adaptive DSs. Indeed, VS is a mature technique and widely used to control robots from visual information. Furthermore, its formalism can be easily interpreted as a DS and, thus, conveniently combined with this class of IL. In a previous work, we have shown that the integration of DS and VS brings advantages to both techniques [START_REF] Paolillo | Learning stable dynamical systems for visual servoing[END_REF]. On the one side, it allows to include visual exteroception into the DS, thus adding adaptability to environmental changes. On the other, it allows to avoid specific programming of additional tasks into the VS. We here propose a novel strategy for DS-based Imitation Learning Visual Servoing (ILVS) built on the large nullspace projection formulation originally proposed in [START_REF] Marey | A new large projection operator for the redundancy framework[END_REF]. After reviewing the related literature (Sec. II) and recalling the basic concepts used in our work (Sec. III), we describe the details of our method in Sec. IV. Results carried out with simulations and experiments are analyzed in Section V; Section VI concludes the papers.

II. RELATED WORK

DS-based IL methods effectively generate stable kinematic motion of robots [START_REF] Ijspeert | Dynamical Movement Primitives: learning attractor models for motor behaviors[END_REF]- [START_REF] Urain | Imitationflow: Learning deep stable stochastic dynamic systems by normalizing flows[END_REF]. The simplest way to retrieve stability is the use of a clock signal to suppress possibly unstable dynamics [START_REF] Ijspeert | Dynamical Movement Primitives: learning attractor models for motor behaviors[END_REF]- [START_REF] Saveriano | An energy-based approach to ensure the stability of learned dynamical systems[END_REF]. Alternatively, stability constraints (e.g., derived from a proper Lyapunov function) can be enforced during the learning process using constrained optimization [START_REF] Khansari-Zadeh | Learning stable non-linear dynamical systems with gaussian mixture models[END_REF], [START_REF] Khansari-Zadeh | Learning control Lyapunov function to ensure stability of dynamical system-based robot reaching motions[END_REF]. In other approaches, the learned inputoutput mapping is forced to be a diffeomorphism, and a stable DS can be obtained using its Jacobian [START_REF] Perrin | Fast diffeomorphic matching to learn globally asymptotically stable nonlinear dynamical systems[END_REF], [START_REF] Urain | Imitationflow: Learning deep stable stochastic dynamic systems by normalizing flows[END_REF].

DSs implementing closed-loop strategies with proprioception have adaptation skills as they are capable of passive interactions [START_REF] Kronander | Passive interaction control with dynamical systems[END_REF]. However, their adaptation potential can be fully reached using exteroception like vision. In this regard, VS laws [START_REF] Chaumette | Visual servo control. Part I: Basic approaches[END_REF], [START_REF]Visual servo control. Part II: Advanced approaches[END_REF] can be interpreted as DSs and used to integrate visual feedback into IL schemes [START_REF] Paolillo | Learning stable dynamical systems for visual servoing[END_REF]. The combination of VS and DS allows to generate complex motion, adapt to environmental changes, and preserve stability.

On the VS side, the basic law performance can be augmented with higher adaptability and flexibility resorting to, e.g., planning [START_REF] Chesi | Global path-planning for constrained and optimal visual servoing[END_REF], [START_REF] Mezouar | Path planning for robust image-based control[END_REF], model predictive control [START_REF] Sauvee | Image based visual servoing through nonlinear model predictive control[END_REF], [START_REF] Allibert | Predictive control for constrained image-based visual servoing[END_REF] or optimization-based control [START_REF] Agravante | Visual servoing in an optimization framework for the whole-body control of humanoid robots[END_REF]- [START_REF] Paolillo | Interlinked visual tracking and robotic manipulation of articulated objects[END_REF]. Planning approaches allow to perform complex VS tasks, but they are normally computed offline and have limited adaptation capabilities. In predictive controllers, instead, one has to balance the anticipatory behavior with the computational cost. This balance can be loosened using a set of pre-computed solutions [START_REF] Paolillo | A memory of motion for visual predictive control tasks[END_REF]. These optimization-based solutions need however to be specifically designed. Exploiting the IL paradigm in a VS scheme would instead allow to easily learn complex visual tasks from few demonstrations. Examples are: Gaussian Mixture Regression (GMR) to learn the pixel to camera motion mapping [START_REF] Pignat | Bayesian Gaussian mixture model for robotic policy imitation[END_REF]; deep neural networks to imitate the surgeon in eye operations [START_REF] Kim | Autonomously navigating a surgical tool inside the eye by learning from demonstration[END_REF]; random forests queried to generate commands for clothes manipulation [START_REF] Jia | Cloth manipulation using random-forest-based imitation learning[END_REF]. Visionbased IL is effective in learning fine operations [START_REF] Paradis | Intermittent visual servoing: Efficiently learning policies robust to instrument changes for highprecision surgical manipulation[END_REF] and can generalize the learned motion to different surfaces [START_REF] Dometios | Vision-based online adaptation of motion primitives to dynamic surfaces: application to an interactive robotic wiping task[END_REF]. In contrast to black-box approaches, our method exploits the controller structure to accurately imitate a trajectory while ensuring stability. Stability guarantees can be obtained if the controller analytical structure is maintained and visual IL is used to provide a reference trajectory [START_REF] Johns | Coarse-to-fine imitation learning: Robot manipulation from a single demonstration[END_REF], [START_REF] Phung | Tool centered learning from demonstration for robotic arms with visual feedback[END_REF]. In [START_REF] Vakanski | An image-based trajectory planning approach for robust robot programming by demonstration[END_REF] the reference of a classical image-based VS is inferred from demonstrations and further refined by an optimizer.

In our previous work [START_REF] Paolillo | Learning stable dynamical systems for visual servoing[END_REF], we have integrated 3 existing DS-based IL methods in the VS scheme. Here, we exploit a specific task priority control paradigm, using the large projection formulation originally introduced in [START_REF] Marey | A new large projection operator for the redundancy framework[END_REF]. The proposed approach conveniently separates the data imitation from the maintenance of closed-loop stability. In practice, the higher priority task is designed so as to guarantee closed-loop stability; the data imitation is handled as secondary task. On the learning side, the scheme greatly simplifies the problem as it allows to use consolidated inference approaches without worrying about stability constraints. On the control side, a modified version of the large projection formulation grants enough degrees-of-freedom to perform complex VS tasks.

III. BACKGROUND A. Dynamical systems for imitation learning

In IL, control affine DS can be used to map a state vector x ∈ R n (e.g., the end-effector position), a control input u ∈ R n , the time t ∈ R + to the state time derivative ẋ ∈ R n :

ẋ = f (x) + u(x, t). (1)
When DSs are used to generate point-to-point motions, one needs to guarantee the stability of (1). The approach in [START_REF] Khansari-Zadeh | Learning control Lyapunov function to ensure stability of dynamical system-based robot reaching motions[END_REF] learns f in (1) to reproduce a set of demonstrations and then compute a stabilizing controller u using a control Lyapunov function. Other methods constrain f to be stable and thus set u = 0 [START_REF] Khansari-Zadeh | Learning stable non-linear dynamical systems with gaussian mixture models[END_REF], [START_REF] Neumann | Learning robot motions with stable dynamical systems under diffeomorphic transformations[END_REF]. The work in [START_REF] Ijspeert | Dynamical Movement Primitives: learning attractor models for motor behaviors[END_REF], [START_REF] Saveriano | Incremental skill learning of stable dynamical systems[END_REF] assume linear and stable f , and use u to reshape the linear motion and follow the demonstrations. Making u vanish with time allows to retrieve the stability. In this work, we use a control law inspired by [START_REF] Paolillo | Towards the integration of a pointing-based human-machine interface in an industrial control system compliant with the IEC 61499 standard[END_REF] to guarantee convergence and follow demonstrated trajectories together.

B. Visual servoing using the large projection operator VS enables robot control using the feedback encoded in sets of visual features [START_REF] Chaumette | Visual servo control. Part I: Basic approaches[END_REF], [START_REF]Visual servo control. Part II: Advanced approaches[END_REF]. In image-based VS, such features are directly extracted from the images, using image processing [START_REF] Marchand | Feature tracking for visual servoing purposes[END_REF] or more sophisticated methods, e.g. neural models [START_REF] Paolillo | Visual servoing with geometrically interpretable neural perception[END_REF]. With eye-in-hand configurations, fixed targets and non-moving scenes, the classic VS writes as v = -λ L+ e. This law generates the velocity v ∈ R 6 that drives the camera to a desired pose, zeroing the error e ∈ R f between the visual features and their corresponding targets. In the VS law, λ is a positive scalar gain tuning the exponential decay of the error. Instead, L+ ∈ R 6×f is the pseudo-inverse of an approximation of the so-called interaction matrix [START_REF] Chaumette | Visual servo control. Part I: Basic approaches[END_REF] that relates the features motion to the camera velocity. The approximation, denoted with the hat over the matrix, is in general due to unknown 3D quantities such as the feature depths. Other tasks can be achieved in parallel to VS using the null-space projection [START_REF] Siciliano | Robotics: Modelling, Planning and Control[END_REF]:

v e = -λ L+ e + P σ (2)
where P = (I 6 -L+ L) is the projector in the null-space of the VS task, being I 6 is the 6 × 6 identity matrix; σ is the velocity command needed to realize a secondary task. The number and the type of features used in the VS scheme are chosen in accordance with the desired task to be executed. In general, the dimension of e needs to be large enough in order to control all the 6 camera degrees of freedom. Therefore, the interaction matrix needs to have full rank and, thus, f ≥ 6. However, full rank of L implies no room for projection in its null-space. A possibility to circumvent this problem is to consider the regulation of the error norm η = ∥e∥ in lieu of the full vector e in the control law, as proposed in [START_REF] Marey | A new large projection operator for the redundancy framework[END_REF], [START_REF] Spica | Coupling active depth estimation and visual servoing via a large projection operator[END_REF]:

v η = -λ η L+ η + P η σ (3)
where both L+ η and P η can be computed in closed-form, i.e.:

L+ η = η L⊤ e e ⊤ L L⊤ e and P η = I 6 - L⊤ ee ⊤ L e ⊤ L L⊤ e . ( 4
)
The advantage of (3) is that the primary task has dimension 1, and the remaining 5 degrees of freedom can be exploited to achieve a secondary task. Indeed, the VS task is actually achieved since, at steady state, e converges to zero for η → 0. However, during the transient, the camera motion is not fully constrained, leaving space to other tasks to be accomplished. The main drawback is that (3) is singular if e is null. Thus, a switching strategy is needed to replace (3) with the classic law (2) in the neighborhood of the singularities [START_REF] Marey | A new large projection operator for the redundancy framework[END_REF]:

v = α(η) v η + 1 -α(η) v e (5)
where the switching function α(η) is defined as follows:

α(η) =    1 if η 1 < η ᾱ(η) -α 0 / (α 1 -α 0 ) if η 0 ≤ η ≤ η 1 0 if η < η 0 (6) with ᾱ(η) = (1 + exp(-ᾱ 0 η-η0 η1-η0 ) + ᾱ1 ) -1 .
The hyperparameters ᾱ0 , ᾱ1 , α 0 , α 1 , η 0 , η 1 are set to design a smooth transition of α in [0, 1]. See [START_REF] Marey | A new large projection operator for the redundancy framework[END_REF] and Sec. V for more details.

IV. PROPOSED APPROACH

We build on the VS law (5)-and thus (3)-to realize our ILVS strategy, achieving the desired visual task while following demonstrated trajectories. We posit that the rationale behind ( 5) is particularly suitable to integrate DS and VS, since it allows to fulfill the visual task at steady-state and imitate the demonstrations during the transient thanks to the large redundancy granted by [START_REF] Khansari-Zadeh | Learning stable non-linear dynamical systems with gaussian mixture models[END_REF]. In practice, the projected quantity P η σ can result in complex camera trajectories as evoked by the DS paradigm. Indeed, the velocity term σ can be inferred from the demonstrations using a proper regression technique and inserted in (3) to realize our DS-based ILVS.

We illustrate how to imitate complex camera trajectories from demonstrations using the VS law with the large projection formulation. In particular, we show the adaptation needed to accurately reproduce the stored trajectories, while preserving the controller stability. Thus, we introduce both the control and learning components of our framework, along with the dataset required to perform the imitation strategy.

A. Error norm tracking with velocity imitation

As explained in Sec. III-B, the space left by the VS task using the scalar error norm η allows to realize a large class of secondary tasks specified by the velocity σ. Thus, we propose to use this secondary task to imitate the demonstrations. However, very complex demonstrations, having a nonexponential (or even non-monotonic) decrease of η, could not be accurately followed as they would conflict with the primary VS task that, instead, forces a "simple" exponential decay of η. To overcome this issue, we propose to track a desired norm trajectory η d (t) that also needs to be learned:

v η = ( ηd -λε) L+ η + P η σ, (7)
thus allowing to implement more complex behaviors for the norm trajectory η d (t). In the proposed law, we define ε = η -η d . The desired norm tracking [START_REF]Visual servo control. Part II: Advanced approaches[END_REF] does not affect the stability of the original controller (3) proposed in [START_REF] Marey | A new large projection operator for the redundancy framework[END_REF]. In fact, by choosing V = ε 2 as Lyapunov candidate, it is

V = 2 ε ε = 2 ε (L η v η -ηd ) = 2ε( ηd -λε)L η L+ η + 2 εL η P η σ -2 ε ηd . (8)
Assuming (as usual) a good approximation of the uncertain 3D parameters in interaction matrix, it follows that L η L+ η = 1, L η P η = 0, and (8) reduces to V = -2λε 2 that is strictly negative for η d ̸ = η. Thus, ε = 0 is an asymptotically stable equilibrium point in its neighborhood: [START_REF]Visual servo control. Part II: Advanced approaches[END_REF] remains locally stable as the norm regulation (3) [START_REF] Marey | A new large projection operator for the redundancy framework[END_REF] and the classic VS [START_REF] Chaumette | Visual servo control. Part I: Basic approaches[END_REF]. Note that σ does not impact the stability and can be learned

3: Execution phase (e, η) ←-compute error(image) σ ←-r σ (e|D) ηd ←-r η (τ |H) η d ←-numerical integration( ηd ) v η ←-compute proposed law(η, η d , ηd , σ) v e ←-compute classic law(e, σ) v ←-compute switching law(v e , v η ) send commands(v)
Repeat until e → 0 with no particular constraint, while the learned η d (t) needs to converge to a proper target, as explained in the next section.

B. Learning VS tasks

We consider data from D demonstrations of VS tasks composed of N samples, consisting in pairs of visual error e d n and the corresponding desired camera commands v d n :

D = e d n , v d n N,D n=1,d=1 . (9)
In computing the visual error, we assume that the target features are known and decided in advance. Its dimension, chosen to meet the requirements detailed in Sec. III-B, is also decided beforehand. Thus, using D, a camera velocity can be simply inferred through the current measurement of the visual error. Recalling that our control paradigm [START_REF]Visual servo control. Part II: Advanced approaches[END_REF] requires also information about the error norm, the following dataset

H = τ d n , ηd n N,D n=1,d=1 (10)
is derived from D by computing the desired error norm η d n = ∥e d n ∥ and its time derivative ηd n for each value of n and d. The dataset H also contains τ d n , a time-like variable that linearly increases from τ d 0 = 0 to τ d N = 1 for each demonstration d. The law [START_REF]Visual servo control. Part II: Advanced approaches[END_REF] greatly simplifies the learning of closed-loop stable VS, as it decouples the velocity imitation from the law convergence. Thus, one can learn separately desired camera velocity and error norm with standard regression techniques.

Using a proper regression function, the velocity σ required to follow the demonstrations can be inferred from the dataset (9) and realized by our strategy (7) as secondary task. Such a regression function is actually a map that, given the dataset D, returns an approximation of the velocity imitating the demonstrations for a certain value of the visual error:

σ = r σ (e | D) . (11)
The regression in [START_REF] Saveriano | Incremental skill learning of stable dynamical systems[END_REF] allows to retrieve the desired camera velocity σ in closed-loop, i.e., starting from the measured Fig. 2. Simulation analysis: demonstrated motion samples contained in the dataset (gray lines) and an exemplar execution of our method (black line); each sample is composed of camera (left) and visual features motion (right). For visualization purposes, only the top view of the camera motion is shown.

error e. This aspect is beneficial as the learning algorithm can adjust the desired velocity to cope with an unexpected variation of the environment (e.g., a change of the target location) or disturb occurring on the camera motion. Furthermore, the proposed control scheme [START_REF]Visual servo control. Part II: Advanced approaches[END_REF] converges to the desired features even if σ has a stable equilibrium at e ̸ = 0 or it diverges. The only requirement is that r σ is smooth. Thus, any regression techniques can be used to learn and retrieve it. Notably, we train directly on D without additional stability constraints that significantly increase the training time and the possibility to fall in poor local minima. In this work, we choose GMR to implement r σ because it produces smooth outputs, is fast enough for online control loops, and the only hyper-parameter (the number of Gaussians) is easy to tune. Similarly, the desired value of the time derivative of the error norm ηd is inferred from data as

ηd = r η (τ | H) , (12)
and η d is obtained by numerical integration. In principle, ηd could also be algebraically retrieved from the learned σ, but this choice cannot ensure convergence as the output of ( 11) may diverge or stay in local minima. Note that r η , instead, takes as input the monotonically increasing time-like signal τ ∈ [0, 1], enabling reliable learning contrary to considering η d as input. Also, being η d a scalar, it is easy to obtain multimodal behaviors, i.e., different values of ηd for the same η d , which are difficult to learn. The use of τ solves this issue as there is a unique pair (η d , ηd ) for each τ . Also for r η , any proper regressor may be used and we consider GMR.

C. Imitation Learning-based Visual Servoing framework

The whole framework implementing our ILVS strategy is summarized in Algorithm 1. The procedure consists of three main phases; the first two are carried out offline, whereas the latter is executed online. The initial phase regards the data collection and processing required to build the datasets D and H used in the learning steps. Thus, as explained in Sec. IV-B, two regression functions (GMRs in our work) are trained to provide the framework with the required demonstrated information. The execution phase is the actual online control loop that realizes our ILVS. Initially, the visual error (and its norm) is computed from the acquired camera image, and used to infer the velocity σ imitating the demonstrations. Similarly, also ηd is inferred from data, and η d computed through numerical integration. These values are used to compute the velocity command as in the law [START_REF]Visual servo control. Part II: Advanced approaches[END_REF]. Using the visual error e, also the classic VS law (2) is computed. Thus, the two laws are combined together as in (5) using the switching strategy [START_REF] Chaumette | Visual servo control. Part I: Basic approaches[END_REF]. Finally, the commands are executed and the loop runs again until convergence.

V. RESULTS

A. Simulations

Our simulations are performed on the LASA dataset [START_REF] Khansari-Zadeh | Learning stable non-linear dynamical systems with gaussian mixture models[END_REF]. This dataset contains 26 uni-modal handwriting motion samples, used as camera lateral trajectories (the vertical displacement is implemented as a constant velocity motion) and augmented with corresponding visual features trajectories as in [START_REF] Paolillo | Learning stable dynamical systems for visual servoing[END_REF]. The visual features are 4 points, corresponding to the image of the 4 vertexes of a squared visual pattern placed at a known position in front of the camera. Each motion sample contains 7 demonstrations; each demonstration has 1000 points. We heuristically set the switching strategy with ᾱ0 = 12, ᾱ1 = 6, η 0 = 10, η 1 = 25, α 0 = 0 and α 1 = 1; the number of Gaussians of the GMRs to 7; the VS gain λ to 3. Each execution is given a length of 1050 samples. In our study, we consider only the linear motion of the camera, which means that both control and learning steps of our framework deal with only position variables.

The qualitative performance of our method in imitating the dataset trajectories is shown in Fig. 2 for 25 of the samples. For each one, the planar camera motion (on the left) and the visual features motion on the image (right) are shown. In each plot, the black line refers to the trajectory executed by our method; the gray traces are the demonstrations; the black dot indicates the starting point (taken as the average of the demonstration starts); the cross is the final point. The plots show that our strategy succeeds in imitating the demonstrations while converging to the target in most samples, with the same learning and control parameters. For two of them (the last plots in Fig. 2), the performance deteriorates due to the presence of loops in all the demonstrated features motion.

For the sake of clarity of the presentation, plots concerning the remaining 26th sample of the dataset are reported in Fig. 3456. More into details, Fig. 3 shows that the proposed method allows to achieve the VS task, bringing the camera to convergence, while imitating the dataset trajectories. As explained in Sec. IV, this is made possible thanks to information extracted from the dataset, as shown in Fig. 4 and5, where the gray, black dashed and black solid lines denote the demonstrations, the GMR inference and the execution, respectively. Note that we run the GMR for 1000 samples, which essentially set the maximum time given to the imitation strategy. After about 900 samples, the system switches to the classic VS law (see the plot of the switch variable α in Fig. 6) and there is no room for null-space tasks. The trajectory of the norm is properly tracked, being it the primary task. Instead, as expected, the inferred velocity is followed with less accuracy since it is considered into the secondary task. At around the 900th sample, the tracking performance degrades as the framework starts switching to the classic VS law. From there on, the behavior converges to the target using the classical VS exponential trend.

Crucially, the tracking of the error norm and the consistent execution of demonstrated velocities as secondary task allow the achievement of the full desired task. In Fig. 3 and4, the green and red lines are the executions without the secondary task, i.e., [START_REF]Visual servo control. Part II: Advanced approaches[END_REF] with σ = 0, and without the norm tracking, i.e., simply (3), respectively. As shown, the combination of these components is essential for the imitation of complex trajectories, i.e. having a non-exponential decay of the η. Fig. 7 shows the capabilities of the framework to cope with external disturbances. During the execution of the task, between samples 300 and 400, a smooth perturbation of maximum amplitude of 0.1 m/s is applied to the three components of the velocity, as visible in the plot on the right in Fig. 7. Nevertheless, our method well handles the unexpected event converging to the target while imitating the demonstrations, properly recovering from the disturbance.

B. Experiments

We consider the task of inserting a screwdriver in an hole with the 7-axis manipulator by Franka Emika [START_REF] Haddadin | The Franka Emika robot: A reference platform for robotics research and education[END_REF], see Fig. 1. The robot is equipped with a camera streaming images of 640 × 480 pixels. A pattern of 4 points is attached to the insertion location, so that the robot can use the image of these points to drive the inserting motion. A classic VS would bring the screwdriver, hold by the robot gripper, to collide with the environment. Classic DSs using proprioception, instead, would fail if the insertion location changed from the demonstrated one. Our experiments show the benefit of integrating VS and DSs using the proposed strategy.

We collect 3 demonstrations by kinestetically showing the robot right insertion moves. Each demonstration has 200 samples; the online execution is given 300 iterations. A greedy search to obtain enough accuracy suggests to use 11 GMR components for ηd , 5 for σ. Also for the experiments, we only consider the camera position. The λ gain is set to 0.5; the switching rule has η 0 = 25, η 1 = 50 and the other parameters as in Sec. V-A. The computed camera velocity is converted to joint velocities commands using the end-effector Jacobian and the camera extrinsic calibration (to transform the velocity from the camera to the end-effector frame).

Fig. 8 shows the performance of our approach in reaching the convergence while imitating the demonstrations. The accompanying video shows the robot performing this task along with the online camera stream; the comparison with the simple norm regulation; and the execution of our framework with the change of the insertion location.

VI. CONCLUSIONS

We presented a dynamical systems-based imitation learning for visual servoing. Our method is built on the large projection formulation, where the scalar norm of the visual error is considered as feedback. Thus, the null-space of the visual task is large enough to project other tasks. We use this paradigm to establish priorities between the realization of stable visual servoing and the imitation of manual demonstrations. Such a strategy enables complex visual tasks with no explicit implementation of complicated trajectory generators.

The effectiveness of our approach has been shown through a simulation analysis, presenting the performance on a dataset of different motion samples. Experiments with a robotic manipulator have challenged the framework in a realworld scenario. These results have shown the benefits of integrating dynamical systems with visual servoing. On the one side, it allows to easily realize complex visual trajectories without explicit coding. On the other, the visual feedback enables adaptation capabilities to environmental changes.

Future work will extend and validate the generality of our approach. Indeed, the large projection formulation might be also applied to classic kinematic dynamical systems. Furthermore, the learning process could be improved by collecting better demonstrations, e.g. using a planner. More complex experiments, including the orientation control, larger motion, and multiple insertions, will be designed. Finally, we aim at realizing adaptive manipulation tasks, to pave the road to the deployment of our strategy in industrial environments.

Fig. 1 .

 1 Fig. 1. In complex robotic tasks, such as tool insertion in narrow locations, dynamical systems allow to easily program the desired motion, and visual servoing allows to adapt the learned motion to changes of the environment.

Algorithm 1 :

 1 DS-based ILVS 1: Data processing phase D, H ←collect and process demos 2: Learning phase r σ ←train regressor(D) r η ←train regressor(H)

Fig. 3 .

 3 Fig. 3. Simulation analysis of one motion sample: lateral camera (left) and visual features (right) trajectories.

Fig. 4 .

 4 Fig. 4. Simulation analysis of one motion sample: norm of the error (top) and its time derivative (bottom).

Fig. 5 .Fig. 6 .

 56 Fig. 5. Simulation analysis of one motion sample: camera velocity.

Fig. 7 .

 7 Fig.7. Simulation analysis on one motion sample with disturb: lateral camera (left) and visual features (center) trajectories, and camera velocity (right).

Fig. 8 .

 8 Fig.8. Experiment with the robotic manipulator: norm of the error (top), camera velocity components (middle) and switching variable (bottom).


domain: phys chim

text:
Giovanni Pireddu
email: benjamin.rotenberg@sorbonne-universite.fr

Frequency-Dependent Impedance of Nanocapacitors from Electrode

published or not. The documents may come

Frequency-dependent impedance of nanocapacitors from electrode charge fluctuations as a probe of electrolyte dynamics

Giovanni Pireddu 1 and Benjamin Rotenberg The frequency-dependent impedance is a fundamental property of electrical components. We show that it can be determined from the equilibrium dynamical fluctuations of the electrode charge in constant-potential molecular simulations, extending in particular a fluctuation-dissipation for the capacitance recovered in the low-frequency limit and provide an illustration on water/gold nanocapacitors. This work opens the way to the interpretation of electrochemical impedance measurements in terms of microscopic mechanisms, directly from the dynamics of the electrolyte, or indirectly via equivalent circuit models as in experiments.

The electrical impedance, which quantifies the current induced by a small oscillatory applied voltage, is a fundamental property of the components of an electric circuit. Electrochemical impedance spectroscopy (EIS) techniques [1,2] are routinely used to characterize batteries [3] and capacitors [4] in the mHz-kHz range, while dielectric relaxation spectroscopy (DRS) techniques allow probing interactions and microscopic dynamics at higher frequencies (MHz-THz) [5,6]. Electrodes in nanofluidic devices also open new avenues for energy conversion from salinity gradients [7,8] or from light [9], for catalysis [10] or sensing charged species and measuring ultra-low flow rates with electrochemical correlation spectroscopy [11,12]. Measurements are usually interpreted in terms of equivalent circuits (EC) [13,14] combining elements with impedance Z(ω) = 1/iCω, R, or iLω, with ω the frequency, for capacitors, resistors or inductors, respectively, or more complex ones such as constant phase elements. The frequency-dependent response highlights the relevant time scales in the dynamics of the polarization of the liquid and the transport of charge carriers within such electrochemical devices. The diffusion and migration of ions can be modeled analytically at the mean-field level [15][16][17] and their coupling with hydrodynamic flows with mesoscopic simulations methods [18], to predict the charging dynamics of a capacitor. The effect of steric and other correlations on the dynamics can be introduced in such descriptions [19] or in classical Density Functional Theory [20,21]. Charge transport within porous electrodes, which are used in many devices, can be modeled using such liquid state theories or simpler constitutive equations, in order to make the link with EC such as the transmission line model (TLM) [22][23][24].

Molecular simulations have greatly contributed to recent progress in the description of electrode/electrolyte interfaces, taking into account the polarization of the electrodes by the ions and solvent molecules of the electrolyte [25]. Several methods have been developed to impose the potential of the metal both in ab initio and classical molecular dynamics (MD) simulations [26][27][28][29][30], which enabled the molecular-scale study of simple and complex electrochemical systems, including supercapacitors [31][32][33][34]. In such constant-potential simulations, the charge of the electrodes fluctuates in response to the dynamical evolution of the electrolyte. The local and global responses of the charge of porous electrodes to an applied voltage step can be monitored to build EC models such as the TLM in order to extrapolate to larger system sizes and longer time scales and bridge the gap with the electrochemical impedance measurements [35,36], or to design optimal charging protocols [37,38]. The capacitance and so-called equivalent series resistance of the system can also be obtained using a recently proposed nonequilibrium "computational amperometry" approach [39].

An alternative strategy to compute the capacitance from MD simulations and to extract information on the interfacial electrolyte is to study the equilibrium fluctuations of the electrode charge [40,41]. The differential capacitance, i.e. the response of the electrode charge Q to a change in voltage ∆Ψ between the electrodes, is related to the variance of the charge distribution via a fluctuation-dissipation relation C diff = ∂⟨Q⟩ ∂∆Ψ = β⟨δQ 2 ⟩, where brakets denote an ensemble average, δQ = Q-⟨Q⟩ and β = 1/k B T with k B the Boltzmann constant and T the temperature. The link between the thermal fluctuations of charge carriers and electrical response was already observed almost a century ago by Nyquist and Johnson, who considered the thermal noise in electronic conductors [42,43], and exploited in the analysis of electrochemical systems [44]. Molecular simulations further allowed to correlate the voltage-dependence of the differential capacitance with changes in the interfacial structure of an ionic liquid or a water-in-salt-electrolyte between graphite electrodes [45][46][47]. More recently the charge fluctuations within electrodes were also linked to the electrode-electrolyte interfacial free energy [START_REF] Scalfi | Proceedings of the National Academy of Sciences[END_REF].

Here, we show that the frequency-dependent impedance can be extracted from the equilibrium dynamical fluctuations of the electrode charge in constant-potential simulations, the above fluctuationdissipation relation for the capacitance being recovered in the low-frequency limit. This opens the possibility to investigate with MD simulations the link between the impedance and the dynamics of the liquid between the electrodes. We then illustrate this on water/gold nanocapacitors. Both experiments and simulations recently highlighted some peculiarities of nanoconfined water, such as the decrease of dielectric constant [START_REF] Fumagalli | [END_REF]50] and its anisotropy close to interfaces [51], as well as its unusual dynamical behavior [52,53] and confinement-induced spectral changes [54,55]. In addition, the dynamics of water on metals proceeds by local reorientation jumps or by rare collective fluctuations, depending on the metal [56][57][58]. Finally, the coupling between the interfacial fluctuations of water and the electronic response of carbon plays an important role in the water-graphite friction [59]. The present approach offers a way to probe the dynamics of the liquid encoded in the electrical impedance measured with electrodes. We consider systems such as illustrated in Fig. 1a, with two rigid electrodes maintained at a fixed potential difference ∆Ψ and separated by a liquid with fixed number of particles and temperature T . The charge of each electrode atom fluctuates in order to satisfy the constantpotential constraint, so that the total charge ±Q of the oppositely charged electrodes (due to global electroneutrality), which is the thermodynamic variable conjugate to ∆Ψ, also fluctuates in response to the thermal fluctuations of the liquid. The complex admittance Y (ω) = 1/Z(ω), with Z(ω) the impedance, quantifies the linear response of the induced electric current, I(t) = Q(t), to a small oscillatory voltage, V (t) = V 0 e iωt , in addition to the voltage ∆Ψ defining the thermodynamic ensemble, from the change in current ⟨∆I(t)⟩ = Y (ω)V (t). Using linear response theory, we show in section I of the Supplemental Material (SM) [60] that the admittance can be expressed using the Laplace transform of the current or charge auto-correlation functions (ACF), as:

Y (ω) = β ∞ 0 ⟨δI(0)δI(t)⟩e -iωt dt = β iω⟨δQ 2 ⟩ + ω 2 ∞ 0 ⟨δQ(0)δQ(t)⟩e -iωt dt . (1)
This fundamental result explicitly links the dynamics of the electrode charge, which reflect the thermal fluctuations of the liquid, with the impedance of the system. It generalizes the above fluctuation-dissipation relation for the capacitance, which is recovered in the low frequency limit ω → 0 since the leading term of Eq. 1 is iωβ⟨δQ 2 ⟩, the admittance of a capacitor with capacitance β⟨δQ 2 ⟩. Importantly, Eq. 1 also offers a practical tool to compute the impedance spectrum from the dynamics of the charge fluctuations over the frequency range corresponding to the time scales accessible to simulations (typically between 10 MHz and 1000 THz for all-atom simulations, or lower if coarse-grained models are used), as shown below.

The fluctuating charges introduced in Ref. [26] as additional degrees of freedom to simulate constant-potential conditions followed an equation of motion with a fictitious mass. Ref. [27] then suggested to determine the set of charges that self-consistently satisfy the constantpotential constraints at each step of the simulation, so that the electrode atom charges are enslaved to the microscopic configurations of the liquid. Such a Born-Oppenheimer (BO) sampling allows to use larger time steps, but also suppresses some charge fluctuations that may contribute to some observables. For example, the suppressed fluctuations lead to a missing term in the capacitance corresponding to that of an empty capacitor, C empty diff [41]. Separating the instantaneous charge as Q(t) = Q BO (t)+Q nBO (t) into the BO and non-BO terms and expanding ⟨δQ 2 ⟩ and ⟨δQ(0)δQ(t)⟩ in Eq. 1 yields

Y (ω) = Y BO (ω) + Y nBO (ω) , (2)
where the first term is obtained by replacing Q by Q BO in Eq. 1 and the second is the remainder. The latter depends on the choice of dynamics for the additional degrees of freedom, which should only manifest themselves at high frequencies since they mimick the electronic response, except for the "static" contribution for ω → 0 where Y nBO (ω) ≈ iωC empty diff . In practice, the non-BO contribution to the capacitance is small compared to the BO one arising from the thermal fluctuations of the liq-uid [41] (see also below). In the following, we drop the BO superscript and refer to the BO charges as simply Q.

As an illustration of the possibilities offered by the analysis of the equilibrium dynamical fluctuations of the electrode charge in constant-potential simulations, we investigate gold/water nanocapacitors (see Fig. 1a). Simulations are performed under periodic boundary conditions in the x and y directions only, with fixed number of molecules and corresponding inter-electrode distance. The latter is determined by prior equilibration at a constant pressure of 1 atm in the direction perpendicular to the electrode surfaces, resulting in d =2. 51, 4.94, 9.76, and 19.42 nm between the first atomic planes of the electrodes. Details on the systems and molecular models can be found in the SM Sec. II [60]. All simulations are performed using the MD code Metalwalls [61], with a potential difference ∆Ψ = 0 V between the electrodes except for the response to a voltage step. Unless otherwise stated, the uncertainty is expressed as standard error computed from 10 blocks of the trajectory.

Fig. 1b shows the BO and non-BO contributions to the differential capacitance as a function of the interelectrode distance d. As indicated above, the former contribution, C BO diff = β⟨δQ 2 ⟩, is much larger than the latter, C empty diff , computed from the properties of the electrodes as described in Ref. 41. The figure also shows that the prediction from continuum electrostatics, ε w ε 0 A el /d, with ε 0 the vacuum permittivity, ε w the relative permittivity of bulk water, and A el the lateral area of the electrodes, only qualitatively reproduces the dependence of the capacitance with d in the considered range. This is due to the failure to capture the potential drop across the interfacial water layers, as previously reported on other surfaces [62]. The prediction can be greatly improved by considering the Dielectric Dividing Surface (DDS), obtained from the interfacial permittivity profiles [63]. As the Gibbs Dividing Surface for density profiles, the DDS locates an equivalent sharp interface between a (water-free) region with permittivity ε 0 and another with ε 0 ε w (see SM Sec. III [60]). In this picture, the capacitance per unit area corresponds to that of three capacitors in series, A el /C DDS = 2w DDS /ε 0 +d DDS /ε 0 ε w , with w DDS ≈ 1.27 Å the distance between the last atomic plane of each electrode and the corresponding DDS, and d DDS = d-2w DDS the width of the equivalent bulk water region, and provides a rather accurate description of the actual capacitance (see Ref. 64 for a related approach for the solvation free energy).

Turning to the dynamics of the charge fluctuations, Fig. 2a shows that the charge ACF displays oscillations at short times and an approximately exponential decay at longer times. Increasing the inter-electrode distance damps the former and slows down the latter. In the linear response regime, ⟨δQ(0)δQ(t)⟩ should reflect the evolution of Q during the (dis)charge of the capacitor. Fig. 2b shows the charge as a function of time before after a change from ∆Ψ = 0 V to 1 V. The plateaus reached after the transient regime are consistent with the capacitance for each d, and the characteristic time to reach the new equilibrium increases with d, consistently with the slower decay of the ACF in Fig. 2a. However, the link between the equilibrium ACF and the out-of-equilibrium response is much deeper: Fig. 2c shows the equivalence between the two when properly rescaled, as expected from the fluctuation-dissipation theorem, indicating that in these systems the response is linear at least up to 1 V (as observed for the capacitance of water/graphite capacitors [62]).

The impedance is usually analyzed in a Nyquist plot, which reports its imaginary and real parts in the complex plane, parameterized by the frequency. The MD results obtained from the ACF using Eq. 1 (see SM Sec. IV [60] for details on the numerical procedure) are shown in Fig. 3, for frequencies ranging approximately from 20 GHz to 300 THz. This range is well above typical EIS measurements [2], which are typically performed for salt solutions. However, no particular dielectric features are expected for pure water in this range [65], while the one we consider overlaps with the ones expected for SPC/E water, the so-called Debye relaxation and librational modes [66]. For a sufficiently large system in the lateral directions, the charge ACF, hence the admittance scales linearly with A el , as shown numerically in SM Sec. II [60]. We therefore plot the impedance multiplied by A el in Fig. 3. As in experiments, we model the data using an EC fitted to the admittance/impedance from MD (see SM Sec. V [60]). The Nyquist plots are typical of capacitors with a vertical asymptote at low frequency. The growing real part with inter-electrode distance reflects an increase in dissipative processes. The corresponding resistance scales approximately as d/A el , as expected for a macroscopic resistor. At high frequency, we observe a sign change in the imaginary part, which reflects the short-time oscillations of the ACF and can be modeled using an inductor in the EC, with L ∝ d/A el .

The frequency-dependence is shown in the Bode plot, Fig. 4, where the norm of the complex admittance is multiplied by the effective water slab width d DDS and divided by the lateral area A el to highlight the scaling relations with system size. The most prominent feature is a common peak, with magnitude inversely proportional to d DDS , located at ω ≈ 1.37 ± 0.03 × 10 14 rad s -1 and corresponding to the short-time oscillations observed in the charge ACF shown in Fig. 2a. Interestingly, this frequency coincides with the one of the previously reported oscillations of the dipole moment of nanoconfined water [52,53] as well as a feature in the bulk permittivity spectrum of water ascribed to librational modes [66]. This points to the link between the water polarization and the electrode charge (see below). The linear scaling of this peak with d DDS suggests that it is mainly due to "bulk" water rather than to an interfacial mechanism.

The broad shoulder at low frequency is well described by a Debye model with a single relaxation time. Indeed, from Eq. 1 the admittance can be approximated at low frequency by

Y (ω) ≈ iωβ⟨δQ 2 ⟩ × (1 -iωτ ) ≈ iωC diff /(1 + iωτ ),
with τ the integral of the normalized charge ACF. The deviations from the scaling with d reflect the interfacial component of the capacitance captured by the DDS model (see the above discussion of Fig. 1). In the EC picture, this low-frequency behavior corresponds to a RC circuit in series, with relaxation time τ = RC. Both estimates of τ coincide and increase with inter-electrode distance approximately as τ ∞ × d/[d + 2(ε w -1)w DDS ] (see inset of Fig. 4 and SM Sec. V [60]). The fitted value τ ∞ ≈ 8.6 ps is close to the characteristic time corresponding to the Debye relaxation in bulk water (τ D ≈ 9 ps) that can be interpreted as aris-ing from the migration of orientational defects in the Hbond network [67], even though the separation between individual and collective motion is challenging [66].

As a first step towards the molecular interpretation of the admittance/impedance, we monitor the components of the dipole moment of the confined water slab in the directions perpendicular and parallel to the electrode surfaces. We show in SM Sec. VI [60] that the former almost perfectly follows the dynamics of the electrode charge, while the latter decays more slowly and is only slightly affected by d. We also show that, under the considered conditions, the perpendicular component of the total water dipole satisfies

M ⊥ wat (t) = Q(t)d Q (t)
, with d Q the instantaneous distance between the charge-weighted atomic positions in each electrode. In the present case of a perfect metal, these positions almost coincide with the first atomic planes for all microscopic configurations (see also Refs. 68, 69 for the charge induced on the same model of gold electrodes), so that M ⊥ wat (t) ≈ Q(t)d, which explains the direct link between the ACF of the electrode charge and that of M ⊥ wat . The situation might be different in the case of imperfect screening of charge and potential within metals, which can be introduced in MD using the Thomas-Fermi model [70]. We use the relation between charge and polarization fluctuations to separate the contributions of interfacial and bulk-like water to admittance. Our analysis reveals that the influence of adlayer water to the librational peak is negligible (see SM Sec. VI [60]).

Overall, the present work offers the possibility to extract the frequency-dependent impedance/admittance of electrochemical systems in molecular simulations from the equilibrium dynamical fluctuations of the electrode charge. This opens the way to the interpretation of impedance measurements in terms of microscopic mechanisms, either directly from the dynamics of the electrolyte, or indirectly by EC models used to analyze the experimental data, provided that all the relevant processes emerge within the time scales probed with atomistic simulations. Beyond the case of pure water considered as an illustration, the next step is to investigate electrolyte solutions. Comparison with previous data (see SM Sec. VII [60]) suggests that ions increase the capacitance and introduce longer timescales in the charge relaxation that can be probed using implicit solvent descriptions [71], allowing to extend simulations to larger time scales and bridge the gap with the lower frequency range probed in DRS and EIS experiments. The theories for the ion dynamics used to model the latter, based e.g. on a frequency-dependent friction [72][73][74] could also benefit from simulations with an explicit solvent and/or electrodes. The frequency range considered here for pure water overlaps with that of THz spectroscopy, which can be used to extract microscopic contributions [75] and to investigate electrode/electrolyte interfaces [76]. Future directions may also include the non-linear dynamical re-sponse to large voltages, using e.g. the statistical tools developed to sample electric current fluctuations in MD simulations [77,78], as well as systems with redox reactions at the surface of the electrodes [79,80].

In this section, we report a more detailed derivation of Eq. 1 in the main text which relates the fluctuations of electrode charges and the admittance. The following derivation is inspired by the one reported in Ref. 1 (see in particular section 7.6), and considers a perturbation of a reference Hamiltonian H 0 written in the form

H tot = H 0 -AF (t), (1)
where F (t) represents the time-dependent control variable perturbing the system and A is its conjugate variable. The response to a generic time-dependent perturbation can be decomposed into its Fourier components, so that it is sufficient to consider the effect of a monochromatic perturbation F (t) = F 0 e iωt . As described in Ref. 1, one should in principle introduce a regularization in the form of a factor e ϵt with ϵ > 0 ensuring the absence of perturbation for t → -∞ and take the limit ϵ → 0 at the end of the calculation, but we drop it in the following for simplicity. Linear response theory then allows expressing the average response of a variable B to the (sufficiently small) perturbation as

⟨∆B(t)⟩ = t -∞ ϕ BA (t -s)F (s)ds , (2)
where ϕ AB (t) is the after-effect function, which can be express as the equilibrium correlation function:

ϕ BA (t) = β⟨δB(t)δ Ȧ⟩, (3)
with δB = B -⟨B⟩ (and similarly for A) and Ȧ denotes the time-derivative or A.

The frequency-dependent electrical impedance Z(ω) is the inverse of the admittance Y (ω), which quantifies the response of the electric current I = Q to a change V in the voltage ∆Ψ between the two oppositely charged electrodes. The conjugate variable of voltage is the total charge ±Q of the latter. We therefore consider the above general results in the particular case A = Q, B = I, and F = V (with the appropriate sign convention to define the voltage). Introducing the monochromatic form of V (t) = V 0 e iωt in Eq. 2, we obtain:

⟨∆I(t)⟩ = t -∞ ϕ IQ (t -s)V (s)ds = V 0 t -∞ ϕ IQ (t -s) exp[iωs]ds =V 0 exp[iωt] t -∞ ϕ IQ (t -s) exp[-iω(t -s)]ds =V 0 exp[iωt] ∞ 0 ϕ IQ (t) exp[-iωt]dt = V (t) ∞ 0 ϕ IQ (t) exp[-iωt]dt .
One can therefore identify the Fourier-Laplace transform of the after-effect function ϕ IQ (t)

with the electrical admittance defined by ⟨∆I(t)⟩ = Y (ω)V (t). From Eq. 3, we obtain the admittance from the autocorrelation function of the current as:

Y (ω) = β ∞ 0 ⟨δI(0)δI(t)⟩e -iωt dt . (4)
An alternative expression as function of the electrode charges can be obtained by noting that ⟨δI(0)δI(t)⟩ = -d 2 ⟨δQ(0)δQ(t)⟩/dt 2 , which follows from Faraday's law. Considering

Laplace transform properties for derivatives, we obtain

Y (ω) = β iω⟨δQ 2 ⟩ + ω 2 ∞ 0 ⟨δQ(0)δQ(t)⟩e -iωt dt , (5)
where we omitted a first derivative contribution d⟨δQ(0)δQ(0 -)⟩/dt. This contribution vanishes since equilibrium autocorrelation functions are even with respect to time, and if we assume that ⟨δQ(0)δQ(t)⟩ is differentiable at t = 0. Eq. 5 allows to determine the frequency-dependent electrical impedance of the system from the dynamics of the equilibrium fluctuations of the electrode charge, which can be sampled in molecular dynamics simulations. Each simulation box is prepared with a preliminary equilibration in the N P z T ensemble, fixing the constant pressure at 1 atm, to determine the equilibrium electrode-electrode distance. At this stage the electrodes act as pistons, and the system box is allowed to elongate or compress along the z direction only. After equilibrating the electrode positions, we run a short N V T simulation to allow the system to fully equilibrate in the conditions used for production. In total, each equilibration includes at least 1 ns of dynamics. In all the simulations, we use a 1 fs time step and the temperature is fixed at 298 K using a Nosé-Hoover chain thermostat [7], with 5 chained thermostats having all the same time constant of 1 ps.

All simulations are performed using the molecular dynamics code Metalwalls [8].

B. Equilibrium and nonequilibrium simulations

The simulations at equilibrium conditions with a voltage ∆ψ = 0 V are conducted for at least 10 ns for production. During the simulations, we sample both the total electrode charge Q and the components of the total dipole moment at each time step. Molecular configurations are sampled each picosecond. From these equilibrium trajectories, we also In order to check the effect of the lateral dimensions of the simulation box L x and L y , we perform a supplemental simulation doubling the dimensions in the two directions, for the system with an electrode-electrode separation d = 4.94 nm. The electrode atoms are replicated in the two directions and the number of water is four times the original amount.

Fig. S1 compares the total charge autocorrelation function for the original and larger systems, scaled by the lateral area. The good agreement between the two shows that the original system is already sufficiently large.

D. Effect of the thermostat

Since our analysis involves the study of dynamical observables, one should consider the possible influence of the thermostat and its parameterization on the frequency-response of the system. We assess this issue in two ways. Firstly, we compare our admittance and impedance results obtained in the canonical ensemble (N V T , with thermostat, using a time constant of 1 ps as in the main text) and an average over microcanonical (N V E, no thermostat) trajectories starting from independent initial conditions generated in the N V T ensemble. Secondly, we investigate the possible effect of the time constant for the coupling of the system with the thermostat, by peforming simulations for a range of values between 0.5 and 8 ps (with the same constant for the 5 elements of the thermostat chain in each case). For both cases, we compare the total charge autocorrelation functions, the admittance with Bode plots and the impedance in Nyquist plots. In order to reduce the computational effort, we performed systematic tests only for the smallest system (d = 2.51 nm), but the conclusions can be extended to the other cases.

Fig. S2 and S3 report the comparison between the electrochemical responses in the N V T and N V E ensemble and the effect of the time constant of the thermostat for the simulations in the N V T ensemble, respectively. The differences are not statististically significant and we can conclude that the presence of a thermostat doesn't affects the dynamics of the system (and corresponding electrochemical observables) in all the considered cases.

III. STATIC PROPERTIES A. Total charge distributions

The total charge histograms are shown in Fig. S4 for all the systems we considered, at 0 V.

All the histograms present a Gaussian shape, with the variance being inversely proportional to the electrode-electrode distance.

FIG. S4. Distribution of total electrode charge for all the systems we considered, at 0 V. The lines are the results from MD simulations, the dots represent Gaussian probability density functions with zero mean and the variance set to be ⟨Q 2 ⟩ as calculated from MD results.

B. Equilibrium density profiles

The equilibrium water density profiles along the z-direction are reported in Fig S5 . All the systems feature the same layered water structure close to the electrode surface. In general, the density profiles differ only in terms of width of the plateau region, which is trivially wider for larger systems, except for the smallest system we considered. In the latter case, the density profile does not reach a plateau, but the layered structure extends through the whole width of the system.

C. Dielectric profile

In order to characterize the local dielectric properties of our systems, we computed the dielectric profiles ε -1 ⊥ (z) in the confining direction (perpendicularly to the electrode planes). We estimated the profiles using the approach described by Schlaich et al. described in Ref. 9.

The results of our analysis are reported in Fig. S6.

All the systems show the same identical charge density profile, shown in Fig. S6a. Based on this observation, we consider the dielectric profile only for the case with d = 4.94 nm without loss of generality. The dielectric profile for interfacial water shows an oscillating behaviour close to the gold surface, as expected. From the dielectric profile, we calculated the position of the Dielectric Dividing Surface (DDS) using the following definition [9]:

z DDS = z v + z l zv ε -1 ⊥ (z l ) -ε -1 ⊥ (z) ε -1 ⊥ (z l ) -ε -1 ⊥ (z v ) dz (6)
where z v is the position of the electrode surface (z = 0 in our reference frame), z l is a position in the water phase. In practice, we calculated ε -1 ⊥ (z l ) as the average value of ε -1 ⊥ (z) in the profile plateau. In our case, we use z DDS to determine the effective water slab width d DDS = d -2w DDS , where w DDS = z DDS -z v ≈ 1.27 Å. The water slab width is used to highlight the scaling relations between the admittance values in systems with varying d. We also define a coarse-grained version of the dielectric profile using a step function based on the position of the DDS [10]

ε ⊥ (z) =      1, for z ≤ z DDS ε w , for z > z DDS , (7)
where ε w is the dielectric constant of bulk water, which corresponds to 70.7 in the specific case of SPC/E water. This simple model is used to rationalize the capacitance results presented in the main text.

IV. CALCULATION OF THE IMPEDANCE FROM MD SIMULATIONS A. Windowing

Time autocorrelation functions estimated from MD results suffer from the presence of noise at large t values, because of the limited sampling. This effect impacts negatively on the numerical calculation of the Laplace transform for the estimation of admittance/impedance.

In particular, we found this noise to particularly deteriorate the results at high frequency.

To suppress this effect, we treated the charge autocorrelation functions by multiplying them with sigmoidal window functions having the following form

W (t, ϵ, τ ) = 1 1 + e ϵ(t-τ ) , (8)
with ϵ and τ two tunable parameters which were adjusted by hand in order to remove the long-time noise while minimizing the impact on the overall shape of the autocorrelation functions. It is clear that this kind of restriction of the original signal suppresses all the slow modes whose period is longer than the region selected by the window function. However, we tuned the window parameters such that only the most significant part of the autocorrelation function was taken into account. A comparison between the original and windowed versions of the autocorrelation functions is reported in Fig. S7 for all the systems. Once the parameters have been adjusted for each system, the same windowing functions were used to treat every data segment used for the estimation of ⟨δQ(0)δQ(t)⟩ and admittance/impedance.

B. Filon-Lagrange integration

The estimation of admittance/impedance involves the numerical calculation of the Fourier-Laplace transform of the electrode charge autocorrelation function ⟨δQ(0)δQ(t)⟩.

For the numerical calculation of the Laplace transforms, and of highly oscillatory integrals in general, general-purpose quadrature methods such as the trapezium or Simpson's rule, may fail in terms of numerical stability, especially at high frequencies. In principle, one could choose a sufficiently small sampling interval until the integrals would converge to a numerically stable solution, but this would require harder computational efforts and the storage of larger amounts of data. A simple but very effective strategy was introduced by Filon, and it is based on representing the integrand using interpolators, for which the integrals can be calculated analytically [11]. In the specific case of the Fourier-Laplace transform of a signal f (t), we can write the following approximation

f (ω) = ∞ 0 f (t)e -iωt dt ≈ N -1 j=0 t kj+k t kj v j (t)e -iωt dt, (9)
where the full integral has been substituted with a sum of integrals over consecutive, nonoverlapping segments (each one embedding k + 1 samples of f (t)), and the signal has been substituted by a function which interpolates the points sampled from f (t). In our case, we chose to interpolate our points using Lagrange polynomials of order 2 (i.e. in the form of at 2 + bt + c) defined on consecutive triplets of data points. With this in mind, Eq. 9 can be rewritten as a sum over consecutive triplets

f (ω) ≈ N -1 j=0 t 2j+2 t 2j v j (t)e -iωt dt = n-1 j=0 t 2j+2 t 2j a j t 2 e -iωt dt + t 2j+2 t 2j b j te -iωt dt + t 2j+2 t 2j c j e -iωt dt, (10)
where the coefficients a j , b j , c j are determined from the Lagrange interpolation of f (t 2j ),

f (t 2j+1 ), f (t 2j+2
). Now the integrals in Eq. 10 can be computed analytically without any numerical error in the quadrature. In this method, we only rely on the assumption that the interpolating functions are sufficiently accurate to represent the original signal. In practice, we first determine the polynomial coefficients by Lagrange interpolation of ⟨δQ(0)δQ(t)⟩ considering consecutive triplets, then use Eq. 10 for the calculation of the Fourier-Laplace transform of ⟨δQ(0)δQ(t)⟩ in the estimation of the admittance. We employ a similar formula for the calculation of the inverse Fourier-Laplace transform, in the analysis of the time response of equivalent circuit models (see next section).

Z(ω) = iωL + R b + 2 iωC el,1 -ω 2 C el,1 R 1 C 1 iωC el,1 R 1 + iωR 1 C 1 + 1 + iωC el,2 -ω 2 C el,2 R 2 C 2 iωC el,2 R 2 + iωR 2 C 2 + 1 -1 . (11)

B. Fitting

We fitted the EC models to admittance and impedance, as calculated from MD results.

Our fitting procedure is entirely defined in the frequency domain and takes both the admittance and impedance into account simultaneously. For fitting a complex function, our approach is to minimize the euclidean distance between the model and the MD results in the complex plane. In the complex plane the distance of a point A to a target value A * is defined as

∥A -A * ∥ = (ℜ[A] -ℜ[A * ]) 2 + (ℑ[A] -ℑ[A * ]) 2 1/2 , (12)
where ℜ and ℑ are the real and imaginary parts, respectively. In order to fit the MD results, we minimize the following loss function

L(x) = w Z ω ∥Z(ω, x) -Z * (ω)∥ ω ∥Z * (ω)∥ + w Y ω ∥Y (ω, x) -Y * (ω)∥ ω ∥Y * (ω)∥ , ( 13
)
where x is a vector containing the EC parameters, Z(ω, x) and Y (ω, x) are the EC impedance and admittance, respectively, and w Z and w Y are weights used to tune the focus on the impedance or on the admittance during the fitting. The weights are manually optimized for each system in order to find a good compromise between the fitting of admittance and impedance. The optimal parameters for the EC models are found as x opt = arg min L(x).

C. Comparison with simpler EC models

The MD impedance results can be fitted with many possible circuit models. We found our EC to accurately represent all the features present in the impedance/admittance spectra, but one might question what happens if we use a simpler model to fit the MD results. In Fig. S9 we compare different EC models used to fit the impedance results, reported in a Nyquist plot. Note that each EC model is fitted independently, using the loss function introduced in the previous section.

As expected, simplifying the original model, labelled as 'EC', reduces the accuracy of the impedance fit. Specifically, removing the inductor and resistor ('EC 2') results in the disappearance of the negative branch at high frequencies, which corresponds to neglecting the oscillations present in the total charge autocorrelation function at short times. In 'EC 3', we further remove one of the two branches from each side of the equivalent circuit. In this case, we obtain a single semicircle in the Nyquist plot, which poorly represents the impedance behaviour at intermediate frequencies.

D. Electrochemical results

Step response in time domain

The equilibrium ⟨δQ(0)δQ(t)⟩ function can be compared with the step response of the EC models. This mainly serves as a validation of the fitted models in the time domain. The charge autocorrelation function is related to the the admittance via

⟨δQ(0)δQ(t)⟩ = L -1 Y (ω) βω 2 , (14)
where L -1 is the inverse Fourier-Laplace transform. Following the approach used by Straube et al. for the definition of the inverse transform [12] this equation can be written in a more explicit form

⟨δQ(0)δQ(t)⟩ = 1 πβ ∞ -∞ ℜ[Y (ω)] ω 2 e -iωt dω. (15)
We numerically calculated the inverse transform using the Filon-Lagrange method explained

Low frequency resistance and cell capacitance

The total resistance at low frequency can be calculated from the real part of the impedance From the equivalent circuit models, the capacitance of the cell can be calculated as

R ω→0 = ℜ[Z(ω → 0)].
C = (C el,1 + C el,2 )/2.
In Fig. S11b we compare the equivalent circuit capacitance with the differential capacitance calculated from the variance of the total charge (β⟨δQ 2 ⟩, from Ref. 13), which corresponds to the capacitance calculated under the Born-Oppenheimer approximation. Overall, we observe a decrease in the cell capacitance as we increase the electrode-electrode separation, as expected for macroscopic capacitors. In all cases, the EC results are in good agreement with the capacitance results from MD simulations.

τ = RC = αε 0 ε w d d + 2(ε w -1)w DDS = τ ∞ d d + 2(ε w -1)w DDS (16)
with τ ∞ the value extrapolated for d → ∞. Fig. S12 shows that this scaling rather accurately describes the evolution of the relaxation time with d. Fitting the results with Eq. 16, we obtain τ ∞ ≈ 8.6 ps. This value is close to the characteristic time corresponding to the Debye relaxation in bulk water (τ D ≈ 9 ps) that can be interpreted as arising from the migration of orientational defects in the H-bond network [14].

Individual parameters

The trends of the equivalent circuit parameters are reported in Fig. S13. In general, we note that all the capacitances tend to decrease as we increase the electrode-electrode separation, as it is expected from macroscopic theory. In contrast, the resistances tend to increase. This is also expected, since wider dielectric slabs are associated with larger values of resistances. We finally note a clear linear scaling the inductance as a function of d.

VI. TOTAL DIPOLE MOMENT

A. Relation with the electrode charge

For a given configuration, the total dipole moment M of the system is

M tot (t) = i∈water q i r i (t) + j∈electrode q j (t)r j , (17)
where the first sum runs over the mobile atoms of the liquid with fixed partial charges, while the second runs over the fixed electrode atoms with fluctuating charges, with q and r denoting the atomic charges and positions, respectively. Note that for the calculation of M, we consider the unwrapped coordinates, to avoid artefacts due to periodic boundary conditions. From the symmetry of the system, one should consider the components of the dipole moment parallel (∥) and perpendicular (⊥) components to the electrode surfaces.

With our choice of axes, the latter simply corresponds to the z component, and we estimate the autocorrelation function (ACF) of the parallel component as the average of the ACF of the x and y components.

Separating the perpendicular component into the water and electrode dipoles, and separating the latter into the contributions of the two electrodes, we have:

M ⊥ tot = M ⊥ wat + M ⊥ el = M ⊥ wat + Q 1 (t)z Q,1 (t) + Q 2 (t)z Q,2 (t) (18)
where in the second line we introduced the total charge of each electrode, as well as the charge-weighted positions of each electrode (which are well defined except when the electrode charge vanishes). Due to the global electroneutrality of the system, at each time step

Q 1 (t) = -Q 2 (t) = Q(t), so that M ⊥ tot = M ⊥ wat -Q(t)d Q (t)
, where we introduced the difference

d Q (t) = z Q,2 (t) -z Q,1 (t).
For a potential difference between the electrodes of 0 V (shortcut conditions), the total dipole in the direction perpendicular to the electrodes vanishes, and we obtain, for each configuration,

M ⊥ wat (t) = Q(t)d Q (t) . (19)
Since the electrodes are treated as perfect metal (for the description of screening within the metal with fluctuating charges in classical MD simulations, see Ref. 15), we expect the charge to be localized mainly in the first atomic plane in contact with the liquid, regardless of the charge distribution within that first plane. Therefore, one can approximate d Q (t) by the time-independent distance d between the first atomic planes of both electrodes. Fig. S14 demonstrates that the resulting proportionality between δM z = M ⊥ wat -⟨M ⊥ wat ⟩ = M z and Q(t) is very well satisfied and that the effective distance d eff obtained from a linear fit of the correlation corresponds to within less than 1% to the the actual distance d between the first atomic planes of both electrodes. This figure shows in contrast no correlation between the total charge of the electrodes and the parallel components of the total water dipole, as expected.

Following Eq. 19 and the observation that d Q (t) ≈ d, one expects that the perpendicular component of the total water dipole moment follows the same dynamics as the total electrode charges. Fig. S15 shows an almost perfect match between the ACFs of Q and M ⊥ wat , while ⟨δM

B. Contributions to the total dipole moment

We run supplemental simulations for the system with d = 4.9 nm, where we split the contributions to the total dipole moment in the perpendicular direction into that of the first layer of water adsorbed on both electrodes (adlayer, m adl ) and that of other water molecules (middle, m mid ). In practice, we divide the system in three contiguous regions along the z direction (with z = 0 and z = d corresponding to the atomic planes of two electrodes in contact with the electrolyte):

• the range {0 < z ≤ z w } denotes the adlayer region of the electrolyte, which fully includes the first density peak of water (we use z w = 4.0 Å);

• the range {z w < z ≤ d -z w } denotes the central region of the liquid (we also refer to this region as bulk-like, even though it also contains some interfacial molecules beyond the adlayer);

• the range {d -z w < z ≤ d} denotes the second adlayer region.

The assignment of each molecule to a given region is done according to the position of its oxygen atom. For each region (adl. or mid.), the polarization is computed as

m region = Nw i=1 z i,O ∈region [q O z i,O + q H (z i,H 1 + z i,H 2 )] , (20)
where the sum runs over all the molecules i with their oxygen atom in the adlayer (resp. middle) region, q is the atomic charge and the labels refer to each element. The dynamics of the two polarization contributions is analyzed by considering the autocorrelation functions ⟨m adl (0)m adl (t)⟩ and ⟨m mid (0)m mid (t)⟩, as well as the cross-correlation function ⟨m adl (0)m mid (t)⟩. For simplicity, we only consider one of the two adlayers in our analysis. ized by their initial value), which all display an oscillatory behaviour at short times that is particularly prominent for the autocorrelation of m adl and its cross-correlation with m mid .

Instead, the oscillations in the autocorrelations of the total polarization and that of the middle region are more damped and are remarkably similar. This is because the number of water molecules residing in the middle part is much larger than the ones in the adlayer (1998 against 161, on average), so that the middle region polarization is almost representative of the global polarization behaviour. This imbalance is also reflected in the magnitude of the polarization fluctuations, ⟨δm 2 adl ⟩ =2.0 e 2 Å2 and ⟨δm 2 mid ⟩ =89.2 e 2 Å2 , respectively. These results already suggest that the fluctuations of the total polarization are dominated by water molecules beyond the adlayers.

In addition, the polarization of the adlayer water molecules oscillates at a higher frequency than the other water molecules. The cross-correlation function ⟨m adl (0)m mid (t)⟩ displays a frequency similar to the autocorrelation of m adl , but its phase appears to be shifted by ≈ π.

The oscillatory behaviour of the total dipole moment and bulk-like polarization is clearly present but due to the stronger damping is not very evident and a deeper investigation is required to treat this feature quantitatively, as discussed below. Note that, in all cases, the time scale of oscillations is much shorter than the characteristic adsorption lifetime of water molecules (≈ 14 ps [17]), thus they are not due to water molecules crossing the boundary between the adlayer and middle regions.

In order to further analyze the oscillatory part, Fig. S16b shows the function ϕ(t) defined as the negative second time derivative of the partial polarization correlation functions, which highlights the change in curvature of the time correlation functions. All the curves display a behaviour reminiscent of a damped oscillator, which we fit with exponentially damped cos (or sin in the case of the cross-correlation) functions to determine the main oscillation frequency. The adlayer water molecules resonate at 1.88×10 14 ± 3.9 × 10 11 rad/s, while those in the middle region resonate with a frequency of 1.33×10 14 ± 1.4 × 10 12 rad/s, matching well with that of the peak in the magnitude of admittance (1.33 × 10 14 rad/s for this specific system). This shift is further illustrated in terms of power spectral density in the inset of Fig. S16b, where we highlight the peak frequency for each polarization contribution.

Our results on the oscillation frequencies are in agreement with the ones obtained for SPC/E water confined between neutral walls [16]. Thus, we follow the same interpretation and assign the oscillations to the librational modes (rotations of water molecules around a fixed orientation). Libration contributes to polarization fluctuations both in the adlayer and in bulk-like water, although the latter is the dominant contribution to the the total polarization. The origin of the shift to higher frequencies of the librational peak at the interface could be due to the stiffer rotational potential felt by water molecules at the adlayer, as already observed experimentally for water/air interfaces [18]. In fact, as seen in a previous study [19], the water molecules in the adlayer on top of the gold surface tend to form a (mostly) flat layer (see Fig. S17). Fig. S17a further indicates that water molecules in contact with the first adlayer exhibit dangling bonds, which may modify their dynamics with respect to bulk water. Performing the same analysis by selecting only the water molecules in the second layer (not shown), we find indeed a peak frequency of 1.65 × 10 14 ± 5.8 × 10 11 rad/s, which is intermediate between that for the adlayer and for "bulk-like" water. Nevertheless, these interfacial modifications of the libration frequency do not significantly impact the total polarization fluctuations.

Overall, the present analysis clearly shows that the prominent peak in the admittance magnitude at relatively high frequencies mainly arises from the librational dynamics of water molecules that reside beyond the adlayer, even though the latter also contribute to a lesser extent, with a relative weight which depends on the distance between the electrodes.

VII. EFFECT OF DISSOLVED IONS

Although being beyond the scope of this work, the contribution of dissolved ions to the electrochemical properties is particularly relevant. As mentioned in the main text, we noted that the presence of ions introduces additional features to the charge dynamics, and, in particular, longer time scales. A first comparison with the case of pure water reported in the present work can be made using the results obtained for a very similar system, using the same microscopic model and a distance d = 4.96 nm, but with a 1 M NaCl solution, reported by Scalfi and coworkers in a previous study on the effect of screening inside the metal on interfacial properties [15]. While this study did not consider the electrochemical impedance, results were shown for the autocorrelation function of the total electrode charge (see Fig. S18). The comparison between pure water and salt solution confirms that the presence of ions leads to additional features, including an increase in the initial value (corresponding to a larger capacitance) and a slow decay at long times, which would translate into additional features at low-frequency in the impedance spectrum. However, the noise at long time in the data reported from Ref. 15 prevents the use of our method to evaluate the frequency-dependent impedance without additional data.

FIG. 1 .

 1 FIG. 1. (a) Simulated system for an inter-electrode distance d = 4.94 nm, with the color of gold electrode atoms indicating their instantaneous charge, and water hydrogen and oxygen atoms shown in white and red, respectively. (b) Differential capacitance divided by the lateral area of the electrodes, A el , as a function of d, with the contributions C BO diff arising from the thermal fluctuations of the liquid and C empty diff arising from those suppressed in the Born-Oppenheimer sampling of the charges, which corresponds to the empty capacitor. Also shown are the predictions of continuum electrostatics without (green dashed line) and considering (symbols) the Dielectric Dividing Surfaces (see text; the error bars on the capacitance arise from that on the positions of the DDS).

FIG. 2 .

 2 FIG. 2. (a) Electrode charge auto-correlation function, ⟨δQ(0)δQ(t)⟩, divided by A el for the 4 considered interelectrode distances; the shaded areas represent the standard error associated with each curve. (b) Response of the charge Q divided by A el (colored solid lines) to a voltage change from ∆Ψ = 0 V to 1 V (black dashed line). For each system, the results are averaged over 10 trajectories. (c) Comparison of the responses with the equilibrium charge ACFs, normalized as (ft -f∞)/(f0 -f∞), where f is either Q recorded during the charging (symbols), or ⟨δQ(0)δQ(t)⟩ (solid lines).

FIG. 3 .

 3 FIG. 3. Nyquist plot of the impedance Z(ω) in the complex plane, parameterized by the frequency ω, for the 4 considered inter-electrode distances (colors). Results from MD simulations using Eq. 1 (symbols) are shown with the predictions of the corresponding equivalent circuit models (solid lines, see text and SM Sec. V [60]).

FIG. 4 .

 4 FIG. 4. Norm of the complex admittance obtained via Eq. 1 for the 4 considered inter-electrode distances (colors). The MD simulation results (solid lines) are shown multiplied by the distance dDDS between the two Dielectric Dividing Surfaces and divided by A el to highlight the scalings with system size. The dashed lines indicate the prediction at low frequency of a mono-exponential model using the capacitance C and the characteristic time τ defined as the integral of the normalized charge ACF. The inset shows τ as a function of d, with a fit corresponding to the DDS picture (see text).

  Each system includes two confining electrodes, separated by different distances d = 2.51, 4.94, 9.76 and 19.42 nm for 1080, 2160, 4320, and 8640 water molecules, respectively. Each electrode consists of 1620 atoms fixed on an FCC lattice (9 × 9 × 5 unit cells with a lattice parameter of a = 4.07 Å corresponding to the unit cell of gold), resulting in ten atomic planes perpendicular to the z direction and facing the inner part of the system with a (100) plane. For all the systems, the box has lateral dimensions L x = L y = 36.63 Å, and periodic boundary conditions are applied only in the x and y directions.All the atoms in the systems interact via electrostatic interactions, computed using a 2D Ewald sum method taking into account the Gaussian distributions of the electrode atoms[2, 3] and truncated and shifted Lennard-Jones (LJ) potentials. Water molecules are modeled with the SPC/E force field[4] and the LJ parameters for the gold atoms are taken from Ref. 5, with the Lorentz-Berthelot mixing rules. Both the electrodes were treated as metals using the fluctuating charge model, in which each electrode atom is equipped with a Gaussian charge of width w = 0.40 Å, with magnitude being determined for each configuration, using the matrix inversion method[6], to fulfill the constraints of constant potential and of global electroneutrality.

extract ten configurations separated by 1

 1 FIG. S1. Total charge autocorrelation functions for the original system (solid lines), and its scaled version (dots), obtained by doubling L x and L y . The results have been scaled by the respective electrode areas A el to facilitate the comparison. The uncertainty is calculated as the standard deviation, obtained splitting the original data in 10 segments.

FIG. S2 .

 S2 FIG. S2. Comparison of electrochemical properties obtained in the canonical ensemble (N V T , with thermostat) or an average over microcanonical (N V E, no thermostat) trajectories starting from independent initial conditions generated in the N V T ensemble. (a) Autocorrelation function of the total charge; (b) magnitude of admittance as a function of frequency; (c) Nyquist plot of the impedance.

FIG. S5 .

 S5 FIG. S5. Water density profiles along the z-direction. (a) Comparison between the different systems. (b) Zoom on the leftmost part of the profiles.

  FIG. S6. Charge density and dielectric profiles along the z-direction. (a) average charge density profiles for all the systems we considered. (b) dielectric profile for the system with d = 4.94. The dashed vertical line indicates the position of the dielectric dividing surface (DDS). The dotted line represents the coarse-grained version of the dielectric profile, defined as a step function and based on the DDS position. The shaded areas indicate the uncertainties, computed as the standard deviation among 10 different samples drawn from the original trajectory.

FIG. S7 .

 S7 FIG. S7. Comparison between the window functions W (see Eq. 8), the raw total charge autocorrelation functions as calculated from the MD simulations, and the same after being multiplied by W . Each plot reports the inter-electrode distance d, and the parameters ϵ and τ used for the window function.

  FIG. S8. Representation of the equivalent circuit used to fit the MD results.

FIG. S9 .

 S9 FIG. S9. Comparison of different EC models (right) and their impedance fitted on MD results, represented in a Nyquist plot (left). MD results are reported as black dots. The impedance of each EC model is represented as a solid line. The circuit structures are reported on the right, with the colors matching the ones used in the Nyquist plot.

  in section IV B, swapping the frequency and time domain. The numerical results for all the systems are shown in Fig S10. In all the cases, the equivalent circuit behaviour is in excellent agreement with the MD results.

FIG. S10 .

 S10 FIG. S10. Comparison between the charge autocorrelation functions as calculated from MD simulations (MD), and as calculated from the equivalent circuit models (EC), for all the systems.

  In our case, the lowest accessible frequency is ω = 1.26 × 10 11 rad/s. The EC and MD results for this property are shown in Fig. S11a. The equivalent circuit results are in excellent agreement with MD, and both the datasets display a linear trend R ∝ d, crossing 0 resistance for d = 0. This result is in agreement with the macroscopic linear relation between the resistance and the width of a resistor.

  FIG. S11. (a) Total resistance at low frequency as a function of the electrode-electrode distance, as calculated from MD (symbols) and EC models (solid line). The results are compared with a linear fit of the MD results to show the extrapolation at d = 0. (b) Comparison between cell capacitance calculated from MD using the electrode charge fluctuations β⟨δQ 2 ⟩ (symbols) and from EC models (solid line) for all the considered inter-electrode distances.

3 .

 3 Fig.S12shows the relaxation time τ of the charge as a function of the inter-electrode distance. τ can be calculated from the integral of the normalized total charge autocorrelation function. In practice, we perform the integral on the autocorrelation function multiplied by the same window function used for the numerical Laplace transforms to mitigate the effect of the noise at long times. We also estimate τ as the product between the total resistance (estimated from the impedance as ℜ[Z(ω → 0)], as calculated from MD results) and the cell capacitance, calculated as β⟨δQ 2 ⟩. The increase of τ with increasing inter-electrode distance reflects the slower relaxation of ⟨δQ(0)δQ(t)⟩. In all cases, we observe a good agreement between the values of τ obtained as RC or from the charge ACF.

FIG. S12 .

 S12 FIG. S12. Charge relaxation time as a function of the inter-electrode distance. τ is estimated either as the integral of the normalized total charge ACF (blue solid line) or as the RC time (orange dotted line), using the total resistance at low frequency and the cell capacitance. The scaling of Eq. 16 is also shown (dashed green line).

FIG. S13 .

 S13 FIG. S13. Equivalent circuit parameters for the four systems we considered.

FIG. S14 .

 S14 FIG. S14. Correlation plots between the total charge and the components of total dipole moment ofwater, for all the considered systems. The total dipole moment values are reported as fluctuations (instantaneous value subtracted by the ensemble average value). Each row represents a system, identified by the electrode-electrode distance (reported on the left). The three columns report the x, y and z component of M wat , respectively. In the third column we also report the effective electrode-electrode distance found from the slope of a line fitted on the data.

∥

  wat (0)δM ∥ wat (t)⟩ decays more slowly and only marginally depends on d, with a slight decrease in the relaxation time with increasing d, consistently with the results reported in Ref 16.

  FIG. S15. Autocorrelation functions of the dipole moment of the confined water slab, with components M ⊥ wat (solid lines) and M ∥ wat (dashed lines) perpendicular and parallel to the electrode surfaces, and of the electrode charge Q (symbols), for the 4 considered inter-electrode distances (colors). The uncertainty, quantified as the standard error from 10 blocks of each trajectory, is shown as shaded areas for M ⊥ wat and as error bars for Q; it is omitted for M ∥ wat for clarity.

Fig.

  Fig. S16a compares the different partial polarization time correlation functions (normal-

  FIG. S16. Dynamical behaviour of the polarization components (see text for their definition). (a) Comparison between the time auto-and cross-correlation functions of the partial and total dipole moments, normalized by their initial value. (b) Negative second derivative, ϕ(t), of the auto-and cross-correlation functions. The inset in (b) shows the power spectral density of ϕ(t) with the same color coding as the main plot.

  FIG. S17. Representative snapshots of the water/gold interface. (a) Side view. (b) Top view (only water molecules in the adlayer are shown). Gold, hydrogen and oxygen atoms are shown in yellow, white and red, respectively.

  FIG. S18. Electrode charge time autocorrelation function ⟨δQ(0)δQ(t)⟩ in the case of pure water (blue), and 1M NaCl solution (orange), as reported in Ref 15, for the system with d = 4.96 nm.

  1, 2, * 1 Sorbonne Université, CNRS, Physico-chimie des Électrolytes et Nanosystèmes Interfaciaux, PHENIX, F-75005 Paris 2 Réseau sur le Stockage Electrochimique de l'Energie (RS2E), FR CNRS 3459, 80039 Amiens Cedex, France

The authors thank Mathieu Salanne, Iurii Chubak, Minh-Thé Hoang Ngoc and Jeongmin Kim for discussions. This project received funding from the European Research Council under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 863473). The authors acknowledge access to HPC resources from GENCI (grant no. A0110912966).

Frequency-dependent impedance of nanocapacitors from electrode charge fluctuations as a probe of electrolyte dynamics

Supplemental Material

Giovanni Pireddu 1 and Benjamin Rotenberg


The Park-and Loop Routing Problem with Parking Selection

Théo Le Colleter 1 , Dorian Dumez 1,2 Fabien Lehuédé 1 , Olivier Péton 1 1 IMT Atlantique, LS2N, Nantes, France {fabien.lehuede,olivier.peton}@imt-atlantique.fr 2 HEC Montréal, CIRRELT, Montréal, Canada dorian.dumez@hec.ca

Mots-clés : logistique urbaine, recherche à voisinage large, park-and-loop

This paper presents a variant of the vehicle routing problem considering deliveries to customers in cities with a combination of walking and driving. The objective is first to offer a better model for actual delivery problems where drivers are usually assumed to have a single parking location to deliver each customer. Second, we show that serving some customers through walking trips is more efficient in congested areas. We introduce the Park-and-Loop Routing Problem with Parking Selection (PLRP-PS), which extends the work of [START_REF] Cabrera | The doubly open park-and-loop routing problem[END_REF][START_REF] Coindreau | Vehicle routing with transportable resources : Using carpooling and walking for on-site services[END_REF]. We investigate the case where parking locations should be selected among a large set of parking areas in a variant of Large Neighborhood Search (LNS).

We consider a set of customers N with non-unitary demand who have to be served using a homogeneous vehicle fleet K. We consider that each vehicle starts from and return to a depot 0. It is assigned to a driver who will perform deliveries by walking from a parking location to the customers. We consider a large set P of parking locations. Each parking location i ∈ P is associated with a parking time pt i . For each pair of nodes i, j in {0} ∪ P , we assume a driving time t d ij and distance d d ij to be known, likewise a walking time t w ij and distance d w ij are defined for each i, j ∈ P ∪ N .

A vehicle route starts from the depot, visits a sequence of parking locations in P , and returns to the depot. The sum of customers demands served by a vehicle route must not exceed the vehicle capacity and a vehicle route duration must not exceed the day duration. A walking-trip starts from a parking location in P , visits a sequence of customers, and returns to the same parking location. It must not serve customers whose cumulative demand exceeds the walking capacity. The overall walking distance of each driver should not exceed a given daily limit. A driver may perform several walking-trips from a single and/or multiple parking locations.

The PLRP-PS consists of designing vehicle routes and walking-trips, respecting vehicle capacities, maximum route duration and the maximum walking distance of each driver, such that each customer is served with its demand by exactly one driver, lexicographically minimizing the number of vehicles and the sum of driving times, walking times and parking times. This problem is illustrated on Figure 1 To solve the PLRP-PS we introduce the Small and Large Neighborhood Search (SLNS) variant of the LNS metaheuristic. We extend the classical destroy and repair operators to efficiently integrate parking selection. Considering the large number of parking locations that may be available in a city, independently of the number of customers, we focus on comparing different parking selection strategies which can be applied in the evaluation of the customer insertion in a vehicle route. We consider two criteria to select which parking to evaluate : the walking time between the parking and the customer, and the ranking of the parking in the list of parking ordered by non-decreasing distance to the considered customer. We then devise two methods to reduce the set of evaluated parkings : filter a priori the parking set from which each customer can be served (keeping the locations below a given threshold on the chosen criterion), and evaluate the parking in non-decreasing order of the chosen criterion with a given probability to stop after each location evaluation.

On realistic instances, we find that combining walking and driving in urban delivery can save 19% of working time on average compared to the classical vehicle routing approach (i.e. CVRP). We also find that the most efficient parking selection method is to evaluate the parking in order of their rank in the list of parking close to the considered customer and to stop according to a geometric distribution on the aforementioned rank. This method proved to be efficient on PLRP benchmarks with respect to the methods of ([2, 1]). A longer description of this work is under review and available on HAL ( [START_REF] Colleter | Small and Large Neighborhood Search for the Park-and-Loop Routing Problem with Parking Selection[END_REF]).

  , with three vehicle routes (green, grey and blue) starting from the black depot on the left. Walking paths are in red. Red points indicate customers and the dark blue points indicate the selected parking locations. Light blue points indicate unused delivery areas.

FIG. 1 -

 1 FIG.1-Visualization of the best found solution for instance "Nantes-50-1" with 50 customers