Victor Hoffmann

Ilias Nahmed
email: ilias.nahmed9@etu.univ-lorraine.frparisa.rastin@loria.fr

Parisa Rastin

Guénaël Cabanes
email: guenael.cabanes@lipn.univ-paris13.fr

Julien Boisse
email: julien.boisse@univ-lorraine.fr

Ensmn

Deep Learning Approaches for Dynamic Mechanical Analysis of Viscoelastic Fiber Composites

The increased adoption of reinforced polymer (RP) composite materials, driven by eco-design standards, calls for a fine balance between lightness, stiffness, and effective vibration control. These materials are integral to enhancing comfort, safety, and energy efficiency. Dynamic Mechanical Analysis (DMA) characterizes viscoelastic behavior, yet there's a growing interest in using Machine Learning (ML) to expedite the design and understanding of microstructures. In this paper we aims to map microstructures to their mechanical properties using deep neural networks, speeding up the process and allowing for the generation of microstructures from desired properties.

Introduction

Designing strong, lightweight polymer-based composites with accurate mechanical behavior requires micromechanical approaches to describe a heterogeneous material as a set of constituents whose morphology, size distribution, and local mechanical behavior affect global behavior. Under dynamic loading, their viscoelastic properties, which are responsible for damping, must be taken into account during the design phase, as they can promote vibration isolation [START_REF] Treviso | Damping in composite materials: Properties and models[END_REF].

Numerous micromechanical codes are used to calculate the mechanical response to various loads applied to Representative Volume Elements (RVEs) of composite materials. Among these codes, Finite Element (FE) or Fast Fourier Transform (FFT) methods are widely used to compute local stress and strain fields as a function of time, frequency, or temperature under various loading conditions and to average them across RVE's for comparison with experimental observations [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF].

In the context of high performance composite design, Dynamic Mechanical Analysis (DMA) is one of the best ways to study and characterize the viscoelastic behavior of polymer matrix composites. DMA is known to be a privileged tool for the study of materials (especially polymers and rubbers) whose rheological behavior is viscoelastic by nature, i.e. involving irreversible dissipation of mechanical energy into heat. In its frequency mode DMA is based on harmonic steady-state strain excitations applied at ω pulsations to a material sample and recording the corresponding output stress signal. Complex algebra from experimental data gives access to the conservative (storage / glassy) or dissipative (loss / rubbery) modulus or compliance, which are the real and imaginary parts of their complex nature: M * (ω) = M ′ (ω) + jM ′′ (ω) for the modulus. This allows a full dynamic characterization of the material [START_REF] Agbossou | Modelling of the viscoelastic behaviour of amorphous thermoplastic/glass beads composites based on the evaluation of the complex poisson's ratio of the polymer matrix[END_REF].

Characterizing the mechanical properties of numerically modeled composites is a complex challenge, usually requiring numerous calculations and considerable computational time when using conventional Finite Element Method (FEM) or Fast Fourier Transform (FFT) based codes [START_REF] Prakash | Simulation of micromechanical behavior of polycrystals: finite elements versus fast fourier transforms[END_REF][START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF]. With the advent of deep learning methods, a new approach to solving these issues has arisen, especially Convolutional Neural Networks (CNNs). However, any deep learning method still require a training database consisting of RVEs and their corresponding mechanical responses.

In this study, we focus on a first generated database consisting of RVEs of simple 2D organized fiber-reinforced composite microstructures built from a home-grown Matlab code, and their associated mechanical DMA signatures (see Fig. 1) computed with the FFT-based code CraFT-Virtual-DMA [START_REF] André | An FFT solver used for virtual Dynamic Mechanical Analysis experiments: Application to a glassy/amorphous system and to a particulate composite[END_REF][START_REF] Boisse | CraFT virtual DMA • GitLab[END_REF]. We present the testing of a convolutional neural network (CNN) algorithm trained on this first database.

State of the Art

Obtaining physical properties of simulated composite materials is a complex challenge that requires numerous calculations and a significant amount of computational time when using classical methods such as finite elements (FEM) or Fast Fourier Transform (FFT) methods. The emergence of deep learning methods, particularly Convolutional Neural Networks (CNNs), has provided a new approach to addressing these issues. A Convolutional Neural Network (CNN) [START_REF] Shea | An introduction to convolutional neural networks[END_REF] is a type of artificial neural network specifically designed for processing structured grid data, such as images and videos. CNNs are particularly well-suited for tasks involving visual recognition and understanding. They have become the foundation of wide range of applications such as Image Classification [START_REF] Guo | Simple convolutional neural network on image classification[END_REF], Medical Imaging [START_REF] Kayalibay | Cnn-based segmentation of medical imaging data[END_REF], Environmental Monitoring [START_REF] İnik | Cnn hyper-parameter optimization for environmental sound classification[END_REF], Self-driving Cars [START_REF] Nugraha | Towards self-driving car using convolutional neural network and road lane detector[END_REF], Natural Language Processing (NLP) [START_REF] Li | Application of convolutional neural network in natural language processing[END_REF] and many other. Over the last few years, CNNs have demonstrated notable efficiency in tackling these problems.

In 2018, Zijiang Yang & al. [START_REF] Yang | Deep learning approaches for mining structure-property linkages in high contrast composites from simulation datasets[END_REF] implemented 3D CNNs for the first time to model elastic homogenization linkages of 3D images of composite materials. The deep learning model introduced in the paper is utilized to understand the nonlinear relationship between the three-dimensional microstructure of the material and its overall stiffness. The research showcases the capability of this comprehensive approach to predict the stiffness of high-contrast elastic composites, considering a diverse array of microstructures. This model showed remarkable accuracy while requiring minimal computational resources for new assessments.

In 2020, Yixing Wang & al. [START_REF] Wang | Mining structure-property relationships in polymer nanocomposites using data driven finite element analysis and multi-task convolutional neural networks[END_REF] introduce a data-driven and deep learning model to establish a part of the relationship between the structure and properties of polymer nanocomposites. This model enables the exploration of structures and the identification of qualitative relationships between microstructure descriptors and mechanical properties. The authors introduce a novel deep learning method that combines convolutional neural networks with multi-task learning to create quantitative correlations between microstructures and property values.

Other techniques have been implemented to improve or simplify the CNN models used. For instance, A. Mann and S. Kalidindi [START_REF] Mann | Development of a Robust CNN Model for Capturing Microstructure-Property Linkages and Building Property Closures Supporting Material Design[END_REF] worked on a novel CNN architecture with significantly fewer trainable parameters than current benchmarks. It achieves this by avoiding fully connected layers and utilizing 2-point spatial correlations of microstructures as input. This approach enhances robustness and enables cost-effective property closures by exploring the space of valid 2-point spatial correlations, which is a convex hull. It allows the generation of new valid 2-point spatial correlations from existing sets through convex combinations. This work highlights the advantages of using 2-point spatial correlations over voxelated discrete microstructures in current benchmarks. [START_REF] Xu | Learning viscoelasticity models from indirect data using deep neural networks[END_REF][START_REF] Kim | Exploration of optimal microstructure and mechanical properties in continuous microstructure space using a variational autoencoder[END_REF][START_REF] Li | A transfer learning approach for microstructure reconstruction and structure-property predictions[END_REF][START_REF] Xu | Use of machine learning methods in syntactic foam design[END_REF], however, as far as we are aware, none have sought to predict the complete modulus curve encompassing both storage and loss moduli. Consequently, there is no established benchmark for comparing our model, except for the conventional methods employed in mechanical simulations.

Numerous papers have been published in this context

In the next section, we will present the dataset used to validate our algorithm and introduce the proposed algorithm.

Experimentation

Dataset

The input data consists of 15 000 distinct configurations of representative volume elements (RVEs) for unidirectional fiber-reinforced composite materials. All RVEs adhere to periodic boundary conditions. The fibers all have the same radius and vary from 1 to 150 per RVE, with volume fractions ranging from 0.05 (5%) to 0.75 (75%). A total of 100 random configurations have been simulated for each volume fraction (with a step of 0.005) Typically, RVEs are represented in three-dimensional matrices, but since we are studying unidirectional / 2D organized fiber materials, a cross-section of each matrix is sufficient. As a result, they will be represented by two-dimensional matrices of size 256 × 256, where each element is either 0.5 (filler) or -0.5 (matrix). Examples of RVE representations are given in figure 2. Each component (matrix and fibers) follows an isotropic viscoelastic behavior decribed by a 3D Standart Linear Solid (SLS) model with one relaxation time. The local constitutive law, i.e. stress as a function of strain, of each voxel is computed with the following complex linear relation in the frequency (or pulsation) domain

