N Adnet

I Bruant

F Pablo

L Proslier

A mixed hexahedral finite element to model the scattering of a deformed microstrip antenna

Keywords: Boundary integral method, electromagnetic/mechanical coupling, finite element method, Maxwell's equations, microstrip antenna, scattering, vector finite element

Microstrip antennas have a major interest in aeronautical applications due to their low profile. This paper deals with the impact of mechanical strain on the scattering properties of these antennas. Considering a weak coupling between electromagnetism and mechanical behavior, the same 3D hexahedral finite element discretization is used to solve both problems. A node-based approximation is used for mechanical displacement, while for the determination of the electromagnetic fields, a vector finite approximation is implemented to ensure a better consideration of electromagnetic boundary conditions. The weak electromagnetic formulation inducing integrals on open infinite domains, a Boundary Integral Method is used.

Introduction

Due to increasing needs in navigation, communication and surveillance, aircrafts are equipped with more and more antennal structures. The antennas implementation have drastic backlashes on aircrafts design. Indeed, antennas reduce aircraft aerodynamic efficiency, lead to structural weakening (drillings), cannot be placed anywhere and induce a weight increase through wiring. Thus, the decrease of antennas number and locations on future aircrafts is a crucial stake in aeronautics.

Many researches focused on embedding antennas on aircraft structural surfaces (fuselage, wings, tail, ...) have been performed over the last 15 years in order to improve both structural efficiency and antenna performances [START_REF] You | Design of load-bearing antenna structures by embedding technology of microstrip antenna in composite sandwich structure[END_REF]. These studies especially led in the development of microstrip antennas which are low profile and conformable to planar and warped surfaces. Hence, these antennas can easily be bounded on aircraft structural surfaces. Moreover, the emergence of metamaterials (new electromagnetic materials) will soon allow us to elaborate multifrequencies microstrip antennas.

In flight conditions, planes are subjected to aerodynamic loads leading to the structure distortion. Hence, the influences of the structural surfaces distortions on the microstrip electromagnetic response have to be studied, independently of antenna developments. The knowledge of mechanical strains could then be used for smart tuning of the antennas.

A few papers dealing with curved antennas can be found in the literature. In Antilla and Alexopoulos [START_REF] Antilla | Scattering from complex three dimensional geometries by a curvilinear hybrid finite element integral equation approach[END_REF], a 3D hybrid finite element-integral equation method has been developed in general curvilinear coordinates. In Jacobs [START_REF] Jacobs | Analysis of planar and curved microstrip antennas[END_REF], planar and curved radiative microstrip antennas are, for example, analyzed. The influence of the curvature on the radar cross section (RCS) is moreover shown for such antennas in Volakis et al. [START_REF] Volakis | Finite Element Method for Adn12 Electromagnetics. Antennas, Microwave Circuits and Scattering Application[END_REF]. Nevertheless, for cavities antennas, only simple configurations CONTACT Dr. F. Pablo fredericrf.1sirap-u@olbap.

Laboratoire d'Energétique, Mécanique et Electromagnétisme, Université Paris Ouest Nanterre La Défense, EA like cylindrical surfaces are considered in these works. In operating conditions, the strained surfaces of the antenna being complex, specific tools have to be developed. On an other hand, an overview of industrial coupling softwares (such as ANSYS, COMSOL) has shown that these tools are not well suited to solve this kind of coupled problem [START_REF] Reig | Monaco: Modélisation nano composants[END_REF] and [START_REF] Su | Overview of Electromagnetic Modeling Software[END_REF]. The novelty of this work is to mix the mechanical and electromagnetic problems using the same finite element in order to take into account the physical distortion of the antenna. The present work, focused on the development of a numerical tool able to predict the impact of microstrip distortions on electromagnetic fields, is part of the MSIE project (Intelligent Materials and Structures for Electromagnetism), launched in 2008 by the competitiveness french cluster AStech with the view of answering to the above problematic.

This paper presents the study of electromagnetic/mechanical coupling modeling for the analysis of microstrip antennas scattering. It derives from usual mechanical approach and specific mathematic tools for free space problems [START_REF] Volakis | Finite Element Method for Adn12 Electromagnetics. Antennas, Microwave Circuits and Scattering Application[END_REF], [START_REF] Bossavits | Computational Electromagnetism, Variational Formulation, Complementary, Edge Elements[END_REF]. From a 3D discretization of the structure, a classical mechanical finite nodal element model is coupled with an electromagnetic finite vector element one, developed by Jin [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF]. This leads to a particular finite element having both nodal mechanical and edges electromagnetic degrees of freedom. Hence, the electromagnetic/mechanical problem is solved using the same mesh. In order to solve mechanical problems with warped or distorted antennas, the considered finite element is an hexahedral one. Using the usual boundary integral method, the electromagnetic fields are obtained in the whole space.

In a first part, the microstrip antenna integrated on a composite structure is described and the governing equations are detailed. The finite element approximation is presented in the second part and deals with the electromagnetic field approximation, using the finite vector element approach and the boundary integral method. In the last section, this recent numerical tool is used through various configurations.

Problem description 2.1. Geometric description

The specific structure to be considered in the present work is illustrated in Figure 1. It consists in a rectangular microstrip patch antenna printed on a dielectric substrate. A ground plane is glued underneath the substrate in order to stop electromagnetic waves propagation. The antenna thickness being far smaller than the substrate one, it will be neglected in the following developments. This structure is then supposed to be integrated on composite aircraft panels.

The cavity model, described in Jin [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF] and illustrated in Figure 2 can be used to represent the above structure. The ground plane is then considered as an infinite metallic plane and the substrate is attended in all the cavity.

In the following developments, the infinity space located above the cavity and the ground plane, and the inside of the cavity will, respectively, be represented by V ∞ and V symbols.

Studying the antenna scattering, the cavity is assumed to be subjected to electromagnetic waves emitted by a source point of V ∞ .

Considering the Cartesian coordinate system (O, x, ŷ, ẑ) presented in Figure 3 r any observation point M(x, y, z) of the infinity space is located by the vector r = -→ OM whose euclidian norm is defined by r = r = x 2 + y 2 + z 2 , (note that M is located by M(r, θ, φ) in the spherical frame ( ρ, θ, φ)); Finally, it is assumed that the cavity is subjected to mechanical strains and loads. The governing equations and variational formulations are presented in the next subsections.

Mechanical governing equations and weak form

As previously underlined, the structural cavity V will be bounded on composite panels which are distorted by external loads. Hence, this distortion will affect V . The study of its mechanical strains is then essential.

Considering a static approach, the mechanical governing equations consists in the equilibrium equation, the constitutive equation in V and the strain-displacement relationship [START_REF] Fung | Foundations of Solid Mechanics[END_REF]:

-→ div T + f = 0, (1) T = a S( u), ( 2
)
S( u) = 1 2 (grad u + grad T u), (3)
where T, S( u), f and a are, respectively, the stress tensor, the strain tensor, the prescribed body forces vector applied to V and the elastic stiffness tensor. In addition to these equations, one has to take into account the mechanical boundary conditions: a first one where mechanical vector displacement u is imposed, and the second one where a surface force vector F is applied on the boundary ∂V F .

The classical weak form of the mechanical problem is then given by the following:

Find u ∈ U ad such that

V S * ( u * ) : T( u) dV = V u * . f dV + ∂V F u * . F dS, ( 4
)
for any admissible virtual displacement u * ∈ U * ad .

This variational principle is the usual starting point for any finite element approximations in mechanical modeling.

Electromagnetic governing equations

Considering a microstrip antenna, the electromagnetic governing equations in the V ∞ space and the cavity V consist in the Maxwell's equations and the constitutive relationships. In addition, the boundary conditions, radiation conditions, jump conditions and electromagnetic excitation are considered for the mathematical formulation.

The usual time harmonic approach is considered in the following developments. As a consequence, all field quantities are supposed to be harmonically oscillating functions with a single frequency ω [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF], [START_REF] Silvester | Finite Element for Electrical Engineers[END_REF].

The Maxwell equations

Considering the time-harmonic modeling, the Maxwell's equations can be written in a simplified form in V ∞ and V :

-→ rot E = -jω B, ( 5
) -→ rot H = jω D, ( 6
)
div D = ρ, ( 7
)
div B = 0, (8)
where E, D, ρ, H and B are, respectively, the electric field intensity, the electric flux density, the electric charge density, the magnetic field intensity and the magnetic flux density.

The constitutive relations

They describe the macroscopic properties of the mediums.

D = E, (9)
B = μ H, ( 10
)
where the constitutive parameters and μ are, respectively, the permittivity and permeability of the medium. These parameters are scalars for isotropic media.

In the present study, the cavity backed patch antenna is considered to be situated in free space. Then, V ∞ is characterized by the free-space permittivity ( 0 ) and permeability (μ 0 ). Moreover, the cavity is supposed to be filled with an homogeneous material represented by its relative permittivity r = / 0 and relative permeability μ r = μ/μ 0 . These two properties are here scalars but would have been tensors in the case of inhomogeneous filling material.

The boundary conditions

Metallic materials being considered as perfect conductors, the following conditions have to be verified on the ground plane and on the surface of the patch antenna:

n ∧ E = 0, ( 11
) n . D = 0, ( 12
)
where n is the unit normal vector of the considered surface.