σ(x, jω) = K(x, jω) - 2 3 G(x, jω) Tr [ε(x, jω)] I + 2G(x, jω)ε(x, jω) (1)
where σ(x, jω) and ε(x, jω) are local stress and strain respectively as a function of position x and pulsation ω. K(x, jω) and G(x, jω) are local viscoelastic bulk and shear moduli which are computed as follow

M (x, jω) = M ∞ (x) + jωτ M (x)M 0 (x) 1 + jωτ M (x) with M = K, G (2)
The local material mechanical parameters M ∞ (x), M 0 (x) and τ M (x) are relaxed modulus, unrelaxed modulus and relaxation time respectively.

For each RVE, a mechanical response to an applied mean shear loading ε xy = ε 0 cos(ωt) have been computed using the CraFT Virtual DMA code. Average of the stress responses over the RVEs as a function of the pulsations ω are obtained directly from the simulations. DMA curves, i.e. storage and loss shear moduli, are real and imaginary parts of the complex shear modulus

G * xy (ω) G * xy (ω) = ⟨σ xy (x, jω)⟩ RV E ε 0 = G ′ xy (ω) + jG ′′ xy (ω) (3)
In the context of our study the Mechanical parameters given in table 1 for each material have been chosen to enlighten the two distinctive contributions of the matrix and fibers to the global averaged behavior of RVEs characterized by the DMA curves (Fig. 1).

The output data consists of a set of two dynamic mechanical analysis (DMA) curves for each RVE: the shear storage modulus and the shear loss modulus, both provided in MPa as ) and relaxation time τ K,G have been chosen considering a typical polypropylen material. For the fibers: unrelaxed young modulus (E 0 ) and poisson coefficient (ν 0 ) have been taken from typical glass mechanical properties [START_REF] André | An FFT solver used for virtual Dynamic Mechanical Analysis experiments: Application to a glassy/amorphous system and to a particulate composite[END_REF]. A ratio of 10 between unrelaxed (0) on Relaxed (∞) moduli as well as a relaxation time of 10s have been chosen artificially. For both materials: Poisson coefficient between unrelaxed (0) and relaxed states (∞) have been kept constant, and, bulk (K) and shear (G) moduli have been calculated starting from Young moduli and Poisson coefficients.

K 0 (GPa) K ∞ (GPa) τ K (s) G 0 (GPa) G ∞ (GPa) τ G (
a function of the given angular frequency in rad s -1 . Each curve is represented by a collection of 30 data points provided on a logarithmic scale of angular frequency. The objective of our model is thus to predict these two curves solely based on the representations of the RVEs. It turns out that under shear stress, both the shear storage modulus and the shear loss modulus remain invariant under a 90°rotation and horizontal and vertical symmetries of the RVE.

In order for the model to understand periodic boundary conditions, the data has been augmented with 4 random (2 horizontal and 2 vertical) translations of each RVE. The computer's used RAM (32 GB) did not allow for further data expansion. Regarding the invariance under 90°rotation, horizontal, and vertical symmetry, we opted for an equivariant convolutional neural network [START_REF] Weiler | General E(2)-Equivariant Steerable CNNs[END_REF] instead of data augmentation (cf. 3.2). 70% of the data was used for training the model, 15% was used for model validation, and the remaining 15% was used for model testing. The batch size chosen for our model is 40.

Equivariant Convolutional Neural Networks

Convolutional Neural Networks (CNNs) have gained substantial attention in recent years due to their effective performance in image processing tasks. However, traditional CNNs are not equipped to deal with transformations or symmetries in the input data. This limitation prompted the development of Equivariant Steerable CNNs (ESCNNs) [START_REF] Weiler | General E(2)-Equivariant Steerable CNNs[END_REF], a form of convolutional neural network designed to handle specific transformations, such as rotations and symmetries without the need to see all possible transformed versions of the input data during training.

Let's denote the input to a layer as x, the transformed version of this input as T (x) where T is the transformation (for instance, rotation), and let f denote the layer's operation. Then, the layer is equivariant if there exists a transformation T ′ such that for all inputs x:

f (T (x)) = T ′ (f (x)) (4)
(4) means that applying the transformation to the input and then applying the layer operation is the same as applying the layer operation first and then transforming the output. Such equivariances are achieved through the use of steerable filters. Steerable filters are special kinds of filters whose responses can be computed for any rotated version from a limited set of basis filters. This property gives them the power to detect features irrespective of their orientation in the image, making the network more robust and efficient. Specifically, in the context of ESCNNs, the steerability of these filters helps in achieving rotational equivariance.

Therefore, ESCNNs have allowed us to make our model remain invariant under a 90°rotation and horizontal and vertical symmetries of the RVE without data augmentation. However, this has a computational cost since for each equivariant convolutional block, we have 8 parallel convolutional blocks (1 for each unique combination of 90°rotation and symmetry).

Architecture of the model

The figure 3 shows a scheme of the architecture's model.

As our objective is to concurrently predict both the shear storage modulus and the shear loss modulus, we have opted for a multi-task model. This model consists of a common equivariant convolutional component, which subsequently diverges into two branches: one for the shear storage modulus and another for the shear loss modulus. Each branch is composed of at least one equivariant convolutional block, and several layers of Recurrent Neural Networks (RNNs). A Recurrent Neural Network (RNN) [START_REF] Medsker | Recurrent neural networks[END_REF]is a type of artificial neural network designed for sequential data processing. Despite their utility, basic RNNs have limitations, including the "vanishing gradient" problem, which makes it difficult to capture long-term dependencies. To address this, more advanced RNN variants, such as Long Short-Term Memory (LSTM) [START_REF] Graves | Long Short-Term Memory[END_REF] and Gated Recurrent Unit (GRU) [START_REF] Dey | Gate-variants of gated recurrent unit (gru) neural networks[END_REF], have been developed. These variants are better equipped to capture and propagate information over longer sequences and are widely used in practice for various tasks.

Each convolutional block is composed of an equivariant convolutional layer with a kernel size of 3 × 3 and circular padding, a batch normalization layer to which we apply the ReLU activation function and an Antialiased MaxPooling layer. The purpose of the Antialiased MaxPooling is to apply MaxPooling to the kernel which has undergone a Gaussian filter that slightly blurs the kernel [START_REF] Zhang | Making convolutional networks shift-invariant again[END_REF]. This enables the ESCNN to be invariant to slight translations of the image, allowing it, in conjunction with circular padding, to better comprehend the notions of periodic boundary conditions.

At the end of the convolutional blocks in both branches, a group pooling is applied. Group pooling aggregates the eight different convolutional blocks executed in parallel for each rotation/symmetry combination in order to extract features invariant to such transformations [START_REF] Weiler | General E(2)-Equivariant Steerable CNNs[END_REF].

The two branches are not identical. The branch for the shear storage modulus comprises two equivariant convolutional blocks and three layers of Gated Recurrent Units (GRUs), between which we have applied a Dropout rate of 0.2. It turns out that the shear storage modulus is more sensitive to short-term dependencies and requires more feature extraction than the shear loss modulus. Conversely, the branch for the shear loss modulus has only one equivariant convolutional block and two layers of Long Short-Term Memory (LSTM) units, between which we have implemented a Dropout rate of 0.3. At the end of each RNN block in both branches, a Fully Connected layer is applied to obtain values in R. Indeed, in Pytorch, the output of RNNs falls between 0 and 1, due to a sigmoid activation function that is implicitly applied [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF].

Performance metric

The scale of both the shear storage modulus and the shear loss modulus varies according to the volume fraction of the microstructures given as input. The lower the volume fraction, the lower the minimum and maximum of each modulus. To avoid the model giving too much importance to error minimization on moduli of microstructures with high volume fractions and not enough importance to low volume fractions, a relative loss function is necessary. The Mean Absolute Percentage Error (MAPE) is a classic choice for these types of situations. However, it poses problems for several loss modulus values that are very close to 0 for some low volume fraction cases, thus causing an explosion of this latter. To address this issue, we have opted for the weighted Mean Absolute Percentage Error (wMAPE) which is defined for each batch and each modulus as:

wMAPE = 1 N seq. j∈seq. i∈batch |y ij -ŷij | |y ij | (5)
where seq. stands for sequence, which is the discrete representation of the storage or loss modulus, N seq = 30, y ij is the j th point of the i th sequence in the batch and ŷij is the prediction of y ij . As shown in [START_REF] André | An FFT solver used for virtual Dynamic Mechanical Analysis experiments: Application to a glassy/amorphous system and to a particulate composite[END_REF] for each point of the moduli, we calculate the average absolute difference between the real and the predicted value of the batch, then we divide it by the average real value of the batch. Finally, we take the average batch error of each point of each modulus and sum the two error modulus. Therefore, the performance metric is equal to wMAPE storage + wMAPE loss . In the section 4, the Mean Absolute Error (MAE) will also be provided for reference, which is defined for each batch and each modulus as:

MAE = 1 N seq. j∈seq. 1 N batch i∈batch |y ij -ŷij | (6)
with the same notations as in [START_REF] André | An FFT solver used for virtual Dynamic Mechanical Analysis experiments: Application to a glassy/amorphous system and to a particulate composite[END_REF] and N batch = 40.

Regarding the optimizer, we have chosen Adam, with a learning rate of 0.001 and an L 2 regularization term with α set to 0.0001.

Results

The table 2 shows the MAE and wMAPE of the model for both shear storage modulus and shear loss modulus. Each mean error also includes its standard deviation. As mentioned in 2, to the best of our knowledge, no prior literature exists that endeavors to predict the complete set of Dynamic Mechanical Analysis (DMA) curves for both shear storage modulus and shear loss modulus as functions of frequency. Consequently, we lack a benchmark against which to compare our findings. As depicted in Figure 4, the challenge in predicting the loss modulus primarily lies within the range of intermediate frequencies, with values at the extremes being nearly negligible. Regarding the storage modulus, the difficulty of prediction often arises for high-frequency values. 4: Prediction example of shear storage and loss moduli. Please take note that this concerns the Mean Absolute Percentage Error (MAPE), not the previously used weighted MAPE (wMAPE). Indeed, when there is just one prediction curve compared to a real curve (which is equivalent to a batch size set to 1), the wMAPE is equivalent to the MAPE.

With the purpose of assessing the robustness of our model, we conducted an analysis of its performance on each RVE as a function of the volume fraction. As depicted in Figure 5, the model demonstrates robustness in predicting storage moduli for RVEs which volume fraction ranges from 0.07 to 0.73. However, the model's robustness appears to diminish when predicting loss moduli for RVEs with volume fractions near 0.3 or 0.5. Additionally, the uncertainty of the loss modulus is higher than the uncertainty of the storage modulus. The MAPE of the loss modulus for volume fractions between 0.05 and 0.1 is challenging to interpret, because the loss modulus exhibits values very close to 0, which leads to a divergence in the MAPE values. 4, please take note that this concerns the Mean Absolute Percentage Error (MAPE), not the previously used weighted MAPE (wMAPE). To generate such curves, the model was assessed on each RVE separately, employing a batch size of 1, as opposed to the previous procedure.

Figure 1 :

 1 Figure 1: Modeled composite and corresponding storage and loss shear moduli. left: Representative Volume Element (RVE) of a 2D organized fiber reinforced composite with volume fractions 20% fibers (in red) -80% matrix (transparent) satisfying periodic boundary conditions (PBC). Right: Storage G ′ xy (blue) and loss G ′′ xy (red) shear moduli as a function of pulsation ω obtained from the computed DMA response of the RVE to a harmonic average shear strain load ε xy = ε 0 cos(ωt).

Figure 2 :

 2 Figure 2: Examples of 2D Representative Volume Elements (RVEs) with volume fraction of 20%, 70% and 50% fibers (from left to right). RVEs are depicted as 256 × 256 matrices (as shown on the right), where each element takes on a value of either 0.5 (filler) or -0.5 (matrix). All RVEs respect periodic boundary conditions (PBC).

Figure 3 :

 3 Figure 3: Architecture of the multi-task model. It consists of a common equivariant convolutional component, which divides itself into two branches: one for the shear storage modulus and another for the shear loss modulus.

Figure 5 :

 5 Figure 5: MAPE as a function of volume fraction. The average value of MAPE for each volume fraction has been plotted as a solid line, and the standard deviation of this value is depicted by the interval around the mean value. As mentioned in figure4, please take note that this concerns the Mean Absolute Percentage Error (MAPE), not the previously used weighted MAPE (wMAPE). To generate such curves, the model was assessed on each RVE separately, employing a batch size of 1, as opposed to the previous procedure.

Table 1 :

 1 Mechanical parameters for matrix and fibers. For the matrix: unrelaxed Young modulus (E 0 ), relaxed Young modulus (E ∞ ), poisson coefficient (ν 0

	s)

Table 2 :

 2 wMAPE and MAE metrics of the model evaluated on the test set. The model has been trained on the wMAPE, the MAE is provided for reference.

		Shear Storage Modulus Shear Loss Modulus
	wMAPE	0.83% +/-0.13%	1.85% +/-0.25%
	MAE	0.0134 +/-0.0027	0.0025 +/-0.0005


George Bokun

Dung Di Caprio

Myroslav Holovko

Vyacheslav Vikhrenko

The system of mobile ions in lattice models: Screening effects, thermodynamic and electrophysical properties

Keywords: charge screening, lattice fluid model, electric double layer, differential capacitance, mean spherical approximation

The lattice fluid model of the system with short range and long range Coulomb interactions is suggested. In the framework of the collective variables method, the screening of the Coulomb interactions in the bulk is considered. It is shown that the Debye length includes additional concentration dependence inversely proportional to the square root of the mean concentration of vacant sites like what is known at the plane boundary. The Coulomb interaction contribution to the free energy of the system is calculated in the approach close to the mean spherical approximation and is given in an analytical form.

The influence of the variation of the crystal field near the system boundary on the structure and characteristics of the electric double layer is investigated. As compared to the system with equal crystal potentials at the lattice sites throughout the system the pronounced difference for the electric capacitance appears at low absolute values of the surface potential and it is more pronounced for negative electric potentials. The capacitance diverges as the potential values at which the electric field tends to zero and attains negative values in regions of the surface potentials depending on their polarity and values of the surface crystal potential. Negative values of the capacitance may indicate the thermodynamic instability of the system that can result from neglecting the short range interaction contribution.

Introduction

Solids with high ability of ionic transport are widely used in devices for energy transformation and storage [1][2][3][4][5], electrochemical analysis of compounds, light transformation [6][7][8], etc. The examples of such solids are superionic conductors, solid electrolytes, intercalation compounds, to name a few. The most prominent feature of such solids is the presence of the subsystem of mobile ions that can move inside of the matrix created by the hosting one, which creates the neutralizing background for moving ions is frequently considered as a continuum that is subject to hydrodynamic equations of motion [9][10][11][12], and can be regarded as a specific fluid or liquid. On the other hand, the mobile ions move in the potential field created by the host subsystem that can be imagined as the host subsystem potential relief. The mobile ions perform thermally activated hops between the potential relief minima that in many cases are distributed like sites of a periodic lattice. The hops can be performed to vacant sites. The concept of vacancies first introduced for explaining diffusion in ionic solids and simultaneously extended to liquids [13,[START_REF] Frenkel | Kinetic theory of liquids[END_REF] afterwards was widely used for investigating numerous phenomena in physics of condensed matter including e.g. hole transport in semiconductors.

In the first approximation, it can be assumed that the role of the host is just in the creation of the potential relief, which does not change under the influence of mobile ions movement. Nevertheless, real solid electrolytes have a complex structure. For example, ceramic electrolytes consist of grains and intergrain layers. Numerous models for describing intergrain regions were suggested [START_REF] Maier | [END_REF][16][17][18][19][20][21][22]. They include space charge distribution due to impurities segregation at the grain boundaries [START_REF] Maier | [END_REF][16][17], Schottky-type potential barriers at the grain boundaries [18][19][20], a linear diffusion model [21,22] and others. The Monte Carlo simulation shows an unexpected charge distribution near the grain boundary modeled by potential barriers [23]. Like a grain boundary, the system boundary violates its symmetry and can result in variation of the energy well depths on the potential relief close to the boundary. Below the one component subsystem of mobile ions in the bulk and at a system boundary is considered in the framework of lattice models.

Lattice fluid models are widely used for describing the subsystem of mobile ions in concentrated electrolytes, ionic liquids, solid electrolytes and superionic conductors [24][25][26][START_REF] Dutta | Proc. Nat. Inst. Sci. India[END_REF][START_REF] Kröger | The Chemistry of Imperfect Crystals[END_REF][START_REF] Yu | [END_REF][30][31][32][33][34]. The main efforts have been concentrated on investigating in the mean field approximation the structure and properties of the electric double layers arising at the influence of an applied external potential. The lattice models permit accounting of spatial restrictions for particles distribution due to their size. These restrictions remind the requirements of the Fermi-Dirac statistics for the occupation numbers. Thus, Fermi-type distributions for ions in the electric double layers were deduced in particular forms, which were model dependent. The less known fact is that the Debye length in the region of the electric double layer contains the difference between the maximal and current occupation numbers in the bulk, as a correction for its concentration dependence [30].

Less attention has been paid to bulk properties of electrolytes especially beyond the mean field approximations although they are important for investigating their electro-physical characteristics and possible phase transitions. The method of collective variables [START_REF] Yukhnovskii | Statistical Theory of Classical Equilibrium Systems[END_REF][START_REF] Golovko | Approaches to the many-body theory of dense ion-dipole plasma. Application to ionic solvation[END_REF][START_REF] Yukhnovskii | Phase transition of the second order: Collective variables method[END_REF] is a promising approach for accounting of correlations in the homogeneous bulk as well as in the non-homogeneous electrolytes of different types. On the other hand, the statistical method of conditional distributions [START_REF] Rott | Statistical Theory of Molecular Systems[END_REF][START_REF] Rott | [END_REF][40] can help rationalize accounting for the restrictions on spatial distribution of ions. In the framework of these methods, the models of interface boundaries can be constructed for accounting of e.g. specific adsorption of ions.

Charge screening in the bulk

Let us consider a lattice fluid model of N charged particles (ions) that can move over M>N sites of the lattice created by the host system, which plays a role of neutralizing background as well. The particles of the host system are fixed at their given positions. As examples, we can mention superionic conductors (e.g., AgJ), where cations can move on the background of anions, or Yttrium stabilized Zirconium (YSZ), where the vacancies in the oxygen sublattice can move on the background of the host YSZ. For the moving particles, the host system creates the potential energy relief with minima of -u i at lattice sites (the subscript i indicates the position of the corresponding lattice site). For a homogeneous system all the minima are equal u i = u.

The distribution of particles over the system volume is described by the distribution functions. The two-particle (binary) distribution function F 2 (q i ,q j ) can be represented [START_REF] Rott | Statistical Theory of Molecular Systems[END_REF][START_REF] Rott | [END_REF][40] through the one-particle (unary) functions F 1 (q i ) and the binary correlation function g(q i ,q j ) 2 1 1

( , ) ( ) ( ) ( , ),

i j i j i j F F F g  q q q q q q (1)
where q i ,q j determine positions of particles near the corresponding lattice sites. For solid electrolytes, the unary distribution function is a periodic function with sharp extremes at the lattice sites. In the method of collective variables [START_REF] Yukhnovskii | Statistical Theory of Classical Equilibrium Systems[END_REF] it was shown that in the approximation corresponding to the Debye description of ionic systems, the correlation function can be written through the dimensionless screened Coulomb potential ( , )

ij u qq : B ( , ) exp[ ( , )](1 ...), 1/ ( ), i j i j g u k T      q q q q (2) 2 B D B 0 ( , ) ( , ) ( / ) exp( r), | |, 1/ , / 4 , i j i j i j u u r r r r r e            q q q q q q (3)
where u and u are the Coulomb and dimensionless Coulomb interaction potentials, respectively, k B the Boltzmann constant, T the absolute temperature, r B and r D are the Bjerrum and Debye lengths, correspondingly, e is the charge of particles,  and  0 are the electric constant and the medium dielectric constant, respectively.

The screening of the Coulomb interaction is determined by the binary distribution function, and in the framework of the method of collective variables, it is shown [START_REF] Yukhnovskii | Statistical Theory of Classical Equilibrium Systems[END_REF] that the Fourier transformation of the screened potential can be written as

2 1 ( / ) ( ) ( ) , 1 ( / ) ( ) ( ) M V k uk M M V k m k    (4)
where

2 B 2 2 0 ˆ( ) / , ( ) ( , ) , k r k m k m         kk kk (5)
(k) is the Fourier transform of the Coulomb interaction, c=N/M is the lattice concentration of particles, V is the system volume, m 2 (k) is the second cumulant of the density fluctuations, 0 designates the averaging over the reference distribution. In [START_REF] Yukhnovskii | Statistical Theory of Classical Equilibrium Systems[END_REF] it was the averaging over the system of hard spheres; in our case it is the averaging over the ideal crystal host system. The lattice version of the collective variables [41] generalized to take into account particle displacements from the lattice sites reads as

,0 0 1 1 ˆˆêxp( ) , , ,
M i i i i i i i n i M c c n M           kk k q q R r (6)
where ˆi n is the occupation number operator ( ˆi n =1 if the lattice cell v i is occupied by a particle, ˆi n =0 if the lattice cell is vacant; multiple occupation of a cell is forbidden), c is the mean lattice concentration of particles, ,0  k is the Kronecker's -symbol, R i and r i are the radius-vectors of the lattice sites and the particle position with respect to its lattice site, correspondingly.

For calculating the second cumulant it is useful to note that

1 0 ˆexp( ) exp( ) ( ), ( ) exp( ) ( ) d , , i i i i i i i i i v n i c i f f i F v        k q k R k k k r r r r (7)
where normalized to 1 unary distribution function is used

3 1 ( ) d 1, / , i i i i v F v v V M     rr (8)
v i =v is the volume of a lattice cell. The Fourier transform of the concentration fluctuation correlation

1 2 1 2 2 1 2 ,0 ,0 0 11 0 1 exp[ ( )] . MM ij ij i Mc M           k k k k k q + k q (9)
The facts that f(0)=1 and ,0 1 exp( ) (10) were used deriving Eq. ( 9). The first term in the r.h.s. of Eq. ( 9) can be split into two sums resulting in

M i i iM      k kR
12 11 0 1 exp[ ( )] MM ij ij i M      k q + k q 1 2 1 2 1 1 1 0 0 11 exp[ ( )] exp[ ( ) ] . M M M i j i i j i ji ii MM            k q + k q k + k q (11)
The first term in the last expression

2 1 2 1 2 1 2 1 1 1 1 0 11 exp[ ( )] ( ) ( ) exp[ ( )] M M M M i j i j i j i j j i j i i c f f i MM              k q + k q k k k R + k R 1 2 1 2 2 1 2 ,0 ,0 ,0 ( ) ( ) . c f f M         k k k k kk (12)
The first term in the r.h.s. of Eq. ( 12) cancels with the last term of Eq. ( 9) and the final result is

1 2 1 2 1 2 ,0 0 (1 ) ( ) ( ) , c c f f       k k k k kk (13)
and the second cumulant 2 2

(1 ) at 0, ()

(1 ) ( ) at 0.

c c k mk c c f k k       (14)
In the approximation, where the unary distribution function is a spherically symmetric function we find that

( ) ( ) ( ) 1 f f f k Jk Wk        kk . ( 15
)
where the second and fourth moments of the distribution function are

2 3 1 1 [( / ) ] ( ) , 2 l v J k F r d   k r r ( 16
) 3 1 1 [( / ) ] ( ) . 4! l v W k F r d   k r r (17)
Substituting Eqs. ( 14) and ( 15) at k0 into Eq. ( 4) for the Fourier transform of the screened Coulomb potential one arrives at the expression

2 2 4 1 ( ) ( ) , 1 ( / ) ( )[(1 ) (2 ( ) ) ] k uk V N V k c c Jk J W k          (18)
that means the Debye length in the k0 limit in the lattice version includes in the denominator an additional multiplier (1-c)

3 0 23 ,, (1 )
D h N Mc c r VV c c e b      ( 19
)
where b is the lattice spacing in the case of a simple cubic lattice.

In accordance with Eq. ( 19) the Debye screening length is symmetric with respect to c and (1-c) as it should be because the description of solid electrolytes can be performed for real ions of concentration c or the Kröger-Vink notations [START_REF] Kröger | The Chemistry of Imperfect Crystals[END_REF][START_REF] Kröger | Relations between the concentrations of imperfections in crystalline solids[END_REF] can be used, when lattice defects (vacancies) of concentration 1-c are considered. The physical quantity, the Debye screening lengths must not depend on the description used. The Debye screening length tends to infinity in both limits c0 and c1. The concept of vacancies was used to explain the increase of the screening length in concentrated liquid electrolytes when solvent molecules are considered as specific defects that make inaccessible to ions corresponding regions of the system [START_REF] Lee | [END_REF]. This is an example of transfer of solid state representations into understanding the behavior of liquids that were originally introduced in [13,[START_REF] Frenkel | Kinetic theory of liquids[END_REF].

It is worth to note that given in Eq. ( 19) the nonmonotonic concentration dependence of the screening length is conditioned by the symmetry properties of the ideal entropy contribution in the mean field approximation. This is conceptually different from such an effect in concentrated electrolytes and ionic liquids where it is concerned with creation of cation/anion complexes, or due to accounting of interparticle correlations [44][45][46].

In a more general case as it follows from Eq. ( 18) the Fourier transform of the screened Coulomb potential in the denominator includes the wave vector dependent terms that account for the displacement of particles from the lattice sites.

Thermodynamics of the system

For the configuration integral of the system the identity

C0 0 0 exp( ( ) , N N N N N Q Q U U U       ( 20
)
can be used where

(| |) N N i j ij U       qq is the energy of short range interactions, C1 B N N ij ij U r r    
is the Coulomb interaction energy; 0 1 ()

MM N j i i j i U     q
and Q 0N are the energy and configuration integral of the reference system, respectively, that in our case is taken as an ideal crystal [START_REF] Rott | Statistical Theory of Molecular Systems[END_REF][START_REF] Rott | [END_REF] and thus

1 0 0 , Q exp[ ( )] , i M M c N j i i ji v Q Q c d                    qq (21) where () ji 
 q is the mean potential exerted by a particle in the lattice cell v j on a particle in the lattice cell v i . The leading term in the expression for the binary distribution function is calculated in accordance with Eqs. ( 1) and ( 2). The averaging over reference distribution results in the equation

  C 0 C 0 0 exp( ( ) 1 exp( ) , M N N N ij N ij U U U f U           (22) exp{ [ (| |) ( ) ( )]}, ij j i j i f         q q q q (23)
where f ij are the renormalized Mayer functions. The mean potentials are calculated at the condition that the two-vertex diagrams in the expansion of the configuration integral are equal to zero [47]. This condition leads to results that are equivalent to the quasichemical (or Bethe-Peierls) approximation.

Averaging over the reference distribution at the mentioned condition leads to the expression

C 0 C 00 exp( ( ) exp( ) N N N N U U U U        . ( 24
)
Using the notation C 0 ln exp( )

cN AU    ( 25
)
and Eq. ( 20) the free energy A of the system per lattice site can be represented by the expression

1 0 (ln ln ) , c A c c Q A          (26)
where the sum represents the ideal entropy contribution and the lattice part of the free energy due to the short range interactions. The last term accounts for the correlated part of the Coulomb interaction and in accordance with [START_REF] Yukhnovskii | Statistical Theory of Classical Equilibrium Systems[END_REF] can be written as

0 , cc A u d    (27)
where c u is the density of the mean Coulomb interaction energy. For a homogeneous bulk, due to local electroneutrality conditions, it is determined by the correlation part of the electric potential

1 2 B 1 1 1 1 1 1 ( ) ( ) ( , ) d d . j M c j j j j vv r u c F F h r      q q
q q q q (28)

The correlation function can be taken in the form that was used for developing the mean spherical approximation [48,49] ( ), () 1, u r r d hr rd

      ( 29
)
where d is the hard core diameter of the particles.

The sharp maxima of the one-particle distribution functions can be approximated by the Gaussian distribution

3/2 2 2 1 ( ) ( / ) exp( ), 1, . F q d d b        q (30)
Then the correlated part of the mean Coulomb interaction is represented as follows

2 1 2 , M cj j u c u    (31) 3/2 1 2 B 1 1 () exp( ) , j j jj j v hr r u q dv r         (32)
where 1 j r is the distance of j-th particle from a chosen lattice site, q j is the displacement of j-th particle from its lattice site.

The integral in Eq. ( 32) can be calculated in the bipolar coordinate system where R 1j is the distance between the first and j-th lattice sites, r 0 =b/2 is the radius of a lattice cell. At the condition 2 1 d   , the shape of the lattice cell has no importance because the main contribution to the integrals in Eqs. [START_REF] Kröger | The Chemistry of Imperfect Crystals[END_REF], (32) and (34) comes from the regions close to the lattice sites.

Taking into account that h(r) is a smooth function of r and is multiplied by the function with a strong maximum at the center of the lattice cell the integral in Eq. ( 34) can be evaluated by the Laplace's method thus giving

1 1 1 () . j j j hR I R  ( 35
)
Passing to a description in terms of a continuous medium and performing spatial integration in accordance with Eq. ( 29) for h(r) we arrive at the following expression for the correlated part of the Coulomb interaction energy in the mean spherical approximation

Finally, the Coulomb interaction contribution to the free energy follows after integration over temperature

indicates (Fig. 1) that the homogeneous distribution of mobile charges becomes unstable already at rather low concentration of particles where the concentration derivative of the chemical potential becomes negative. The region of the system stability can be clarified due to more consistent accounting of the short-range interactions and correlations among the mobile ions and especially with the particles of the neutralizing background. The role of charged defect interactions in phase transitions of ionic crystals was discussed in [50,51].

Charge distribution at a plain boundary

In the first approximation the statistical sum Q M of a subsystem of mobile (e.g. positive) ions on the lattice of M sites can be calculated from the expression

  B 1 ln ln (1 ) ln(1 ) 2 MM i j i j M i i i i i i j ij ij c c c c r Q c c c c rr                 e 11 , 2 2 2 M M M M M ij i j ij i j ij ij i j i i i i i j i j i j i i h c c J c c J h c c e c u c                      (41)
where c i and j c are the lattice concentrations of ions and counterions, respectively, and ij ij rr  are the distances between the corresponding ions, J ij the short range interaction potential between ions, e i  the external potential, u i the crystal potential distribution, i.e. the potential energy well depths that are created by the host system for mobile particles.

In Eq. ( 41) the interparticle correlations are taken into account through lattice version of the binary correlation function h ij , which is defined through the expression , 1,

ij i j ij ij ij F c c g h g    (42)
and its modification for Coulomb interactions

/. ij ij ij h h r  (43)
Due to the condition for the weakening of correlations the terms of Eq. ( 41) containing functions h ij and ij h describe the short range contributions. The chemical potential follows from differentiation over the lattice concentration B 1 ln ln 1

M jj i M i j i ij ij i ji cc c Q r c r r c                    e 1 1 1
.

M M M ij j ij j ij ij j i i j j j j i j i j i h c J c J h c e u                    (44)
It is assumed here that the concentration dependence of the binary distribution functions F ij is mainly represented by the concentration multipliers in Eq. ( 42) and the concentration dependence of the correlation functions can be neglected.

In the bulk of the system at zero electric field

B 1 1 1 1 ln . 1 M M M M j ij ij ij ij j j j j ij ij j i j i j i j i c cc r c h c J c h J u c r r                                (45)
The equilibrium condition  i = results in the expression for the deviation c i =c i -c of the ion concentration from its bulk value (the concentration of counterions is considered homogeneous through the whole system)

1 ln 0, 1 i i i i i cc c w e u c c c              ( 46
)
where the total electric field e B 1

M j ii j ij ji c r r         . ( 47
)
and the short range interparticle interaction contribution 1 ( ) , .

M i ij ij ij j i i j ji w g J h c u u u           ( 48
)
The short range interactions account for the interparticle correlations that are similar to that considered in a recent paper [52] where the lattice gas model of liquid electrolytes was used.

The variation of the potential energy well depths u i can appear near the system boundary and result in special surface states for ions.

The solution of Eq. ( 46)

1 exp[ ( )] 1 [ / (1 )]exp[ ( )] i i i i i i i c u e w c c c u e w                  ( 49
)
can be rewritten for the concentration distribution

1 , 1 [(1 ) / ]exp[ ( )] i i i i c c c u e w          ( 50
)
that looks like a Fermi-Dirac distribution [24][25][26][START_REF] Dutta | Proc. Nat. Inst. Sci. India[END_REF][START_REF] Kröger | The Chemistry of Imperfect Crystals[END_REF][START_REF] Yu | [END_REF][30][31][32][33][34] corrected for inclusion of the interparticle interaction contribution and concentration dependent multiplier in front of the exponent. This is a consequence of the occupation number definition to take values 0 or 1.

In fact, Eq. ( 50) is not an explicit solution for the concentration distribution because it contains the concentration dependence in the r.h.s. through the correlation term w i representing the difference analog of an integral equation. This is a fundamental difference of Eq. ( 50) from the conventional Fermi-Dirac distribution derived for independent particles or in the mean field approximation. It can be solved iteratively. The approach suggested in Ref. [52] to solve such type of nonlinear equations through their linearization with subsequent renormalization of the control parameters can be used.

In the limit of small concentration this distribution approaches to the Maxwell-Boltzmann distribution, again corrected for including interparticle correlations exp[ ( )].

i i i i с c u e w        (51)
In the Debye-Hückel approximation 0, 0, 1

i i i u w e        (52)
Eq. ( 49) reduces to (55) that means the Debye length is again given by Eq. ( 19) and the width of the double layer diverges at low (c0) and high (c1) ion concentrations.

Similar expression for the Debye length was deduced [30] on the basis of thermodynamic arguments taking into account that the concentration derivative of the chemical potential has to be proportional to the product c(1-c) [more precisely to c(1-c max )] that was explained by the symmetry of low particle and low lattice defects concentrations.

In the general case the equation for the electric potential (56) can be rewritten for the electric field E=-d/dz, axis z is perpendicular to the system boundary

3 0 1 exp[ ( )] (1 ) 1 exp[ ( )] i i i i i i i u e w c c e c c u e w h                   
0 22 0 3 0 (1 ) 1 exp[ ( )] 2, 1 exp[ ( )] c c e u e w E E d c c u e w h                        (57)
where E 0 and  0 are the electric field and potential at the system boundary.

The double layer differential capacitance can be evaluated from the expression [30,33]

2 0 0 0 0 0 0 0 . 2 dE dE C d E d      (58)
At the conditions =0 and E=0 in the bulk it follows from Eq. (57) that

2 0 0 0 0 3 0 0 0 0 0 1 exp[ ( )] (1 ) 2, 1 exp[ ( )] dE u e w c c e d c c u e w h                    (59) 0 1/2 00 3 0 0 (1 ) 1 exp[ ( )] 2 sign . 1 exp[ ( )] c c e u e w Ed c c u e w h                           (60)
At u=0 and w=0 Eqs. ( 50), (60) coincide with those of Ref. [30]. As is mentioned above the well depths u i can vary in the region close to the system boundary that will result in additional ion redistribution close to the boundary. Positive u i values correspond to the enrichment of the surface layers by moving ions.

For qualitative estimation of the influence of this surface potential on the double layer characteristics we adopt that u differs from zero in the first surface layer of thickness b only and consider the mean field approximation when the binary distribution function is represented by the product of unary ones and then w=0.

The integral in Eq. ( 60) can be split into two parts from 0 to  1 and from  1 to  0 , where  1 is the value of  at z=b/2, the position of ion centers in the first surface layer. u contributes in the latter integral and the result of integration can be represented as follows is taken into account. The nonlinear system of Eqs. ( 61), (63) can be solved by the method of successive approximations taking e.g.  1 =0.8 0 as an initial trail value. Fig. 2 demonstrates that the solution of Eqs. ( 61), ( 63) for E 0 is a smooth function of the electric potential. Thus, the method of successive approximations provides with the required solution. The differences between the solutions for B1 and B=1 are not very pronounced although they increase with increasing the difference of B from one. At  0 close to zero the values of E 0 are close to zero as well. In a narrow region of surface potential values close to zero for B>1 and  0 <0 or B<1 and  0 >0 the solution for E 0 attains imaginary values that means there is no charge distribution in the electrolyte that satisfies Eqs. ( 49), (61) with the value of  1 found selfconsistently from these equations. The width of this region increases with increasing the difference of B from one. The differential electric capacity as a function of the surface electric potential at particular values of the parameters for different values of the variation of the crystal surface potential is shown in Fig. 3. The surface crystal potential variation leads to strong deviation of the differential capacitance at low absolute values of the surface potential and for the negative potentials the difference is more pronounced (Fig. 4). The capacitance diverges at values of  0 at which the electric field tends to zero because at u0 the derivative of square of the electric field over the electric potential is not equal to zero. At the same time the derivative of square of the electric field over the electric potential changes its sign at e 0 =u (see Fig. 4; the inset displaces the behavior of this derivative at small values of the surface potential) that results in negative values of the capacitance in the corresponding ranges of the electric potential and can result in thermodynamic instability of the system. This means that there are no solutions for the charge distribution in the system that corresponds to the surface potential values at which the capacitance is negative.

  1 2 00 2 D 21 ln 1 e () e E e c c c er              0 1 3 0 0 2 1 1 c e ln , exp (
From physical point of view, the negative values of the capacitance and the system instability are due to the competition of the external electric field and internal crystalline field. It should be noted that according to Eq. ( 48) the short range interparticle interaction contribution w i depends on the concentration distribution and Eqs. ( 49) and ( 50) are nonlinear equations with respect to the concentration that have to be solved selfconsistently. The account of the influence of the charge distribution in the near surface region through the short range interparticle interaction contribution can significantly change the behavior of the system at low absolute values of the surface potential. On the other hand, the possibility of a negative differential capacitance of the double layer and its divergences were discussed in numerous papers (see [53][54][55][56] and references therein) and the problem of system stability in the case of the model under consideration requires additional consideration.

Conclusion

The lattice fluid model of the system with short range and long range Coulomb interactions is suggested. In the framework of the collective variables method, the screening of the Coulomb interactions in the bulk is considered and it is shown that the Debye length includes additional concentration dependence inversely proportional to the square root of the mean concentration of vacant sites. Additionally, the contribution of the thermal ion displacements from the lattice sites to the Fourier transform of the screened Coulomb interaction is estimated. This contribution results in a more complicated wave vector dependence of the screened Coulomb interaction.

The Coulomb interaction contribution to the free energy of the system is calculated in the approach close to the mean spherical approximation and is given in an analytical form. The analysis of the concentration dependence of the chemical potential shows that this contribution can lead to thermodynamic instability of the system already at rather low mobile ion concentrations. The direct account of the interaction with the immobile neutralizing impurities and short range interion correlations can significantly enlarge the stability region of the system.

In the mean field approximation, the mobile ions near the system boundary obey the Fermi-Dirac type distribution. When the interparticle correlations are taken into account, this distribution reminds the difference analog of the integral equation for the concentration distribution. The influence of the variation of the crystal field near the system boundary on the structure and characteristics of the electric double layer is investigated. As compared to the system with equal crystal potentials at the lattice sites throughout the system, the pronounced difference for the electric capacitance appears at low absolute values of the surface potential. For positively charged particles this difference is more pronounced for negative electric potentials.

The capacitance diverges as the potential values at which the electric field tends to zero and attains negative values in the regions of the surface potentials depending on their polarity and values of the surface crystal potential. From mathematical point of view, the capacitance diverges at values of the surface potential at which the electric field tends to zero while the derivative of the square of the electric field over the electric potential is not equal to zero. From physical point of view, the negative values of the capacitance and the system instability are due to the competition of the external electric field and internal crystalline field. Negative values of the capacitance may indicate the thermodynamic instability of the system although more consistent accounting of the short range interaction contribution can considerably change the behavior of the system at low absolute values of the surface potential.

24

 24

  three terms of its expansion in the first surface layer can be used

Fig. 1 .

 1 Fig. 1. The chemical potential versus concentration at =40, T=1000 K, b=0.4 nm, d=0.35 nm, r D =1.1 nm, r B =0.42 nm.

Fig. 2 .

 2 Fig. 2. The electric field in units

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. The differential electric capacity versus electric potential for positively charged particles. C 0 =/r D , B=exp(u), =40, T=1000 K, b=0.4 nm, c=0.01, r D =1.1 nm.

Fig. 1 .

 1 Fig. 1. The chemical potential versus concentration at =40, T=1000 K, b=0.4 nm, d=0.35 nm, r D =1.1 nm, r B =0.42 nm.

Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .

 2345 Fig. 2. The electric field in units

Acknowledgement

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skodowska-Curie grant agreement No 734276, the Belarusian Republican Foundation for Fundamental Research (grant №Ф16К-061), the State Fund for Fundamental Research of Ukraine (grant №Ф73/26-2017), and the Ministry of Education of Belarus.


Samuel Hidalgo

Alvaro Cassinelli

Emmanuel Fort

Matthieu Labousse
email: matthieu.labousse@espci.fr

Frugal random exploration strategy for shape recognition using statistical geometry

Very distinct strategies can be deployed to recognize and characterize an unknown environment or a shape. A recent and promising approach, especially in robotics, is to reduce the complexity of the exploratory units to a minimum. Here, we show that this frugal strategy can be taken to the extreme by exploiting the power of statistical geometry and introducing new invariant features.

We show that an elementary robot devoid of any orientation or observation system, exploring randomly, can access global information about an environment such as the values of the explored area and perimeter. The explored shapes are of arbitrary geometry and may even non-connected.

From a dictionary, this most simple robot can thus identify various shapes such as famous monuments and even read a text.

One-Sentence Summary:

Blind robots, with no compass and moving randomly, can read and recognize complex shapes based on a statistical approach.

Main Text:

Image, object or environment recognition is a highly complex operation which involves acquiring a global information based on spatial correlations. The standard strategy is centralized using multiplexed sensors coupled to complex computational capacities and data processing which can be possibly embarked in high-tech robotic systems. Recently, an alternative fruitful strategy has emerged based on the simple cost-effective basic units collecting unsupervised information, no direct localization capacities and eventually working in parallel in an unmanned manner [START_REF] Rubenstein | Programmable self-assembly in a thousandrobot swarm[END_REF] [START_REF] Ben Zion | Morphological computation and decentralized learning in a swarm of sterically interacting robots[END_REF]. This strategy is more robust and well-suited for analysis and exploration with restricted computational and hardware resources, and in the case of harsh or complex environments for which part of the information is not accessible. The most elementary system one can imagine is limited to the acquisition of a local information only with no possible communication, no localization capacity or orientation awareness.

The unavailability of localization capability and the associated randomness of the exploration make the use of statistical geometry crucial to move from local information to a global knowledge of the environment. In this context, the statistical properties of chord length distributions across various kinds of geometrical shapes play a central role for the characterization of size and shape of the intercepted objects with applications such as in acoustics, ecology [START_REF] Mcintyre | Estimation of plant density using line transects[END_REF]- [START_REF] Frangipane | Invariance properties of bacterial random walks in complex structures[END_REF] image analysis [START_REF] Matheron | Image analysis and mathematical morphology[END_REF], [START_REF] Mazzolo | Probability density distribution of random line segments inside a convex body: Application to random media[END_REF], stereology [START_REF] Underwood | Quantitative stereology[END_REF]- [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF] and reactor design [START_REF] Dirac | Approximate rate of neutron multiplication for a solid of arbitrary shape and uniform density[END_REF]- [START_REF] Mazzolo | On the properties of the chord length distribution, from integral geometry to reactor physics[END_REF]. The moments of the chord length distributions are related to mathematical invariants and give access to global geometrical parameters such as the volume, the surface or the perimeter of objects. As an historical example, Cauchy formula states that the mean chord length is proportional to the ratio of the volume of an object to its surface area for a 3D object, or to the ratio of its area to its perimeter for a 2D object [START_REF] Cauchy | Mémoire sur la rectification des courbes et la quadrature des surfaces courbes[END_REF], [START_REF] Santaló | Integral geometry and geometric probability[END_REF]. Recent generalization to random motion [START_REF] Blanco | An invariance property of diffusive random walks[END_REF]- [START_REF] Shukla | Surprising variants of Cauchy's formula for mean chord length[END_REF], to closed trajectories [START_REF] Mazzolo | Invariance properties of random curves: an approach based on integral geometry[END_REF], [START_REF] Hidalgo-Caballero | Mean arc theorem for exploring domains with randomly distributed arbitrary closed trajectories[END_REF], and to wave propagation properties [START_REF] Bardsley | The average transport path length in scattering media[END_REF]- [START_REF] Savo | Observation of mean path length invariance in light-scattering media[END_REF] have renewed interest for such invariants. Unfortunately, most of theorems of this field of mathematics are restricted to convex shapes and with no inclusions which restrict the statistical approach for exploration and pattern recognition.

Here, we introduce statistical invariants that extend to any arbitrary, possibly non-convex, nonconnected or with inclusions shapes. This generalization allows us to revisit and push the limits of frugal exploration strategies to its most minimalistic version by fully integrating the assets of statistical geometry. We show experimentally that an elementary robot devoid of any localization, orientation and observation systems can assess the surface and perimeter of a region in which it evolves with random displacements. Then, with the help of a dictionary, we show how it recognizes these shapes. We also discuss the benefits of global over sequential exploration strategies in the case of non-connected shapes and apply it to reading strategies.

Generalized statistical geometry invariants

Statistical geometry establishes a link between probabilities arising from local properties and global geometric quantities. Cauchy's theorem plays a central role in relating the average chord length of a random distribution 〈ℓ〉 chords to the ratio of the area 𝐴 to the perimeter 𝑃 of the explored domain. The average is realized over all possible chords of the domain. In 2D, it reads 〈ℓ〉 chords = 𝜋𝐴/𝑃. The Crofton-Hostinský theorem uses the third moment, and in 2D satisfies 〈ℓ 3 〉 chords = 3𝐴 2 /𝑃 so that combining these two theorems can yield 𝑃 and 𝐴 separately. While Cauchy's theorem can be extended to non-convex shapes [START_REF] Mazzolo | Properties of chord length distributions of nonconvex bodies[END_REF] unfortunately, the Crofton-Hostinský theorem is only valid for convex shapes, which in practice, prohibits its relevance for exploration and pattern recognition purposes. A generalization of this theorem to arbitrary forms is therefore essential. To this end, we introduce a new statistical invariant. We first label the intersections between a straight trajectory 𝐿 and the boundaries of the explored shape according to whether they are ingoing 𝑖 or outgoing 𝑜 from the shape. We then derive all the distances between pairs of intersections by distinguishing the chords depending on their ends: 𝑙 𝑖𝑜 , 𝑙 𝑖𝑖 or 𝑙 𝑜𝑜 . We define the geometric function of order 𝑛 ∈ ℕ associated to the line 𝐿, ℒ 𝑛 (𝐿

) = ∑ ℓ 𝑖𝑜 𝑛 𝑖,𝑜 -∑ ℓ 𝑜𝑜 𝑛 𝑜 -∑ ℓ 𝑖𝑖 𝑛 𝑖 .
We find that the area of the explored shape 𝐴 satisfies

𝐴 = 3 𝜋 〈ℒ 3 〉 lines 〈ℒ 1 〉 lines (1)
with the average being made with respect to all straight lines crossing the shapes (see Supplementary Materials 1 for the demonstration and graphical illustration of Eq. 1). This relation, combined with Cauchy's theorem, allows statistical geometry to be used for exploration and pattern recognition by measuring 𝑃 and 𝐴 independently. It applies to any arbitrary shape, whether nonconvex, with inclusions, and even non-connected. This establishes an autonomous strategy for obtaining global geometric features without any a priori assumptions.

Implementation with an elementary robot.

We experimentally implement this frugal exploration strategy to evaluate the robustness of Eq. 1 under real conditions. We use a commercial toy-robot called Sphero which is a self-propelled rolling sphere (see Fig. 1A and Supplementary Materials 2). It possesses an internal clock, two encoded motors for motion and a light sensor. The shape to recognize such as the Eiffel Tower (Fig. 1B), is projected using a video-projector on the arena explored by the robot. The robot performs a series of uncorrelated straight linear motion at 1.0 cm.s -1 , changing direction randomly as it reaches the edge of the arena (Fig. 1B, see Supplementary Movie 1). Crucially, the robot possesses no information about neither its absolute position nor its direction. From the light received by the sensor, the robot can only detect if it is inside or outside the projected shape. From the times of entries and exits of the shape along the line 𝐿, the chord lengths are deduced and the functions ℒ 1 (𝐿) and ℒ 3 (𝐿) are calculated on board. From the statistical averages, the area 𝐴 and the perimeter 𝑃 of the projected shape are deduced using Eq. 1 and Cauchy's theorem. Note that the memory required to perform the calculations is minimal, containing only the time of intersection events along the latest line to update the previously calculated mean values. Shape recognition is then performed using a dictionary containing a list of area and perimeter entries for a given shape, associated in the present case to famous monuments (Figs 1B-E, see Supplementary Materials 2). Although two different objects can have the same entries, we found in practice that this was hardly ever the case, with up to hundreds of entries. Figure 1B shows the accumulated experimental trajectories during the robot exploring the Eiffel Tower for 𝑁 = 10, 100 and 1000 random lines. An external camera was used to follow the on-board led which blinks when it is inside of the domain. The final exploration for the other shapes is shown in Figs 1C-E are geometrical invariant (while 〈ℒ 1 〉 lines and 〈ℒ 3 〉 lines are not). The trajectories followed by the estimated parameters and areas as the number of lines 𝑁 increases (increasing line contrast) can be plotted in a Perimeter-Area representation space for various explored shapes (Fig. 1G, see Supplementary Movie 2). The statistical confidence can be estimated by Monte Carlo simulations for 𝑁 = 1000. The increasingly small ellipses in Fig. 1G show the 75%, 95% and 99% confidence boundaries associated to each shape. This proves that the robots discriminate each registered shape in the dictionary without any statistical ambiguity even though some shapes are very similar, such as in the case of the triangle and the Chichen Itzá pyramid (chosen with an identical area)

The extreme frugality and efficiency of the recognition pathway relies on iterative averaging of the statistical features. Figure 2A shows the distributions of the chord lengths normalized by the maximal length associated to various explored shapes. The full chord length distribution obviously holds some additional geometric information from the shape as shown by distinct profiles. However, in the quest for the simplest exploration strategy, the robot does not need to store all these the data. It only iteratively computes and save the mean values 〈𝑙〉 chords , 〈ℒ 1 〉 lines and 〈ℒ 3 〉 lines . This is performed after each exploration line. Note that even the full knowledge of the chord PDF would not be sufficient to establish a strict correspondence with all possible 2D shapes [START_REF] Mallows | Linear-intercept distributions do not characterize plane sets[END_REF]. For a non-convex shape, a line contains at most 𝑘 ingoing and outgoing intersection points (𝑘 = 1 for convex shapes). The robot thus needs to temporarily store at most 𝑘(𝑘 -1) +4 values.

For a complex shape such as the statue of liberty (Fig. 1D), 𝑛 ≤ 6 and the robot needs to store at most 34-scalar values. While memory is not restrictive for large robots, it is essential applications at small scales.

We now investigate the convergence to assess the efficiency of the probabilistic method. Figures 2B and2C show the evolution of the surface area 𝐴 and the perimeter 𝑃 respectively as a function of the number of lines 𝑁. The values are computed by the robot during experiments associated to the five shapes in the shape dictionary (colored contrasted solid lines). The robot is programmed to send these updated geometrical estimates by wireless Bluetooth connection only for visualization purposes as it performs an on-board surface and perimeter estimations. Additionally, we perform Monte Carlo simulations by generating uniformly distributed random lines (colored dimmed solid lines) to compare with the experiments. The convergence towards the expected values of 𝑃 and 𝐴 for the experimental and simulated estimates are similar. It is found that only a few hundreds of lines are typically required to measure the area and perimeter with an error of less than 10%. This is sufficient to distinguish the shapes listed in our dictionary. However, in the general case, the recognition speed would strongly depend on the L2-norm in the representation space {Perimeter-Area} between shapes listed in the dictionary. The standard deviation over many independent realizations can be computed from the numerical simulations as a function of exploration lines 𝑁. Figures 2D and2E show the evolution of the standard deviations of 𝐴 , 𝜎 𝐴 (𝑁), and of 𝑃, 𝜎 𝑃 (𝑁), respectively for the five shapes of the dictionary. It is found that the convergence, in both cases, is independent of the complexity of the explored shape and scales as 1/√𝑁. This can be rationalized by noting that the measure converges towards a normal distribution as an ensemble averaging of independent random variables with the same mean and standard deviation.

The influence of the shape appears only in the prefactors 𝜎 0,𝐴 and 𝜎 0,𝑃 associated to 𝐴 and 𝑃, respectively. The latter increases from the circle, to the square, to the triangle, to the Chichen Itzá pyramid, to the Eiffel Tower, and finally, to the Statue of Liberty. The prefactors thus seem to increase with the complexity of the shape. For regular shapes, complexity can be related to the decrease of the symmetries. For irregular shapes, complexity needs to be assessed quantitatively.

The Kullback-Leibler divergence 𝐷 KL which is a measure of the relative entropy between two probability distributions provides such a parameter. Figure 2F shows the prefactors 𝜎 0,𝐴 and 𝜎 0,𝑃 as a function of the 𝐷 KL of the chord length distributions associated to the various shapes, with the chord length distributions of a circle taken as a reference. The prefactors increase monotonically with 𝐷 KL confirming the role of complexity. However, the variation is only a factor of about 4 for a logarithmic variation between distributions of about 3. The dependence on shape complexity is therefore rather modest.

In practice, the exploration must be stopped, once the likelihood of determining a shape is larger than a preset value here 95%. The corresponding likelihood landscape in the representation space {Perimeter-Area} depends on the given dictionary, 𝜎 0,𝐴 and 𝜎 0,𝑃 all of the entries, their relative distance in the {Perimeter-Area} space and the duration of the exploration. Figures 2G and2I illustrate the influence of the number of crossing lines, N, where each colored domain delineates a 95% likelihood of finding the targeted shape (for N=200 and 1000, see Supplementary Materials 3). Figure 2H and 2J show quantitatively the relative number of lines contributing to a 95% likelihood of recognition. One observes that about a hundred lines are enough to recognize a circular shape with this level of certainty which is due to its isolation from the other shapes in the representation space {Perimeter-Area}. In contrast, it takes more than 10 3 lines to achieve this level of confidence for the other shapes that are more closely packed in the representation space {Perimeter-Area}.

Eq. 1 remains valid for arbitrary shapes with inclusion or even non-connected shapes.

Disconnected objects can thus be explored and gathered within a common global shape with area or perimeter equal to the sum of those of its parts. Since the area and perimeter are geometrical invariants, they are conserved under rigid body transformation such as translations, rotations, mirrors, and any combinations of them (see Supplementary Materials 4). Two alternative strategies can be chosen to recognize unconnected shapes: a global strategy considering all the shapes as a single clustered entity or a local strategy considering the addition of each sub-entity separately.

We assess these two strategies through the emblematic case of reading, based on a word or letterby-letter recognition strategy. The local strategy consists of mapping the capital alphabetic letters into the {Perimeter-Area} representation space which depends on the font. Figure 3A shows the position of the projected capital letters for a Helvetica font and the associated domain where the likelihood of recognition is larger than 95% for N=3 × 10 4 . Most common fonts do not overlap in the dictionary partly because letters do not transform into each other by symmetry. Figure 3B specifies the number of lines required to distinguish a given letter with a likelihood larger than 95%. Letters can thus be distinguished without ambiguity. For the global strategy, the dictionary is composed of words in the {Perimeter-Area} representation space. Figure 3C shows the dictionary with words sampled from the preamble of the United Nations charter as a dictionary [START_REF] Assembly | Universal declaration of human rights[END_REF] with the associated 95% likelihood domains. All the words can be distinguished. One necessary condition for non-overlapping in the dictionary is the absence of anagrams. In order to compare with the local strategy (Figure 3B), Figure 3D specifies the number of lines, normalized by the of letter in the word, required to distinguish a given letter with a probability larger than 95%. For this moderately-small pool of words, the two strategies are of about the same efficiency.

Figures 3E and3F illustrate the principle of these strategies applied to the word "FREEDOM" with successive identification of each letter (Fig. 3E) and the global estimate (Fig. 3F) of the word perimeter and area. We measure in Fig. 3G and 3H the perimeter and the area of the word FREEDOM. With this moderately-small pool of words, the two strategies are of about the same precision.

The best strategy is a compromise between the size, the choice of the word dictionary and the accumulation of independent errors in the letter-by-letter strategy. On one hand, if the word dictionary is larger than our current one, the words accumulate in a closed region of the Perimeter-Surface space, and as a consequence, a letter-by-letter strategy is more efficient. On the other hand, if we choose a letter-by-letter strategy, the perimeter and surface area computation accumulate the uncertainties corresponding to each letter. In contrast, the direct word recognition combines the chords topology such that the overall noise is smaller. So, for smaller pool of words, measuring directly the word FREEDOM converged quickly by a "word" rather than a letter-by-letter recognition.

The strategy, present here, does not require any knowledge about neither absolute nor relative position as well no information about the direction of motion. On the contrary, the randomness of orientation and absolute position is crucial. At every single stage of the physical computing, the robots acquire no global information about the number of holes, the convexity or any global information. However, by being ignorant about any global information and by measuring only local quantities, the robot is able to recognize shapes listed in a dictionary by measuring the perimeter and the area of any given 2D shapes. This was only possible thanks to the introduction of a new generalized statistical invariant which allows us to estimate the area from a random exploration of any 2D domain, regardless of its complexity. Finally, we note that the dictionary is not a necessary prerequisite and is used here for convenience. A pre-classification phase, for example with k-means or Gaussian mixture strategies, can be added to our protocol to perform an autonomous and unsupervised recognition. Extending our technique to a swarm of robots, to take advantage of distributive computation strategies, is a natural step where each can compute these invariants locally and share them with the others to improve convergence. However, the inevitable contact between individuals and the possible formation of aggregates makes swarm computation a real challenge, for which new invariants for more disordered paths need to be discovered.

  (𝑁 = 1000). The explorations are found homogenous and statistically isotropic. The chord distribution for the Eiffel Tower is shown in Fig 1F, its average value 〈ℓ〉 chords and 3 𝜋 〈ℒ 3 〉 lines 〈ℒ 1 〉 lines

Figure 1 Fig. 1 .

 11 Figure 1

Figure 2 Fig. 2 .

 22 Figure 2

Figure 3 Fig. 3

 33 Figure 3

Funding: European Union's Horizon 2020 (Marie Skłodowska-Curie grant No 754387) AXA Research Fund CONACYT Mexico (mobility Grant CVU-849369) French National Research Agency (ANR-10-LABX-24).

Supplementary Materials

Supplementary Text 1-5

Figs. S1 to S2

Tables S1

References ( [START_REF] Mazzolo | On the properties of the chord length distribution, from integral geometry to reactor physics[END_REF]) Movies SM1 to SM2