The radiation conditions

If the outer boundary of a domain recedes to infinity, this domain is called unbounded or open. In this case, to make sure that the solution of the problem is unique, an additional condition has to be specified at this outer boundary. It is referred to a radiation condition as the Sommerfeld gauge [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF]:

lim r→∞ r [ -→ rot E + jk 0 r ∧ E] = 0, (13)
lim r→∞ r [ -→ rot H + jk 0 r ∧ H] = 0, ( 14
)
where k 0 is the free-space wavenumber (k 0 = 2π/λ 0 = ω √ 0 μ 0 ).

The jump conditions

In the present problem, the conservation laws imply to satisfy the following jump conditions across the interface between the two medias V ∞ and V [START_REF] Demkowicz | Modeling of electromagnetic absorption/scatering problems using hp-adaptive finite elements[END_REF]:

n . [[ D]] = 0, ( 15
) n . [[ B]] = 0, ( 16
) n ∧ [[ E]] = 0, (17) n ∧ [[ H]] = 0. (18)

Electromagnetic excitation

Studying the antenna scattering, an electric field intensity E inc is supposed to appear at a source point of V ∞ . It is considered equals to

E inc ( r) = E 0 e -j k 0 kinc . r , ( 19
)
where E 0 and kinc are, respectively, the amplitude and propagation vectors (see Figure 3). The analytical expression of E 0 , for a polarization angle denoted α, is given by

E 0 = R cs inc T E sph 0 ,
where R cs inc is the cross matrix between cartesian and spherical coordinates, and E sph 0 is the electric polarization field in spherical coordinates:

R cs inc = ⎡ ⎣ sin θ inc cos φ inc sin θ inc sin φ inc cos θ inc cos θ inc cos φ inc cos θ inc sin φ inc -sin θ inc -sin φ inc cos φ inc 0 ⎤ ⎦ E sph 0 = (0 cos α sin α) T .
The corresponding magnetic field intensity H inc equals to

H inc ( r) = 1 z 0 ( kinc ∧ E 0 ) e -j k 0 kinc . r , ( 20
)
where z 0 = √ μ 0 / 0 is the free space wave impedance. In the next sections, the usual approach which consists in solving the problem according to E is followed. Once E is established, the magnetic flux density can be obtained using the first relationship of Maxwell's equations.

Weak electromagnetic formulation

The main difficulty in the resolution of an electromagnetic problem for antennas consists in estimating integrals on an open infinite domain. Using numerical methods such as the Finite Element Method (FEM), the discretized domain to be considered has to be extended far from the source region in order to impose the radiation condition. This leads to huge computing times.

Meanwhile, it is possible to avoid this problem using the Finite Element-Boundary Integral (FEM-BIM) Method [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF], first developed in mechanical engineering and later introduced in electromagnetism.

This method consists in introducing a fictitious boundary which encloses the structures to be studied. Interior to this boundary, the FEM is used to formulate the fields, whereas in the exterior region, the fields are represented by a boundary integral. The fields of the two regions are coupled at the fictitious boundary via the field continuity conditions. This leads to a coupled system for solution of the interior and boundary fields.

The analytical developments of this electromagnetic problem are detailed in Jin [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF]. The fields formulation in the V ∞ space using the boundary integral method leads to introduce the scalar Green's function G 0 .

Then, the weak form of the electromagnetic boundary value problem can be presented in the following formulation [START_REF] Adnet | Modélisation numérique du couplage mécanique/électromagnétisme pour l'étude de la sensibilité du comportement électromagnétique d' antennes patch aux déformations mécaniques[END_REF]:

Find E ∈ E ad such that F v ( E, E * ) + F s ( E, E * ) = F ext ( E * ) (21)
for all admissible electric field E * ∈ E * ad , with

F v ( E, E * ) = V -→ rot E * ( r) 1 μ r -→ rot E( r) dV -k 2 0 V E * ( r) r E( r) dV F s ( E, E * ) = 2 S ap S ap div | r ( n ∧ E * s ( r))G 0 ( r, r )... div | r ( n ∧ E s ( r ))dS dS -2k 2 0 S ap S ap G 0 ( r, r )... ( n ∧ E * s ( r)).( n ∧ E s ( r ))dS dS F ext ( E * ) = -2 j k 0 S ap n ∧ E * s ( r) kinc ∧ E inc dS,
where S ap = S aperture denotes the area of the cavity which is not covered by the antenna (see Figure 2), n its unit normal vector, E s the surface electric field and G 0 ( r, r ) is the scalar Green's function, given by [START_REF] Tai | Dyadic Green's Function in Electromagnetic Theory[END_REF]- [START_REF] Conil | Propagation électromagnétique en Milieu Complexe[END_REF]:

G 0 ( r, r ) = e -jk 0 | r-r | 4π| r -r |
Note that the integrals are here reduced to the inside of the cavity and the aperture surface. Nevertheless, the resolution of this weak form, will also give the far fields E f ar and H f ar , using the following relationships [START_REF] Volakis | Finite Element Method for Adn12 Electromagnetics. Antennas, Microwave Circuits and Scattering Application[END_REF], ∀ | r| ≥ 1/k 0 :

E far ( r) = -z 0 kinc ∧ H far ( r) , ( 22
)
H far ( r) = j k 0 z 0 S ap G( r, r ) . n ∧ E s ( r ) dS , ( 23
)
where G( r, r ) is the dyadic Green's function, evaluated for far zone.

Apparently, the functional Eq. ( 21) only depends on the electric field. Consequently, it may be discretized. However, due to the singularity associated with the derivatives of the Green's function G 0 , the evaluation of the surface integral in F s ( E, E * ) leads to a new difficulty.

The finite element approximations

This section is dedicated to the finite element approximations.

In the present work, the mechanical behavior of substrate is supposed nonsensitive to the electromagnetic fields. Thus, a weak coupling between electromagnetism and mechanical behavior is assumed. At each time, the coupled problem can then be processed in the following two steps: r First, the mechanical part (Eq. 4) is solved using usual node-based discretization. The mechanical displacement is obtained by interpolating the nodal values of elements.

r Then, the electromagnetic problem (Eq. 21) is solved, tak- ing into account the mechanical strain of the structure. The specificity of the present method lies in the use of a single finite element having both nodal mechanical and edges electromagnetic degrees of freedom. The mechanical and electromagnetic problems are then solved using the same mesh.

Moreover, the antenna curvatures are time-dependent. This tool allows to predict how the electromagnetic response evolves at each time.

Discussion about the finite element shapes

In order to ensure the solution of the electromagnetic problem to verify the four Maxwell's equations in the three directions, and for a better consideration of electromagnetic boundary condition, the microstrip antenna needs to be discretized using 3D finite elements. For meshes compatibility reasons, the substrate is also discretized using 3D finite elements. Both mechanical and electromagnetic problem are thus solved in the 3D space.

Five main kinds of elements shape can be considered: brick, tetrahedra, prism, pyramid and hexahedra. Brick elements are suited to the study planar structures but they would lead to poor results for warped or distorted antennas, which is our specific study. Indeed, they are not able to automatically mesh arbitrary structures unless staircasing is permitted [START_REF] Volakis | Finite Element Method for Adn12 Electromagnetics. Antennas, Microwave Circuits and Scattering Application[END_REF]. Due to their sharp form, the deformed shape of prisms, pyramids and tetrahedral elements are known to lead to nonrealistic stress concentrations in mechanical studies. Finally, the hexahedral elements (or distorted brick element), able to automatically mesh arbitrary structures, are best suited to model the warped antennas.

The finite element approximation of the coupled problem is described in the following subsections.

The finite element approximation for the mechanical part

Implemented finite element

The element implemented in the software developed for the present study, is similar to the ANSYS SOLID45 element [17].

It consists in eight nodes linear brick element with three degrees of freedom at each node (three translations in the nodal directions).

The elemental matrices

Here, classical methods are used. Using the engineering formalism, the strain and stress tensors can be written as vectors:

S = S} = {S 1 S 2 S 3 S 4 S 5 S 6 } T T = {T } = {T 1 T 2 T 3 T 4 T 5 T 6 } T
Then, the stiffness tensor a can be written as a 6 × 6 matrix a.

Discretizing the studied structure using the FEM, nodal displacements and strains of any element are approximated by {u} = N {q e } and {S} = B {q e }, where the vector {q e } contains the element nodal displacements and N and B are 6 × 24 matrices, respectively, containing the interpolation functions and their derivatives.

Substituting these relationships in each integral of the weak formulation (4) restrained to elementary volume V e and surface S e , the following elemental stiffness matrix and loads vector are obtained:

K e = V e B T a B dV e , {F e } = V e N T { f } dV e + S e N T {F} dS e ,
where { f } and {F} are vectors containing the body and surface loads applied to the element nodes.

The discretized problem

Assembling the elemental matrices over the volume, the mechanical problem can be expressed in the usual following global form:

K {q} = {F mecha }, (24)
where {q} contains the mechanical displacement degrees of freedom of the structure and {F mecha } corresponds to the applied external mechanical loads. The stiffness matrix K being symmetric and the equation being linear, classical numerical schemes can be used to solve the mechanical problem.

The finite element approximation for the electromagnetic part: the use of vector finite elements

Node-based expansions for finite element solutions are suitable for modeling scalar quantities but typically not so for simulating electromagnetic vector fields. As a consequence, a new approach, recently used in electromagnetism, uses the vector finite elements or edge elements which assign degrees of freedom to the edges rather than to the nodes of the elements. These types of elements have been described by Krantz [START_REF] Krantz | Geometric Integration Theory[END_REF]. They have been shown to be free of the shortcomings of node-based expansions [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF], [START_REF] Andersen | Adaptive multiresolution antenna modeling using hierarchical mixed-order tangential vector finite elements[END_REF]. Especially, the occurrence of nonphysical solutions, also called spurious solutions and generally attributed to lack of enforcement of the divergence condition is avoided. Several investigations and applications can be found in the literature [START_REF] Volakis | Finite Element Method for Adn12 Electromagnetics. Antennas, Microwave Circuits and Scattering Application[END_REF], [START_REF] Bossavits | Computational Electromagnetism, Variational Formulation, Complementary, Edge Elements[END_REF], [START_REF] Andersen | Adaptive multiresolution antenna modeling using hierarchical mixed-order tangential vector finite elements[END_REF], and [START_REF] Nédélec | A new family of mixed elements in r 3[END_REF].

Vector finite elements

The electromagnetic problem represented by the weak formulation (21) will be solved using a classical hexahedral edge element presented in Figure 4.

Here, ξ , η and ζ represent the reduced coordinates, varying in [-1; 1] 3 , in the (O, ξ, η, ζ ) system of reference. Each edge is assigned by a number i: r i = 1 to 4 for edges parallel to ξ , r i = 5 to 8 for edges parallel to η, r i = 9 to 12 for edges parallel to ζ .

The electrical field in the element is then approximated by

E e = 12 i=1 -→ W v i (ξ, η, ζ ) E e i . ( 25
)
The functions -→ W v i are the basis functions and are given in Appendix A. The use of these basic functions leads to the following two important features to be underlined: r The continuity of the tangential field across all element edges is guaranteed.

r The divergence condition div -→ W v i = 0 within the element is satisfied. The Eq. ( 25) can thus be written in a following matrix form:

E e = ⎧ ⎨ ⎩ E e x E e y E e z ⎫ ⎬ ⎭ = W v {E e } (26)
with

W v = ⎡ ⎢ ⎣ -→ W v 1 . x -→ W v 2 . x ... -→ W v 12 . x -→ W v 1 . ŷ -→ W v 2 . ŷ ... -→ W v 12 . ŷ -→ W v 1 . ẑ -→ W v 2 . ẑ ... -→ W v 12 . ẑ ⎤ ⎥ ⎦ (27)
and {E e } is the vector of components E e i (size 12 × 1). Considering the weak formulation Eq. ( 21), two kinds of terms can be distinguished: terms depending on volume integrals and terms depending on surface ones. All integrals are estimated using the above hexahedral edge element. Meanwhile, it is necessary to impose one coordinate equals to one for surface integrals. The basis functions -→ W s i so obtained are gathered in the W s matrix.

The elemental matrices

Substituting Eq. ( 26) into each integral of Eq. ( 21) restrained to elementary volume V e and surface S e gives the following elemental matrices for a curved structure:

K v e = V e W T v ROT T μ r -1 ROT W T v dV e , M v e = V e
S e W s T N T k inc w R cs {E sph 0 } dS e ,
where ROT is the rotational operator:

ROT = ⎡ ⎣ 0 -∂/∂z ∂/∂y ∂/∂z 0 -∂/∂x -∂/∂y ∂/∂x 0 ⎤ ⎦ . ( 28
)
N is the matrix associated to n ∧ . with n = n x x + n y ŷ + n z ẑ:

N = ⎡ ⎣ 0 -n y n z n z 0 -n x -n y n x 0 ⎤ ⎦ ( 29
)
and DIV is the divergence operator, defined by

DIV = ∂/∂x ∂/∂y ∂/∂z . ( 30
)
N and DIV are, respectively, associated to the source point M . Their expressions are similar to N and DIV.

The discretized problem

Assembling the elemental matrices over the microstrip antenna volume, the electromagnetic problem can be expressed in a very global compact form as follows:

Y{q E } = {F} emag , ( 31
)
where {q E } contains the electric field degrees of freedom, and

Y = k 2 0 (M v + M s ) + (K v + K s ), {F} emag = jk 0 {F} ext .
M v , M s , K v , K s and {F} emag are the generalized matrix and vector.

This relationship is similar to the classical static mechanical one. In the electromagnetic case, the equivalent stiffness matrix is made of two terms:

r k 2 0 M v + K v
which is the contribution of the antenna structure discretized using the FEM. r k 2 0 M s + K s which is the contribution of the space (containing the electromagnetic sources), discretized using the Boundary Integral Method.

Matrices M v , K v , M s and K s being symmetric and Eq. ( 31) being linear, its resolution seems to be easy. Meanwhile, the evaluation of these last two matrices needs some attention: the green function, which has to be integrated, presents singularities at some particular points of the elements. This difficulty can be avoided using different Gauss points when integrating with respect to surfaces S and S .

Moreover, k 2 0 M s + K s is a full matrix. In order to reduce computing times, it is converted into a sparse matrix using usual numerical technics [START_REF] Volakis | Finite Element Method for Adn12 Electromagnetics. Antennas, Microwave Circuits and Scattering Application[END_REF].

Numerical examples

Our numerical tool was validated in Adnet [START_REF] Adnet | Modélisation numérique du couplage mécanique/électromagnétisme pour l'étude de la sensibilité du comportement électromagnétique d' antennes patch aux déformations mécaniques[END_REF]. In this section, it is used to model rectangular microstrip patch antennas, integrated on curved or strained structures. The first two tests deal with curved microstrip antennas on cylindrical and spherical surfaces. In the two last examples, the structure is subjected to mechanical loads. The impact of strains on the microstrip electromagnetic response is discussed.

Considering a scattering problem, the usual output parameter of interest is the RCS in the far field. Indeed, this parameter is an evaluation of the fields scattered from the structure: higher this coefficient is, more detectable is the structure. The RCS scattering is given by Refs. [START_REF] Volakis | Finite Element Method for Adn12 Electromagnetics. Antennas, Microwave Circuits and Scattering Application[END_REF], [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF]:

σ = lim r→+∞ 4π r 2 | H sca | 2 | H inc | 2 ,
where H sca is the scattered magnetic field intensity obtained by the far zone expression (Eq. 22).

It is usually represented in logarithmic coordinates and can be normalized with respect to the wavelength λ as follows [START_REF] Volakis | Finite Element Method for Adn12 Electromagnetics. Antennas, Microwave Circuits and Scattering Application[END_REF], [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF]: σ = 10 log 10 (σ ), (dBsm), σ = 10 log 10 σ λ 2 , (dB).

In the next applications, the influence of microstrip curvatures on electromagnetic performances will be hereafter analyzed through RCS diagrams.

A curved microstrip antenna on cylindrical surface

This test is about a patch antenna inserted in a metallic cylinder (Figure 5). It is vertically illuminated and two kinds of polarizations (vertical α = 90 • and horizontal α = 0 • ) are considered. This antenna has been studied by Kempel et al. [START_REF] Kempel | Scattering by cavity-backed antennas on a circular cylinder[END_REF] and Volakis The monostatic RCS are presented in Figures 8 and9 for vertical and horizontal polarizations. In spite of the use of different meshes, they are similar with results provided by Kempel et al. [START_REF] Kempel | Scattering by cavity-backed antennas on a circular cylinder[END_REF] and Volakis et al. [START_REF] Volakis | Finite Element Method for Adn12 Electromagnetics. Antennas, Microwave Circuits and Scattering Application[END_REF]. The case R = 2 m is similar to results obtained with a planar antenna. Near the resonance frequency (5.20 GHz for vertical polarization and 3.62 GHz for horizontal polarization), the level of RCS increases when radius decreases. These simulations show that the scattering of a curved antenna on cylindrical surface is upper than the one of a planar antenna. The antenna is thus failed.

A curved microstrip antenna on spherical surface

In this test, the patch antenna is inserted in a spherical metallic surface (Figure 10). It is vertically illuminated by an incident wave with two kinds of polarizations (vertical α = 90 • and horizontal α = 0 • ). Four radius values are considered: R = 0.05 m, 0.1 m, 0.2 m and 2 m. Mesh for R = 0.05 m is given in Figure 11. The RCS mono static are plotted Figures 12 and13. Conclusions are similar to the cylindrical case. At the resonance frequency, the antenna scattering strongly depends on the curvature. The antenna again becomes less adapted when R decreases.

The two previous tests deal with antennas on conformal configurations. The simulations show that the ability of hexahedral elements to take into account each node coordinates allows both planar and conformal configurations to be simulated using the same element. This is not true for brick elements: computing planar and conformal cases needs to, respectively, implement two brick elements: a planar element [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF] and a curved element [START_REF] Volakis | Finite Element Method for Adn12 Electromagnetics. Antennas, Microwave Circuits and Scattering Application[END_REF]. Moreover, the RCS predictions show the influence of the surface curvature on the antenna scattering. From these results, the study of mechanical behavior of the antenna and the evaluation of the impact of mechanical strains on the electromagnetic signals must be done. In the next tests, two antennas are considered, subjected to sinusoidal pressure. These simulations need mechanical solution from Eq. ( 24). The electromagnetism problem (Eq. 31) is solved taking into account the obtained mechanical strain of the structure. As specified previously, the same FE and, consequently the same mesh, are used for mechanical and electromagnetic problems.

A strained planar patch antenna

The following test is about a planar patch antenna, subjected to a sinusoidal pressure (Figure 14). This antenna is supposed to be a 0.03 m × 0.02 m rectangular patch lying on the top of a 0.06 m × 0.05 m × 7.874 10 -4 m dielectric filled cavity. The dielectric constant and permeability of the filling substrate are ε r = 2.17 and μ r = 1.00, and its Young modulus E is about 1 GPa. The structure is clamped at faces defined by x = -0.03 m and x = 0.03 m. Moreover, the structure is supposed to be vertically illuminated.

The mechanical pressure equals to

P(x) = P 0 sin π x L ,
where L is the length of the dielectric. The sensitivity of electromagnetic variables is studied for horizontal and vertical polarizations, with respect to two pressure magnitudes: case 1: P 0 = 0.025 N/m 2 and case 2: P 0 = 0.05 N/m 2 .

The maximal vertical displacement value is given in Figure 15 for case 1 and in Figure 16 for case 2. The RCS prediction are plotted Figures 17 and18. "Def1" and "Def2" means, respectively, results for cases 1 and 2.

The impact of mechanical strains on the electromagnetic behavior cannot be neglected. Especially for case 2, where a 24.7% RCS shift is noticed for a maximal displacement value equals to 0.0042 m. It induces an increasing of the antenna detectability.

A conformal antenna with five patches on a strained ring

This last test deals with a conformal antenna made with fives patches, embedded on a ring and subjected to mechanical loads. The geometry of the antenna is given in Figure 19. The thickness of the substrate is 0.8779 10 -3 m. The angle between two patches equals to 72 • . The dielectric constant equals to r = 2.17. Its young modulus is 1 GPa. The structure is excited by five cylindrical pressures P which are applied periodically on areas without patches (see two of them Figure 19). Their expression is

P(x) = P 0 sin π θ θ 0 , ( 32
)
where θ 0 is defined in Figure 19.

The structure is clamped on the lateral face, where y = -25.5 × 10 -3 m.

The mesh is presented in Figure 20, and the mechanical displacement induced is given in Figure 21. Maximal displacements are located on the free faces of the antenna and equal to 0.0082 m.

The antenna is vertically illuminated by a planar wave with horizontal polarization. Four frequencies are considered: 0.5 GHz, 1 GHz, 2 GHz and 4 GHz. For each case, bistatic RCS are plotted Figures 22232425, for initial and strained configurations. The impact of mechanical strains seems to be negligible. This can be explained by the small influence of mechanical displacement on the electrical field plotted Figures 26272829.

In conclusion, in this application, in spite of the maximal value of displacement, the mechanical loads do not affect the RCS curves and the efficiency of the antenna. The antenna is always in good running.

Conclusion

This work deals with the modeling of the electromagnetic patch antennas subjected to mechanical stresses. It has been motivated by the MSIE project from the competitiveness French Cluster Astech, where the development of metamaterial strip antennas for aeronautical application is evaluated. Such antennas which are thin and conformable can be embedded on any aircrafts surfaces, subjected to mechanical loads. Consequently, the influence of mechanical strains on electromagnetic signals has to be evaluated.

In this way, a recent numerical tool is presented in this paper, taking into account the weak coupling between electromagnetism and mechanical behavior to predict the microstrip antenna scattering. The same hexahedral finite element is developed with a classical nodal formulation for mechanical fields, while the electromagnetic ones are expressed by an edge formulation with vector finite elements.

From the weak formulation, an FEM/BIM approach is used to avoid the discretization of overall space. It consists in introducing a fictitious boundary that encloses the structure to be studied. Classical FEM is used to approximate the fields in the closed domain, whereas the fields in the opened region are evaluated by the BIM.

Then, the 3D finite element developed in this work allows to distort the antenna subjected to mechanical loads and to simultaneously obtain the modified electromagnetic fields (the antenna curvatures are time-dependent).

The numerical simulations presented in this paper show that the distortion of microstrip antenna can affect its scattering. This is the case of the three first examples. Moreover, sometimes, this impact can be negligible even if mechanical displacements are nonzero. For example, in the last test, the maximal displacement of the radius cylinder is close to 5.5%, but simulations show that the antenna is always well running.

Consequently, the study of the coupling between electromagnetism and mechanical behavior is necessary and insightful for antenna scattering. This is also the case for radiative antenna which is the subject of our next work.

Appendix A. The basis functions for the hexahedral edge element

The expression of the shape functions along each edge of the hexahedral element are [START_REF] Volakis | Finite Element Method for Adn12 Electromagnetics. Antennas, Microwave Circuits and Scattering Application[END_REF]

Figure  .

  Figure . The studied structure.

Figure  .

  Figure . The cavity model.

Figure  .

  Figure . The coordinate system. r the E inc and H inc field are emitted from a source point M , located by r = --→ OM , and received by the observation point M.Finally, it is assumed that the cavity is subjected to mechanical strains and loads. The governing equations and variational formulations are presented in the next subsections.

Figure  .

  Figure . The electromagnetic element.

  G 0 ( r, r ) N W s dS e dS e , DIV T G 0 ( r, r ) DIV N W s dS e dS e , {F ext e } = 2

Figure  .

  Figure . A curved microstrip antenna on cylindrical surface.

Figure  .

  Figure . A meshed curved microstrip antenna on cylindrical surface with R =  m.

Figure  .

  Figure . A meshed curved microstrip antenna on cylindrical surface with R = . m.

  .

Figure  .

  Figure . RCS of the antenna on cylindrical surface, vertically illuminated with α = 90 • .

Figure  .

  Figure . RCS of the antenna on cylindrical surface, vertically illuminated with α = 0 • .

Figure

  Figure . A curved microstrip antenna on spherical surface.

Figure

  Figure . A meshed curved microstrip antenna on spherical surface with R = . m.

Figure  .

  Figure . RCS of the antenna on spherical surface, vertically illuminated with α = 90 • .

Figure  .

  Figure . RCS of the antenna on spherical surface, vertically illuminated with α = 0 • .

Figure  .

  Figure . Planar patch antenna under sinusoidal pressure.

Figure  .

  Figure . Vertical displacement for P 0 = . N/m 2 .

Figure  .

  Figure . Vertical displacement for P 0 = . N/m 2 .

Figure  .

  Figure . RCS for P 0 = . N/m 2 .

Figure  .

  Figure . RCS for P 0 = . N/m 2 .

Figure

  Figure . A conformal antenna with five patches on a strained ring.

Figure  .

  Figure . Ring Mesh.

Figure  .

  Figure . Displacement field of the ring.

Figure  .

  Figure . Bistatic RCS for . GHz frequency.

Figure

  Figure . Bistatic RCS for  GHz frequency.

Figure  .

  Figure . Bistatic RCS for  GHz frequency.

Figure  .

  Figure . Bistatic RCS for  GHz frequency.

Figure  .

  Figure . Undeformed case for . GHz frequency.

Figure  .

  Figure . Deformed case for . GHz frequency.

Figure  .

  Figure . Undeformed case for  GHz frequency.

Figure  .

  Figure . Deformed case for  GHz frequency.

8 ( 1 + 8 ( 1 + 8 ( 1 +

 818181 η k η) (1 + ζ k ζ ) --→ grad ξ for edges parallel to ξ, ξ k ξ ) (1 + ζ k ζ ) --→ grad η for edges parallel to η, η k η) (1 + ξ k ξ ) --→ grad ζfor edges parallel to ζ , where ξ k , η k , ζ k denotes the edge's location in (O, ξ, η, ζ ) and l e k denotes its length.


B W Jervis

PhD C Bigan

PhD M W Jervis

BSc M Besleaga

New-onset Alzheimer's Disease and Normal Subjects 100% Differentiated by P300

Keywords: Alzheimer's disease, biomarkers, P300, auditory evoked potential, independent components analysis, artificial neural network, PSFAM, Bayes classifier

Previous work has suggested that evoked potential analysis might allow the detection of subjects with new onset Alzheimer's Disease, which would be useful clinically and personally. Here, it is described how new-onset Alzheimer's disease subjects have been differentiated from healthy, normal subjects to 100% accuracy, based on the back-projected independent components (BICs) of the P300 peak at the EEG electrodes in the response to an oddball, auditory evoked potential paradigm. After artefact removal, clustering, selection, and normalisation processes the BICs were classified using a neural network, a Bayes classifier, and a voting strategy. The technique is general and might be applied for pre-symptomatic detection, and to other conditions and evoked potentials, although further validation with more subjects, preferably in multi-centre studies is recommended.

Introduction

It is important to be able to diagnose new onset Alzheimer's Disease (AD) sufferers both for their care and their personal planning. Evoked potential analysis might provide a relatively inexpensive, quick, and non-invasive technique for this and has therefore been investigated. A method of distinguishing with 100% accuracy between early-stage AD patients and normal, healthy subjects (normals) based upon the non-oscillatory, independent components of the P300 peak in the P300 waveform elicited by an auditory oddball paradigm is described in this paper. Because averaging is not used, potentially significant components, unsynchronised to the stimulus, are not reduced, and results may be obtained using fewer trials per subject. This method could be a useful tool to aid diagnosis, and the selected independent components may be regarded as biomarkers. The method might also be useful for pre-symptomatic testing for Alzheimer's Disease (AD).

Since the review and description of work previous to 2011 [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF] , there have been further publications on the topic in question. In a review of 2011 [START_REF] Olichney | Cognitive Event-Related Potentials: Biomarkers of Synaptic Dysfunction Across the Stages of Alzheimer's Disease[END_REF] it was concluded that the sensitivities of a number of ERP components have great promise for the detection of the stages of Alzheimer's disease. Another review in 2014 [START_REF] Cintra | P300 evoked potential and risk of mild cognitive impairment progression to AD: a literature review[END_REF] was focussed on the progression from mild cognitive impairment (MCI) to AD. All the studies quoted followed changes in amplitude and latency of the P300 peak, but on an averaged basis. In reference [START_REF] Cecchi | A Clinical Trial to Validate Event-Related Potential Markers of Alzheimer's Disease in Outpatient Settings[END_REF] trial averaging and statistical analysis of the peak ERP amplitudes and latencies derived from a three-stimulus auditory oddball paradigm showed that the P3a and P3b peaks produced the most sensitive and reliable measures of the cognitive deficits associated with early Alzheimer's disease. None of this work [START_REF] Olichney | Cognitive Event-Related Potentials: Biomarkers of Synaptic Dysfunction Across the Stages of Alzheimer's Disease[END_REF][START_REF] Cintra | P300 evoked potential and risk of mild cognitive impairment progression to AD: a literature review[END_REF][START_REF] Cecchi | A Clinical Trial to Validate Event-Related Potential Markers of Alzheimer's Disease in Outpatient Settings[END_REF] addressed the analysis on a single trial basis as described here. However, Ouyang et al [START_REF] Ouyang | Reconstructing ERP amplitude effects after compensating for trial-to-trial latency jitter: A solution based on a novel application of residue iteration decomposition[END_REF] have analysed single trials by applying the technique of residue iteration decomposition (RIDE) to identify the latencies of the different ERP peaks in different trials. It seems that this technique, though, does require some averaging of trials to obtain the initial most likely latency of the ERP peaks, after which the individual latencies are found by an iterative method. No application which differentiated between different subject groups was presented. By contrast, in our work we derive the individual components of the individual ERP peaks comprising each individual single trial using independent components analysis (ICA), and apply this knowledge to differentiate between normals and ADs. It seems none of the authors [START_REF] Olichney | Cognitive Event-Related Potentials: Biomarkers of Synaptic Dysfunction Across the Stages of Alzheimer's Disease[END_REF][START_REF] Cintra | P300 evoked potential and risk of mild cognitive impairment progression to AD: a literature review[END_REF][START_REF] Cecchi | A Clinical Trial to Validate Event-Related Potential Markers of Alzheimer's Disease in Outpatient Settings[END_REF][START_REF] Ouyang | Reconstructing ERP amplitude effects after compensating for trial-to-trial latency jitter: A solution based on a novel application of residue iteration decomposition[END_REF] were aware of the previous work by both ourselves and those we quoted [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF] , although the work of Jung et al using ICA is mentioned in one paper [START_REF] Ouyang | Reconstructing ERP amplitude effects after compensating for trial-to-trial latency jitter: A solution based on a novel application of residue iteration decomposition[END_REF] . In another review [START_REF] Vecchio | The Use of Auditory Event-Related Potentials in Alzheimer's Disease Diagnosis[END_REF] it was concluded that MCI patients had prolonged P300 latencies compared to controls, but shortened P300 latencies when compared to AD patients meaning that ADs had longer latencies than normals.

Our research was carried out using the data obtained in earlier work, which has been thoroughly described in two previous publications [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF][START_REF] Jervis | The independent components of auditory P300 and CNV evoked potentials derived from single trial recordings[END_REF] . Thus, only the essentials of that work are repeated here. A selective analysis of that data using an artificial neural network, the Probabilistic Simplified Fuzzy ARTMAP (PSFAM) [START_REF] Jervis | The Probabilistic Simplified Fuzzy ARTMAP (PSFAM)[END_REF] , and a voting strategy is presented.

The undulatory P300 waveform includes a number of positive and negative peaks [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF] . By using independent components analysis (ICA) and back-projecting the nonoscillatory, independent, source signals to the scalp electrodes, it was found that the peaks in the P300 waveform consisted of many short duration, randomly occurring, and randomly positive or negative half-sinusoidal pulses [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF] . Here, attention is focussed upon the positive back-projected independent components (BICs) centred on the P300 peak, because the shape of the peak is primarily determined by these, and the latency of this peak is delayed in ADs compared to normals [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF][START_REF] Vecchio | The Use of Auditory Event-Related Potentials in Alzheimer's Disease Diagnosis[END_REF] . Therefore, these BICs were deemed the most likely to be useful for differentiating between ADs and normals.

Theoretical aspects

The voltage measured at each electrode depends upon the contributions there from all the independent cortical signal sources. These depend upon the unknown source signals and their unknown transmission paths from the sources to the electrodes. Fortunately, the individual source signals may be computed from the measured scalp voltages using ICA [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF] , where it was explained that if be a matrix of temporally independent source signals and be the matrix of measured signals at the electrodes, which are assumed to consist of linear sums of the source signals ), which have passed through an unknown, linear transmission system characterised by an m x m mixing matrix, , then we may write Thus, the estimated source signals may be found since can be found. Selected estimated source signals may then be multiplied by the mixing matrix to obtain their estimated contributions at the measurement electrodes, . These are referred to as the back-projected independent components (BICs). Thus, (3) The BICs are correct in both magnitude and sign, and so may be compared.

The PSFAM, used to classify the data, consisted of a Simplified Fuzzy ARTMAP (SFAM) and a Bayes classifier [START_REF] Jervis | The Probabilistic Simplified Fuzzy ARTMAP (PSFAM)[END_REF] . The latter produced the Bayes posterior probability that the test vector belonged to the class AD or class normal.

Measurements

Six male and three female normal, healthy subjects and two male and seven female newly diagnosed, early stage, mildly cognitively impaired AD subjects participated in auditory-evoked oddball P300 recordings as fully detailed previously [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF][START_REF] Jervis | The independent components of auditory P300 and CNV evoked potentials derived from single trial recordings[END_REF] . The ADs were under various drug treatments [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF] , where age effects are also discussed.

Scalp voltages were recorded at 27 standard electrode sites (see below). The voltage waveforms were sampled at 1024 Hz, the high pass cut-off frequency was 0.016Hz and the low pass cut-off was at 60Hz. A notch filter eliminated the mains frequency of 50 Hz. There were 40 target tones of 2 kHz and 160 non-target tones at 1 kHz. The inter-stimulus interval was 1.29 s. The subjects had closed eyes, were relaxed, and responded to the target tones by button-pressing. For each subject 360 target stimuli were recorded, with 600 pre-stimulus samples and 700 post-stimulus samples.

Procedures

The following signal processing was performed as fully detailed before [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF][START_REF] Jervis | The independent components of auditory P300 and CNV evoked potentials derived from single trial recordings[END_REF] .The independent components (ICs) of the P300 waveforms were obtained by applying Principal Components Analysis (PCA) first and then Independent Components Analysis (ICA) [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF] . These ICs were then back-projected to the measurement electrodes as the BICs. These were separated into separate bins centred around the P300 peaks. The highest variance BICs were selected for further processing. The BICs in each bin were clustered in two stages using the k-means clustering algorithm [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF][START_REF] Jervis | The independent components of auditory P300 and CNV evoked potentials derived from single trial recordings[END_REF] . In the primary stage clustering was by amplitude and latency; in the secondary stage by the scalp topographies [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF][START_REF] Jervis | The independent components of auditory P300 and CNV evoked potentials derived from single trial recordings[END_REF] . Noise components were eliminated by filtering out ICs according to the number of zero crossings in their waveform and their largest and smallest amplitudes [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF][START_REF] Jervis | The independent components of auditory P300 and CNV evoked potentials derived from single trial recordings[END_REF] . Within each bin the peak amplitudes, latencies, and the scalp topographies of the BICs were saved for analysis [START_REF] Jervis | Waveform Analysis of Non-Oscillatory Independent Components in Single-Trial Auditory Event-Related Activity in Healthy Subjects and Alzheimer's Disease Patients[END_REF][START_REF] Jervis | The independent components of auditory P300 and CNV evoked potentials derived from single trial recordings[END_REF] .

In the previous paper 1 , the BIC results obtained at this stage of processing were discussed, and have also been briefly reviewed in the above Introduction. Here we describe that data in detail and how it has been processed further to allow identification of the individuals with newly diagnosed Alzheimer's disease. This data may be requested from the corresponding author.

The data spreadsheet was 39 columns wide and contained 5,302 rows. The data in each row included subject details, subject class (AD or normal), trial number, BIC information (which bin, which cluster, positive or negative, amplitude, latency), and the voltages of the BICs at the twenty seven measurement electrodes used which were Fp1, Fp2, F7, F8, F3, F4, FC5, FC6, FC1, FC2, T7, T8, C3, C4, CP5, CP6, CP1, CP2, P7, P8, P3, P4, O1, O2, Fz, Cz, and Pz. These twenty seven BIC voltages, taken in the above order, comprise the BIC topology vector. This spreadsheet was divided into separate spreadsheets for the AD and normal subjects.

Data preparation for PSFAM

It was intended to use the data in the above two spreadsheets to train classifiers to distinguish between the AD and normal subjects. These 5,302 x 39 data contained some personal details which were irrelevant to this training, and so these columns were ultimately deleted. Since it had been established 1 that those positive BICs associated with the P300 peak, and their latencies, were the most significant in distinguishing the two classes, the amount of data could be considerably reduced by using only that for positive BICs found close to the P300 peak, i.e. those in bin 5. This reduced the data arrays to 520 and 581 rows for the AD and normal subjects respectively. Since latency was more important than amplitude, the amplitude column was also deleted. The training and test vectors then consisted of the latency vector and the topology vector, making a 1 x 28 row vector per trial.

The PSFAM neural network required input data normalised to between 0 and 1. The latency column was normalised by dividing all values by the largest. Some of the elements of a topology vector could be positive, while others were negative, reflecting the scalp voltage topography. A simple formula was applied to convert the values over the positive and negative voltage range of the topology vector array to values between 0 and 1. A spreadsheet function was used to set all data to numerical values, because the classifiers could not accept exponentials. A sample section of a resulting data array up to column F4 is shown in Table 1. The values are as calculated, but of course were not measured to the accuracy shown. The contents of the first column define the bin and cluster numbers of the BIC in a particular row. Subj is the subject number. Trial is the trial number for that subject. Latencies are in ms, and voltages in mV. Files in training and validation data formats were derived from the two data arrays by deleting columns 1 to 3, and replacing column 4 (trial) by the subject class; 0 for a normal, and 1 for an AD subject. Test vectors may be formed by deleting the first four columns.

A training file was constructed consisting of the data for the first five normal subjects and for the first five AD subjects. A validation file was constructed from the last four normal subjects and the last four AD subjects. Thus all the data was used. The structures of these two files were identical, and so their training and validation roles were reversible. These files contained the data obtained from all 40 trials. In clinical practice, where patients may not be sufficiently co-operative, it is necessary to use fewer trials. For this reason we also explored the use of just 10 trials and of just five trials. We had observed that five trials were necessary to ensure that a BIC was found (remember the random nature of the appearance of the BICs). Thus additional training and validation files were produced by eliminating the row data for trial numbers greater than 11 and greater than 6 using functions of the spreadsheet.

PSFAM procedures and considerations

A number of tests were carried out using the training and validation files to determine the optimal values of the PSFAM parameters. These were found to be: vigilance, ρ = 0.65; global smoothing parameter, σ = 0.02; the remaining parameters set to 1. In the training mode the data were presented in random order using the "shuffle" button. The normal subjects were assigned to the reference class 0, and the AD subjects to class 1. Each subject was represented by several row vectors, which

Insert Table 1 contained the values calculated for different trials and clusters. The classification of an individual row for a given subject could be correct or incorrect. When the overall accuracy per subject was investigated the SFAM was found to result in the higher accuracy for the normal subjects, but the AD subjects were most accurately classified by the Bayes method. The reason for the difference lies in the different classification techniques used by these classifiers. In the SFAM the degree of fuzzy membership of the test input vector to the fuzzy power set of the weight vector is calculated using the match function [START_REF] Jervis | The Probabilistic Simplified Fuzzy ARTMAP (PSFAM)[END_REF] where (4) and ˄ is the fuzzy logic operator. Thus, each input vector's similarity to each weight vector is determined and the test input vector is assigned to the class of the weight vector for which , i.e. for which it is acceptably close. In the Bayes classifier the summed differences of the test input vector to all the training vectors of the normal class is compared to those of the AD class, according to (5) where are the numbers of training vectors for the normal and AD subjects respectively, the indices refer to the normal and AD classes, is the assumed normal input test vector, the and the are the training vectors, is the smoothing parameter, and denotes the vector transpose. If the equality is satisfied, the test input vector is assigned to the class normal. This difference in the classification methods explains why in general the two classifiers do not always assign the same class to a test vector. The discrepancies between the two classifiers can be expected to be greater when the training and test vectors have more random properties as in the case for the BICs of both the normal and the AD subjects obtained in this work. Further, it was also suspected that the ADs' BICs were more random than those for the normal subjects. In fact it was found that the SFAM gave the higher classification accuracy for the normal subjects, while the Bayes classifier gave the higher classification accuracy for the AD subjects. The correct class for a subject was indicated by the classifier which classed the most input vectors for that subject as being of the same class. The classification of a subject was therefore decided by adopting this voting strategy. When the number of test vectors assigned to the two classes was equal, they both indicated the correct class.

Another test was undertaken to establish whether fewer, carefully chosen electrodes might be used. Thus classification using only columns Lat, Fp1, Fp 2 , P3, P4, Fz, and Pz was attempted, but the results were worse than when all 28 columns were used. In another test, the training and validation files were exchanged, when equally good results were obtained.

PSFAM results

Forty trials

Table 2 gives the overall percentage classification accuracies for training and testing the PSFAM with the full set of data from 40 trials when the training and validation files were prepared as described above. It is seen that the higher percentage correct classification of normal subjects is achieved by the SFAM (77%:45%), and that of the AD subjects by the Bayes classifier (62%:40%).

The percentage classification accuracies of the individual subjects, by their individual input test vectors, are shown in Table 3. It is seen that, if the percentage correct in the SFAM column is greater than that in the Bayes column, the subject is normal; if the converse is true, the subject is in the AD class. This is the basis of the voting strategy, which yields 100% correct classification of the newly diagnosed AD subjects in this study from the positive BICs at the 27 EEG electrodes centred on the P300 peak. Of course, slightly different numerical values are obtained when the training is repeated, owing to the shuffling of the training data, but the conclusions remain unaltered.

Ten trials

Similar results were obtained as for the case of the 40 trial data when only the first 10 trials were used. Table 4 gives the subject classification accuracies. The same conclusions apply as for 40 trials, with the addition that in two cases, N14 and AD40, the percentage accuracies are the same for both classifiers, but they both predict the same correct class. We also see that the accuracy of classification of the normal subjects by the SFAM, and that of the AD subjects by the Bayes classifier has increased with the reduced number of trials. Some of this may be attributed to the shuffling of the input vectors, but it seems more likely it could be owing to a reduction in the degree of randomness associated with using fewer vectors.

Table 5 presents a sample of the numbers of 0s and 1s output by the classifiers to represent the classes of the input vectors for the subjects, normal and AD respectively, and the classification voted, which is correct in each case.

Insert Table 2 Insert Table

Five trials

Similar results were obtained again in the case of 5 trials as shown in Table 6. Apart from subject N15, the classification accuracies were again improved by this further reduction in the number of trials included. However, the use of fewer trials still is likely to lead to more classification failures when there may be too few input vectors for reliable testing or even the absence of any vector.

Discussion

It has been clearly demonstrated in this research that the positive voltage BICs associated with the P300 peak may constitute an excellent biomarker for new onset AD, since 100% accurate differentiation between new onset ADs and normals was achieved. It is quite possible that they could indicate AD pre-symptomatically. This could be tested by making measurements on subjects at-risk of AD, such as carriers of the apolipoprotein E4 gene with a family history of AD, or subjects for whom synaptic dysfunction has been detected by elevated CSF phosphor-tau 2 . The technique is non-invasive, requires a reduced number of trials, is inexpensive, and can be employed in any hospital EEG department. Because of the small sample size it is desirable that far more subjects be tested, and preferably in multi-centre studies, to validate it, if it be considered useful. The digitised multi-centre recordings could be processed centrally to ensure conformity. This validation could take place during the clinical studies. Such studies might also be used to investigate the effects of drug treatment, the relationship of the results to those of similar studies on subjects with other neurological diseases to reduce the risk of misdiagnosis, and the possible usefulness of the BICs associated with other peaks in the waveform. Extension to other conditions such as Parkinson's disease and other evoked potentials is also a possibility.

Conclusions

New-onset ADs were differentiated from normals to 100% accuracy by classifying the positive voltage, back-projected, independent components centred on the P300 waveform peak response to an oddball task, auditory evoked potential using a simplified fuzzy ARTMAP neural network, a Bayes classifier and a voting strategy. It may also be possible to detect AD pre-symptomatically, but more, preferably multicentre research on more subjects is necessary to validate the technique. Extension to other conditions and evoked potentials is a possibility.

Table 1

 1 Sample section of data array

	Insert Table 5
	Insert Table 6

Table 2

 2 Overall percentage classification accuracies by the two classifiers for normal and AD subjects.

	Subject	SFAM	Bayes
		% correct	% correct
	N14	61	53
	N15	69	41
	N16	75	39
	N20062	84	41
	AD36	52	90
	AD38	39	86
	AD40	34	82
	AD43	43	77

Table 3

 3 Percentage classification accuracies for individual test subjects based on their test vectors using 40 trials.


Sofiane Schaack
email: schaack@insp.jussieu.fr

Philippe Depondt

Fabio Finocchi

H-bond symmetrization in high pressure methane hydrate

First-principle molecular dynamics simulations of methane hydrate MH-III, including the quantum properties of the hydrogen nuclei, were carried out at pressures in the 5-65 GPa range, in order to observe the H-bond symmetrization at high pressure and at room temperature. According to our simulations, the symmetrization transition takes place around 40 GPa and is little dependent on isotope substitution. We find that, consistently with the rather complex crystal structure of MH-III, the transition is much more convoluted in the hydrate than in ices VII and X. In methane hydrate, due to the presence of non equivalent O ions with distinct O-O distances, the dissimilar H bonds symmetrize in a pressure domain rather than at a single critical pressure.

Introduction

Methane hydrate (MH) is a clathrate compound in which methane molecules are trapped within a solid cage of water [START_REF] Sloan | Clathrate hydrates of natural gases[END_REF][START_REF] Demirbas | Methane gas hydrate: as a natural gas source[END_REF]. MH could have occurred in the nebula from which ice giants Uranus and Neptune were formed and is still present on Titan [START_REF] Loveday | [END_REF]4]. Significant deposits of methane clathrate hydrates have been also found on Earth, in cold upper sediments on the sea and ocean floors, which could be destabilized as the oceans warm up and release large quantities of greenhouse methane gas [5]. The previous facts motivated many researches on the stability of methane hydrates under various conditions.

The interaction between the host molecules and the cage is essentially repulsive; this is why methane hydrates were initially thought to be stable only at relatively modest pressure, below a few GPa. However, quite recently [6,7], several phases of MH, differing by structure and composition, were found to be stable at very high pressures, beyond 10 GPa and up to ∼ 100 GPa. At those pressures, not only could the interaction between the cage and the methane molecules be significantly altered as ice cages shrink, but the ice cage itself could undergo deep modifications. Indeed, the hydrogen bonds between the water molecules that form the cage may become symmetric and the ice cage eventually evolves towards a non-molecular form, which modifies the mechanical properties of the hydrate, as well as the molecule-cage interaction.

The hydrogen bond symmetrization under pressure was first observed in ice phase X. Along the transition from phase VII (H-bonds, with proton disorder) to phase X (symmetric bonds, proton ordered) the quantum properties of the hydrogen nuclei play an important role [8,9]; the most striking one regards the computed transition pressure at room temperature, which is reduced from approximately 100 GPa within the classical frame down to 65 GPa, in agreement with the experimental findings, when nuclear quantum effects are included. Among the latter, the zero-point energy (ZPE) is crucial in two respects: first, at the onset of proton tunneling between the two potential wells that correspond to the asymmetric H bond; second, in lowering the critical pressure for symmetrization with respect to the classical frame. In methane hydrate, a similar process is expected in the high-pressure MH-III phase, which is the stable phase of the system at ambient temperature and above 2GPa.

The paper is organized as follows: after a brief description of the simulation methods, in section 3 we first describe the order parameter that is commonly used to describe the H-bond symmetrization; then, we analyze the displacements of the proton along the O-O axis or normally to it (section 4); finally, we investigate the role of the interaction between the protons (deuterons) in the cage with the guest methane molecules in Section 5. All along the paper, we compare the behavior of the protons in MH-III to prototypical pure ice.

Computational details

Constant volume, constant temperature (NVT) molecular dynamics simulations at room temperature were carried out using the Quantum Thermal Bath method (QTB) [10] in order to include nuclear quantum effects. The QTB method is implemented in the Quantum Espresso [11] package. It is based on a Langevin equation in which the random force R(t) is tailored to fulfill the quantum fluctuation-dissipation theorem [12]: instead of a white noise, one uses a colored noise with the following power spectrum: R(ω)

2 = 2mγhω   1 2 + 1 exp hω k B T -1   ( 1
)
where m is the nucleus mass, k B the Boltzmann's constant, T the temperature, ω/2π the frequency, and γ is the friction coefficient entering also the dissipation term in the Langevin equation. QTB is strictly exact only for harmonic systems. While dealing with anharmonicity, semi-classical methods such as QTB are submitted to zero-point energy leakage (ZPEL) which tends to equilibrate the energy distribution following the classical equipartition theorem.

Adjusting the friction coefficient γ (eq 1) allows to lower the effect of the ZPEL [13] as the larger γ is the more the classical system is coupled to the quantum bath. Despite the widening of the vibrational peaks, the QTB has proven valuable in several systems [9,14,15,16] and the QTB simulations presenting ZPEL have shown that they provide correct mode frequencies [17]. Rather importantly in the present context, a systematic comparison with the classical case can be carried out by simply replacing the colored QTB noise in equation ( 1) with the standard white noise R(ω) 2 = 2mγ k B T , while all the other computational ingredients remain identical [START_REF]hω 2 accounting for the zero point energy when hω k B T → 0[END_REF].

The initial structure was taken from ref. [7]. The simulated samples consist of 16 methane molecules and 32 water molecules in a 2 × 1 × 2 orthorhombic supercell. The methane molecules are hydrogenated (CH 4 ) while the ice cage is considered with hydrogenated (H 2 O) or deuterated (D 2 O) water molecules. The electronic structure and atomic forces were obtained within the Generalized Gradient approximation (GGA) [START_REF] Perdew | [END_REF] to the DFT. Van der Waals interactions were added by following the semi-empirical scheme by Grimme [20]; after extensive tests, a better agreement with experimental lattice constants and compressibility could be obtained with respect to the pure GGA. Ultra-soft pseudo-potentials were used to describe the interaction between the ionic cores and the valence electrons: a plane wave expansion cutoff of E cut = 40 Ry ensures total energy and atomic force convergence. The friction coefficient was set to γ = 1 THz limiting the impact of ZPEL effects in the simulations which tend to underestimate NQE. Their impact was estimated to be at maximum ∼ 7% of the temperature of the different degrees XXIXth IUPAP Conference on Computational Physics CCP2017 IOP Conf. Series: Journal of Physics: Conf. Series 1136 (2018) 012018 IOP Publishing doi:10.1088/1742-6596/1136/1/012018 3 of freedom of MHIII. Simulations were run at constant volume, with lattice parameters chosen to sample pressures from 5 up to 65 GPa by obtaining an isotropic stress tensor within our statistical error. Oxygen and carbon atoms were initially set at their crystallographic positions, while cage hydrogen atoms were let to relax during short simulations with large friction coefficients γ in order to explore different configurations and start with variable orientations for the CH 4 molecules. All atoms are free to move during the simulations. The typical duration time of the simulations is 30 ps.

Characterizing the symmetrization transition

Along the VII-X transition of ice under pressure, the effective potential that is felt by the proton in hydrogen bonds is usually described in terms of a one-dimensional double-well potential, along the O-O axis [21,22]. As pressure increases, the oxygen-oxygen nearest neighbor separation continuously decreases: the distance between the double wells and the barrier is accordingly reduced, making it possible for the proton to tunnel back and forth between the two wells. Upon further contraction of the O-O distance, the proton has its largest probability halfway between the two O ions. The simplest order parameter and the most widely used [8,9,14,16] to estimate the location of a hydrogen nucleus between the two nearest oxygen atoms is:

χ = d(O [1] -H) -d(O [2] -H) (2)
where d(O [n] -H) is the distance between the selected proton and its n th oxygen nearest neighbor. So defined, this parameter equals zero for proton-symmetric hydrogen bonds. The order parameter χ was observed to depend on the O [1] -O [2] distance d(O [1] -O [2] ) in several H-bonded systems: below a threshold (d(O [1] -O [2] ) 2.42 Å), the probability density P (χ) is maximum at χ = 0 and the H-bond becomes symmetric [22,21]. As the O-O distances in MH-III are shorter than in ice VII at all pressures, it is reasonable to expect the symmetrization transition to happen at a lower pressure than the ice VII → X transition. Although direct measurements of the symmetrization of the hydrogen bonds in MH-III are still lacking, recent Xray diffraction and Raman spectroscopy experiments [23,24,25] indicate a phase transition from MH-III to an High Pressure (HP) phase between 40 and 55 GPa although a precise description of the transition was not provided.

We show in figure 1 the probability distributions of the order parameter χ, averaged over all O-H..O bonds ( top panel ), and for each different OO distance ( bottom panel ), in MH-III as computed from our simulations. P (χ) evolves from a distribution with two maxima, which is characteristic of asymmetric hydrogen bonds, to a single-peak distribution, which denotes the symmetric hydrogen bonds, around 40 GPa. This points towards a symmetrization transition of the ice cage in MH-III around 40 GPa, therefore at a much lower pressure than in pure ice. One can also note that even at a quite low pressure the distribution does not vanish for χ = 0 (P (χ = 0) = 0) showing that the ice skeleton is proton disordered even at P ∼ 10 GPa.

However, the situation is more complex in MH-III than in pure ice: first, P (χ) behaves rather similarly in both hydrogenated and deuterated cages. This would imply that isotope effects are practically negligible, and that nuclear quantum effects are not very relevant, in striking contrast with the case of ice VII → X transition [22,9,6]. Second, there are three non equivalent O-O distances in MH-III (2.811, 2.807 and 2.783 Å detailed in figure 2a) as experimentally measured at 3 GPa [6]. The mean O-O distances along the dynamics can indeed be collected in three distinct groups, two of which ((a), corresponding to O1-O1 and oriented along c; (b), connecting O1 to O6 in the notation of ref. [START_REF] Loveday | [END_REF]) have the same mean value d(O [1] -O [2] ) while the third group (c, corresponding to O1-O2 and O1-O2+, at the vertex of the 6-fold water rings) shows a slightly smaller mean O-O distance at all pressures. The corresponding probability distributions for the order parameter P deuterated and hydrogenated ice skeletons, respectively, while P (c) (χ) is flat around χ = 0. The difference between the distributions reveals that the transition takes place in two steps: at slightly lower pressures, the hydrogen bonds are symmetrized along the shorter O-O axes in group (c), and those in the remaining groups follow, with a delay corresponding to a few GPa. The picture that emerges from the P (α) (χ) distributions is that of a symmetrization transition that extends over ∼ 5 GPa around 40 GPa rather than happening at a precise critical pressure, as in the case of pure ice [9]. In order to confirm the indications that we obtained from the trends of the order parameter χ, we computed the vibrational spectra through the Fourier transforms (FT) of the time correlation functions. In particular, we examined the velocity-velocity time correlation functions ṙα (t) ṙα (0) , where r α is the α component of vector joining the proton (deuterium) in the ice cage and its oxygen first neighbor O [1] -H. We then considered the projections of the previous vectors parallel or perpendicular (r and r ⊥ , respectively, detailed in Figure 2b) to the O [1] -O [2] axis. The spectra of the parallel components ṙ show a shallow minimum of the mode's frequency around 35-40 GPa, which is consistent with the softening of the shorter H bonds (figure 3). Such a softening is a signature of the symmetrization transition, as discussed in the case of pure ice [9,22], which points towards a critical pressure around 35-40 GPa, which is slightly smaller but consistent with the indications from the distribution of the corresponding order parameter χ. As for the P (χ) distributions, the spectra for hydrogenated and deuterated ice cages little differ once the latter have been renormalized by the mass ratio √ 2, confirming that nuclear quantum effects do not play a crucial role in the symmetrization transition of MH-III. As it downshifts, the OD stretching mode peak intensity flattens, thus the uncertainty on the mode's frequency IOP Conf. Series: Journal of Physics: Conf. Series 1136 (2018) 012018 IOP Publishing doi:10.1088/1742-6596/1136/1/012018 6 prevents the difference of P (χ) distributions in terms of the three OO distances (Figure 2a) to be resolved. The distinction between the proton motion along or normally to the O-O axis translates also in the distributions of the proton (deuterium) localization, which we analyze in the following.

Analysis of the proton delocalization

Using simulation data from [9], we compare the hydrogen distribution with that of ice under pressure in which the VII→X transition shows important nuclear quantum effects, in constrast with methane hydrate.

First, we show (figure 4) the hydrogen distribution for both ice and methane hydrate with the same scales and normalization, for comparison. One observes that the OH r distance is approximately the same in both compounds as are the parallel widths, while the perpendicular extension is larger for the hydrate at all pressures. In order to obtain a more quantitative result one may compute the mean value of r and parallel σ and perpendicular σ ⊥ standard deviations (figure 5). This clearly confirms a broader distribution in the normal direction than along the O-O axis, and suggests that a purely onedimensional model hamiltonian, like those usually adopted [22,21] are unable to catch the complexity of the hydrogen bond in methane hydrate or other highly anisotropic systems. A finer analysis based on the differences of P (r , r ⊥ ) between the hydrogenated and deuterated cages shows that the radial delocalization is larger in the case of (H theoretical analysis of normal modes and accurate Raman measurements, it can be deduced that the interactions between the water cage and the methane molecules play a significant role in the system's dynamics, starting at a pressure of 20GPa [26]. Figure 4 shows the probability P (r , r ⊥ ): the normal contributions to the distribution are broad with a significant amplitude up to r ⊥ ≈ 0.4 Å. The inverted comma shape of the distribution indicates that the hydrogen atoms, as they move about, tend to retain a relatively constant O-H distance. As pressure is increased, the shape of the distribution changes: the parallel contribution is sensitive to increased tunneling and thus broadens until approximately 40GPa while the normal contribution tends to shrink slightly.

Molecule-cage interaction

A step further in the analysis between the ice cage and the methane molecules is taken with the conditional probability ρ(d HH |r ⊥ ) (resp. ρ(d DH |r ⊥ )) which is computed for the cage hydrogen (deuterium) nuclei: d HH (d DH for deuterium) is the distance between a cage hydrogen nucleus and the closest hydrogen nucleus of a methane molecule. The point is now to assess how the interaction of the methane molecules with the cage acts upon the lateral motion of the cage hydrogen atoms. Figure 6 shows a clear correlation between the presence of a methane hydrogen and the sideways displacement of the cage hydrogen (deuterium) atoms for all examined pressures and without clear indication of a change at the alleged transition pressure of 40GPa. When protons in the CH 4 molecules approach the cage hydrogen (deuterium) below a typical distance (∼ 2.0 Å at P=10 GPa, ∼ 1.9 Å at P=35 GPa and ∼ 1.7 Å at P=65 GPa), r ⊥ tends to increase, in order to maximize d HH (resp. d DH ). This provides an insight into the reason why no clear isotopic effect can be seen at the transition pressure: we know that a small asymmetry in a double well potential can severely reduce tunneling between the two wells [14], now our results show that the effective double well in which the cage hydrogen nuclei reside is perturbed by the molecular hydrogen nuclei. This shows that the normal delocalization is a result of the repulsive interaction between the hydrogen (deuterium) in the cage and those in the methane molecule.

Other simulations were carried out, removing the methane molecules while keeping the oxygen atoms fixed to their original positions and letting the cage hydrogen atoms move freely: the result thereof is that these hydrogen atoms escape from the O-O axes, which shows that the repulsive interactions between the hydrogen of the host CH 4 molecules and the cage hydrogen (resp. deuterium) is also necessary for the stability of the H bonds in the latter frame. molecules is to keep the hydrogen atoms close enough to the straight line between the oxygen atoms, but not on it.

Concluding remarks

According to previous suggestions, we confirm that, assuming the MH-III structure remains stable at these pressures, a symmetrization transition should occur at approximately 35-40 GPa, much below that in ice VII→X transition due to the smaller O-O distances in methane hydrate with respect to pure ice at all pressures. Moreover, in MH-III, we observe two main groups of hydrogen bonds corresponding to distinct O-O distances, which symmetrize at slightly different pressures. This behaviour contrasts with that of the ice VII where all hydrogen bonds are equivalent and symmetrize at the same pressure. Appropriate spectral experiments could provide a confirmation thereof. Another difference between ice and MH-III is the isotope effect. While in pure ice the H → D substitution significantly upshifts the symmetrization transition pressure, we found that the symmetrization of the hydrogen bonds occurs at almost same pressures in the hydrogenated and deuterated cages of methane hydrate. Both the lower symmetry of MH-III with respect to pure ice and the interaction between the confined methane molecules and the protons (deuterons) in the H bonds are at the root of a more complex behavior. In particular we observe in methane hydrate an increase of the lateral disorder of the cage hydrogen (deuterium) atoms, as compared with ice VII and X, which points towards a reduction of nuclear quantum effects. As a matter of fact, protons in MH-III cages are more delocalized than deuteriums in both longitudinal and normal directions along the O-O axis. Our results thus exhibit a rather complex and/or subtle

Figure 1 .Figure 2 .

 12 Figure 1. Distributions of the order parameter P (χ) at several pressures, for the hydrogenated (left) and deuterated (right) ice skeleton. Upper panels: P (χ) is computed by integrating on all hydrogen bonds. Bottom panels: Distributions of the order parameter P (χ) for each group of O-O equivalent distances (a, b, c) (see text).

Figure 3 .

 3 Figure 3. Fourier transforms in the frequency domain of the parallel component of the velocityvelocity time-correlation functions ṙ (t) ṙ (0) along the O-O axis. Thick line are for D 2 O-CH 4 , thin lines for H 2 O -CH 4 . The frequencies of the latter spectra were divided by √ 2 in order to take into account the mass difference.

Figure 4 .

 4 Figure 4. Probability density P (r , r ⊥ ) for (H 2 O) 2 : CH 4 (left) (D 2 O) 2 : CH 4 (middle) and ice (right) systems at three pressures : one below the transition, one not very far from it and one above. Pressure values have been selected in order that the mean OO distances < d oo > in the three systems match, as the hydrogen bond symmetrization process is controlled by this parameter.

Figure 5 .

 5 Figure 5. Ice and methane hydrate systems : (left) mean r , (middle) parallel and (right) perpendicular standard deviations.

Figure 6 .

 6 Figure 6. Conditional probability ρ(d HH |r ⊥ ) for a hydrogenated cage (left panels) or deuterated cage (right panels).

  XXIXth IUPAP Conference on Computational Physics CCP2017 IOP Conf. Series: Journal of Physics: Conf. Series 1136 (2018) 012018 IOP Publishing doi:10.1088/1742-6596/1136/1/012018 9 behavior of the cage hydrogen atoms close to the symmetrization transition and should hopefully constitute a strong enticement for experimental endeavours in the same pressure range.

Acknowledgments

We would like to thank L.E. Bove, U. Ranieri and R. Gaal for useful discussions. We acknowledge critical reading of the manuscript by Simon Huppert. This work was granted access to the HPC resources of CINES under the allocation A0010906719 made by GENCI.