Richard Renou
email: richard.renou@cea.fr

Laurent Soulard

Emilien Lescoute

Corentin Dereure

Didier Loison

Jean-Pierre Guin

Silica glass structural properties under elastic shock compression: Experiments and molecular simulations

Introduction

High power lasers are widely used in industrial processes for material design and in optical systems. One major problem is the damage an optical component can encounter under intense irradiation : density changes, spallation or crack initiation. Damaged component could initiate damages in other optical components leading to a global failure of the whole optical system. Understanding the damaging mechanisms occurring under shock loading conditions is crucial. Amorphous materials do not exhibit a structural organization at long range. However, in silica glass, a two-level organization can be observed. Oxygen atoms are structured around silicon atoms in a tetrahedral fashion. The Short-Range Order (SRO) in silica is related to the tetrahedral properties (Si-O bonds, intra-tetrahedral angle). In addition, Oxygen atoms are usually connecting several tetrahedra and rings composed of 3 to 10 tetrahedra are observed. The Medium-Range Order (MRO) corresponds to the rings properties (Si-Si distribution, inter-tetrahedral angle). Free volumes can also be observed in between rings due to the amorphous nature of glass materials. This peculiar two-level organization has been widely studied in the past. However little is known about the evolution of the silica structure under shock loading conditions.

The overall behavior of fused silica glass under shock compression has been studied for decades. First experiments of shock produced by explosive devices with free surface velocity measurement were published by Wackerle in 1962 [START_REF] Wackerle | ShockWave Compression of Quartz[END_REF] . In 1970, Barker published results of plate impact experiments on fused silica 2 . In these two papers, maximum pressure generated into the sample was about 10 GPa. From the early 1980s to the present day, other works have dealt with the shock compression behavior of fused silica under higher pressure condition up to 80 GPa [START_REF] Marsh | LASL Shock Hugoniot Data[END_REF][START_REF] Chhabildas | In Shock Waves in Condensed Matter[END_REF][START_REF] Alexander | Changes to the Shock Response of Fused Quartz Due to Glass Modification[END_REF] .

Accessing the atomic structure inside glasses by experimental means is quite challenging. Usual methods for crystalline phases (for example X-ray diffraction) are not able to give reliable results for amorphous material. Raman spectroscopy however can be used to investigate small rings inside silica glasses and the inter-tetrahedral OSiO angle distribution. Larger rings and many other structural properties remains unavailable though.

To overcome these limitations, numerical simulations are often carried out to investigate properties unavailable by experiments.

Molecular dynamic simulations is a powerful technique to gain insight on the microscopic properties of a material. Several works using non-equilibrium molecular dynamics simulations (NEMD) have been previously reported on fused silica under shock loading.

Early works were published by Kubota et al. in which the authors studied a 1D-shock generated by a rigid wall impacting a fused silica system [START_REF] Kubota | Structural Modifications in Fused Silica due to Laserdamage-induced Shock Compression[END_REF][START_REF] Dávila | Transformations in the Medium-range Order of Fused Silica under High Pressure[END_REF] . The authors focused on the analysis of the ring distribution in the densified silica glass. Su et al. (from 3000 to 5000 m/s) on a silica surface [START_REF] Su | Molecular Dynamics Simulation of Shock Induced Ejection on Fused Silica Surface[END_REF] . More recently, the transition of fused silica and α-quartz toward a high-pressure polymorph (stishovite) was thoroughly investigated by Wang et al. [START_REF] Wang | Atomic Scale Modeling of Shock Response of Fused Silica and α-Quartz[END_REF] . Some alternative methods are able to reproduce shocked states without using direct shocks. Izvekov et al. [START_REF] Izvekov | Mechanism of Densification in Silica Glass under Pressure as revealed by a Bottom-up Pairwise Effective Interaction Model[END_REF] investigated silica shocked states using the adptative Erpenbeck procedure [START_REF] Erpenbeck | Molecular Dynamics of Detonation. I. Equation of state and Hugoniot Curve for a Simple Reactive Fluid[END_REF] where a succession of simulations in the isothermal-isobaric ensemble (NPT) is used to converge slowly toward a shocked state. The authors essentially focused on the structural properties of silica glass in the densification regime and under hydrostatic compression. Barmes et al. 12 studied the angular and ring distribution of densified states using both NEMD and a method called Hugoniostat. In the latter case, the motion equations are rewritten and the Hugoniot relations are used as a constraint. The volume is chosen at the beginning of the simulation and the system will thus converge toward a shocked state. Barmes et al. showed that this Hugoniostat method is a good alternative to NEMD simulations and could be successfully applied to glass materials 12 .

Despite the number of studies on fused silica, data in the elastic regime are somehow scarce. The aim of this paper is four-fold : 1) to provide experimental data on the shocked states of silica glass in the elastic regime, 2) to use these experimental results for molecular model calibration and validation, 3) to ensure that the new Hugoniostat method developed in a recent work [START_REF] Maillet | Sampling Constraints in Average: The Example of Hugoniot Curves[END_REF] is yielding correct results on glassy systems and 4) to investigate the structural silica glass properties under shock compression using molecular simulations. In this work, we report both experimental and molecular dynamic modeling results on the elastic behavior of fused silica glass under laser driven shock conditions.

Experimental Study

A series of experiments have been performed on the Transportable Laser Shocks Generator (GCLT) Facility of the CEA. Figure 1 beam, which properties are as follow: a 1053 nm wavelength, a 4 to 100 ns duration, and an energy up to 50 J. The sample is a 10x10mm parallelepiped rectangle of fused silica glass, with a thickness of 2 mm, cut from a fused silica glass rod (spectrosil R 2000)

whose density is 2.2 g.cm -3 . All surfaces are polished down to the micron level to reduce surface defects and to have a good optical transparency for instrumentation purpose.

A 1 µm-thick aluminum is vapor deposited on the front face of the sample in order to promote a controlled and known interaction between the laser pulse and the matter. The later produces a high-pressure plasma whose expansion generates a shock wave that propagates into the sample. During the experiment, the sample is placed in a secondary vacuum chamber in order to avoid laser breakdown. The time resolved instrumentations used for each shot consists in: a transverse shadowgraphy to observe the shock wave while it propagates across the transparent sample [START_REF] Ecault | Observation of the Shock Wave Propagation Induced by a High-Power Laser Irradiation into an Epoxy Material[END_REF][START_REF] Loison | Soda-lime Glass Behavior under Laser Shock[END_REF] and a Velocity Interferometer System for Any Reflector (VISAR) apparatus 2,16 to measure the velocity of the interface between the aluminum layer and the fused silica glass accelerated by the shock near the front face.

Regarding the transverse shadowgraphy technique, a continuous laser is collimated and crosses through the transparent fused silica glass to reach 2 ICCD cameras (4Picos from Stanford Computer Optics able to take 2 pictures with a minimum interframing about 500 ns). The propagation of the shock wave induces a pressure variation, which modifies the optical refractive index of the glass. The collimated light of the laser beam is therefore deflected in the shocked area, which creates a shadow contrast on the ICCD camera. The velocity of the shock wave is computed from the measure of the displacement of the shock wave within a controlled time gap. This time gap as to be kept as short as possible due to the action of the release waves that can modify the shock. The 4Picos cameras were trigged by a Stanford DG535 to obtain two images of the shock propagation with a small interframe, typically 50 ns. In order to limit motion blur, which is an issue to capture fast displacement, the acquisition time was set to 5 ns. Figure 2 presents several shadowgraph images of a sample irradiated by a 123.7 GW.cm -2 intensity and 10 ns duration laser shot. The first picture is captured before the laser shot and used as a reference image. The two other images are obtained at 25 ns and 75 ns respectively after the beginning of the laser matter interaction.

The diameter of the laser spot is usually chosen greater than the sample thickness in order to neglect two-dimensional effects [START_REF] De Rességuier | Experimental and Numerical Study of Laser Induced Spallation in Glass[END_REF] . In this work, we chose to work with larger samples for convenience and front shock velocities were extracted from a small region of the sample, depicted here as the inset of Figure 2 and on a very short duration of 50 ns. In these conditions, the 2D effects can be neglected and we assume a one-dimensional shock propagation.

Molecular Dynamics Simulations

Silica pair potentials

Two pair potentials were chosen to model silica glass behavior in this work: the Beer-Kramer-Sanden [START_REF] Van Beest | Force Fields for Silicas and Aluminophosphates based on ab initio Calculations[END_REF] (BKS) and the Carré-Horbach-Ispas-Kob [START_REF] Carré | New Fitting Scheme to Obtain Effective Potential from Car-Parrinello Molecular-Dynamics Simulations: Application to Silica[END_REF] (CHIK) model. The BKS model was originally designed by quantum calculations to reproduce the quartz behavior at ambient temperature [START_REF] Van Beest | Force Fields for Silicas and Aluminophosphates based on ab initio Calculations[END_REF] . This potential has been used in many works and is still the reference potential for silica materials. A parametrization of this model was designed by Carre et al. [START_REF] Carré | New Fitting Scheme to Obtain Effective Potential from Car-Parrinello Molecular-Dynamics Simulations: Application to Silica[END_REF] simulations [START_REF] Car | Unified Approach for Molecular Dynamics and Density-Functional Theory[END_REF] . The choice of these two potentials was based on several comparative works.

Cowen et al. [START_REF] Cowen | On Force Fields for Molecular Dynamics Simulations of Crystalline Silica[END_REF] have proved that the BKS potential was able to reproduce the equation of state of the high-pressure silica polymorph (stishovite) in a better way compared to some other potentials. According to Soules et al. 22 , the CHIK and BKS potentials were able to reproduce the glass behavior, whereas the potential developed by Takada et al. [START_REF] Takada | Molecular Dynamics Simulations of Vitreous Silica Structures[END_REF] failed to reproduce the densification mechanisms. Recently, Carre et al. [START_REF] Carré | Developing Empirical Potentials from ab initio Simulations: The case of Amorphous Silica[END_REF] The Journal of Physical Chemistry some more sophisticated potentials obtained by force matching methods were unable to reproduce the silica glass structure compared to the simple CHIK model. Finally, Horbach et al. [START_REF] Horbach | Molecular Dynamics Computer Simulation of Amorphous Silica under High Pressure[END_REF] have shown that the CHIK potential was able to reproduce the experimental equation of state of silica glass along different isotherms until 10 GPa. To conclude, the CHIK and BKS potentials are able to reproduce fairly well the overall silica glass structure, better than many other potentials. Potentials based on reactive force fields could have been more accurate than simple pair potentials. However those reactive potential are known to be resources consuming. In this work we choose to focus on low computational cost potentials to be able to simulate huge systems (several millions particles).

The CHIK and BKS potential are taking into account electrostatic and Van Der Waals interactions, Equation 1. A Buckingham potential is used here to model the Van Der Waals contribution.

U ij = U vdw ij + U elec ij = A ij exp(-B ij r ij ) - C ij r 6 ij + q i q j 4π 0 r ij (1)
U is the potential, A, B, C are Buckingham parameters, q is the atomic charge, r is the distance between particles i, j. One major drawback of Buckingham-based potentials is the divergence for small distance values. This issue was addressed by adding a simple twoorder polynomial as described by Barmes et al. 12 . The new Van Der Waals contribution is therefore given by Equation 2.

U vdw ij =        A ij exp(-B ij r ij ) - C ij r 6 ij r ij > r 0 ij a ij r 2 ij + b ij r ij + c ij r ij ≤ r 0 ij (2)
a, b, c are the two-order polynomial coefficients. r 0 is the cut-off distance between the two mathematical forms. The coefficients were adjusted in order to insure the continuity of potentials and forces and to fit smoothly the original potentials around the cut-off distance.

The BKS and CHIK parameters used in this study are gathered in Table 1. The atomic 7

Page 7 of 52

ACS Paragon Plus Environment

The Journal of Physical Chemistry charges are q Si =2.4 and q O =-1.2 for the BKS model, and q Si =1.910418 and q O =-0.955209 for the CHIK model. The two-order polynomial coefficients for the BKS model were taken from the work of Barmes et al. 12 . The coefficients for the CHIK model were developed in this work. Figure 3 presents the new CHIK model with corrected potentials. The dotted lines show the original non-corrected potentials diverging for small distances.

Table 1: BKS and CHIK parameters used to model silica glass. A, B and C are the Buckingham parameters taken from their original works [START_REF] Van Beest | Force Fields for Silicas and Aluminophosphates based on ab initio Calculations[END_REF][START_REF] Carré | New Fitting Scheme to Obtain Effective Potential from Car-Parrinello Molecular-Dynamics Simulations: Application to Silica[END_REF] . a, b, c are two-order polynomial parameters used to correct the diverging potential at small distances.

BKS CHIK Si-O O-O Si-Si Si-O O-O A (J) 2.88411E-
(J.m -2 ) 5.93250E+02 6.17500E+01 2.39000E+01 7.39000E+02 3.15000E+02 b (J.m -1 ) -1.96148E-07 -3.61561E-08 -1.36192E-08 -2.38283E-07 -1.36590E-07 c (J) 1.61081E-17 5.21439E-18 1.90447E-18 1.88894E-17 1.51907E-17
Electrostatics interactions are usually taken into account with the so-called Ewald summation [START_REF] Ewald | Die Berechnung Optischer und Elektrostatischer Gitterpotentiale[END_REF] . However, another method with an interestingly low computational cost may be used for silica materials. This method, called Wolf method [START_REF] Wolf | Exact Method for the Simulation of Coulombic systems by Spherically Truncated, Pairwise r-1 Summation[END_REF] , is a shifted and truncated coulombic potential given by Equation 3.

U wol f ij = q i q j 4π 0 ( 1 r ij - 1 r 0 w + r ij -r 0 w r 0 w ), r ij ≤ r 0 w (3)
r 0 w is the cut-off distance usually set at 10.17 Å. The potential decreases therefore smoothly to 0 when approaching r 0 w . Despite the significant error, this method has been proved to give acceptable results with silica glass [START_REF] Carré | Amorphous Silica Modeled with Truncated and Screened Coulomb Interactions: A Molecular Dynamics Simulation Study[END_REF] and is a good alternative to the Ewald method.

We have checked (not shown) that the results obtained with both a Wolf method and an Ewald summation were identical under shock loading conditions. The Wolf method yields acceptable results on the pressure range tested.

8

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Modelling glass materials

The initial amorphous structure is usually obtained using a thermal cycle [START_REF] Vollmayr | Cooling-rate Effects in Amorphous Silica: A Computer-simulation Study[END_REF] . This method is able to generate a silica glass material with a structure close to the experimental one with a low percentage of defects [START_REF] Gedeon | Molecular Dynamics of Vitreous Silica -Variations in Potentials and Simulation Regimes[END_REF] . In this work, a β-cristoballite configuration was first generated and equilibrated at 300 K during 50 ps in the canonical ensemble (NVT) using a Langevin thermostat [START_REF] Grest | Molecular Dynamics Simulation for Polymers in the presence of a Heat Bath[END_REF] . The system was then heated at very high temperature (8000 K)

then cooled down to 300 K by decreasing the temperature by 500 K every 50 ps. Finally, the density was adjusted with a simulation in the isothermal-isobaric ensemble (NPT) at 300K using the so-called Nose-Hoover barostat 32 . The final densities were of 2.203 g.cm -3 and 2.260 g.cm -3 for the CHIK and BKS respectively which is very close to the experimental value of 2.200 g.cm -3 29 . According to Bidault et al. [START_REF] Bidault | Deformation of Silica Glass studied by Molecular Dynamics: Structural Origin of the Anisotropy and non-Newtonian Behavior[END_REF] , the unit cell should have a minimum length of 50 Å in order to obtain good silica ring statistics. We choose here to start with a cubic unit cell of 71 Å (24000 atoms) to be able to obtain 14-membered rings.

Simulating shocks with molecular dynamics simulations

Investigating material behaviors under direct shock loading requires systems of huge dimension in order to catch all mechanisms involved. The unit cell of 24000 atoms was therefore replicated 160 times along the shock direction (in this work the x-direction).

The system final dimensions are 1200x7x7 nm [START_REF] Marsh | LASL Shock Hugoniot Data[END_REF] which approximately corresponds to 4 million atoms. Direct shock loading were simulated using Non-Equilibrium Molecular Dynamics simulations (NEMD). In this case the shock was generated with a silica glass system impacting a hard immobile wall with a high velocity u p0 , as depicted in Figure 4.

Periodic boundary conditions in the y and z directions were used to simulate an infinite system and thus a 1D-shock is studied here. Free boundary conditions were imposed along the shock direction. All NEMD shock simulations were performed in the microcanonical ensemble (NVE) and were run during 500 ps. 9

Page 9 of 52

ACS Paragon Plus Environment

The Journal of Physical Chemistry

An alternative to direct shocks is the use of the Hugoniostat method. In this case the motion equations are rewritten so the system is naturally converging toward a shocked state. In this work we used the description established by Maillet et al. [START_REF] Maillet | Sampling Constraints in Average: The Example of Hugoniot Curves[END_REF] where the system is explicitly described by the variables (q, p, T) where q is the particle positions, p the momentum and T the temperature. q and p are updated using a regular Langevin scheme.

The temperature is updated using the Hugoniot relation as described by Equation 4.

T(q) = 1 Nk B 2 4 -c (< E > -U(q) + 1 2 (P(q)+ < P >)(1 -c)V) (4)
The bracket < ... > denotes a statistical average over the different configurations. k B is the Boltzmann constant, E the total energy, U the potential energy, P the total pressure, V is the initial volume and c is the compression rate. The temperature is computed using the Hugoniot relation and the positions and momentum are related to temperature changes by the Langevin scheme insuring a canonical distribution over time. The choice of a given compression ratio at the beginning of the simulation will therefore lead the system to converge toward the corresponding shocked state. Very small systems (3000 particles) may be used with this method reducing significantly the computational cost.

However, one drawback is that the material experiences a thermodynamic path slightly different compared to the one observed in a direct shock. The successful application of this method on another material (or a different phase) may not be guaranteed. Moreover

Results and Discussion

Hugoniot Curves

Assuming a 1D shock propagation, velocities, pressure and densities are related via mass and momentum conservative equations, respectively Equations 5a and 5b.

ρ 1 (u s -u p1 ) = ρ 0 (u s -u p0 ) (5a)
P 1 -P 0 = ρ 0 (u s -u p0 )(u p1 -u p0 ) (5b)
In these equations, index 0 and 1 refer to the state before and after the shock propagation respectively. u s is the shock wave velocity, u p is the velocity of the matter, P is the pressure and ρ the density.

The experimental shock front velocities were calculated using shadowgraphs and the velocity of the matter at the aluminum layer-silica boundary induced by the shock is obtained from VISAR measurements. Figure 5 represents the time evolution of the velocity of the interface for the shot reported in Figure 2. A double-wave structure is observed from both the shadowgraphs and the VISAR measurements. The first wave corresponds to the ramp wave setting the material in motion up to 67 m.s -1 . The second wave is depicted by the peak where the particle velocity reaches a maximum value of 277 m.s -1 . Interestingly, a delay is observed between the two waves as shown by the plateau between the ramp and the shock wave. The origin of this delay is uncertain at the moment. Wave shapes are correlated to several parameters like laser shot duration and sample thickness [START_REF] Davis | Laser-Shock-Induced Spall and the Intrinsic Strength of Glass[END_REF][START_REF] Wang | Laser-induced Decompression Shock Development in Fused Silica[END_REF] and additional data are required to accurately investigate this phenomenon. We will address this question in a future work.

In the elastic regime, silica glasses exhibit a particular behavior as reported by Alexander et al. for fused silica, borosilicate and soda-lime glasses [START_REF] Alexander | Changes to the Shock Response of Fused Quartz Due to Glass Modification[END_REF] . From the velocity profile of 11

Page 11 of 52

ACS Paragon Plus Environment

The Journal of Physical Chemistry the free surface, two waves can be distinguished : an elastic ramp wave and an elastic shock wave. Ramp waves were thoroughly studied theoretically and experimentally by Wang et al. [START_REF] Wang | Laser-induced Decompression Shock Development in Fused Silica[END_REF] and Davis et al. [START_REF] Davis | Laser-Shock-Induced Spall and the Intrinsic Strength of Glass[END_REF] . Decaying waves were reported in the numerical studies of Su et al. [START_REF] Su | Molecular Dynamics Simulation of Shock Induced Ejection on Fused Silica Surface[END_REF] and Luo et al. [START_REF] Luo | The Relation Between Shockstate Particle Velocity and Free Surface Velocity: A Molecular Dynamics Study on Single Crystal Cu and Silica Glass[END_REF] . In the elastic regime, fused silica exhibit a non-linear behavior and the compressibility increases under compression which prevents the formation of standard shock waves. The final wave shapes, ramp wave or decaying wave, depends on many parameters like shock pressure, laser fall time and sample thickness.

The particular double-wave structure is not observed in NEMD simulations. Examples of the time evolution of particle velocities in the material is shown Figure 6.a and Figure 6.b respectively extracted from the 2.07 GPa and 7.07 GPa NEMD simulation with the CHIK potential. A ramp wave is clearly observed in the first case and a decaying wave is seen in the latter case. However, the two-wave structure was not obtained in our simulations.

Su et al. in their work defined two velocities to investigate decaying waves [START_REF] Su | Molecular Dynamics Simulation of Shock Induced Ejection on Fused Silica Surface[END_REF] . The head velocity corresponds to the velocity of the wave propagating at 0.9 u p . The tail is the velocity computed at u p =200m/s. In our work, the tail velocity was chosen to be 0.1 u p .

The head and the tail of the ramp wave are propagating at a constant velocity of 6272m/s.

The front shape does not change over time. In the opposite, the decaying front becomes larger as the compression wave is crossing the material. The head and tail velocities of the decaying wave were estimated at 6317 m/s and 5765 m/s respectively. Pressures and densities can be extracted from molecular simulations, allowing us to use the Hugoniot relations to calculate the shock front velocity. Using Equation 5b the front is propagating at 5949 m/s which is close to the velocity of 5992 m/s computed at 0.5 u p . In the following, we chose to use the Hugoniot relations to estimate the shock velocity as it should yield a fair estimation of the front velocity. However, the definition of the front velocity is questionable in the case of decaying waves.

In order to compare both simulations and experiments, two methods were used to obtain the experimental shock front velocity. In the first method, only the second wave 12

Page 12 of 52

ACS Paragon Plus Environment

The Journal of Physical Chemistry velocity was considered as it corresponds to the maximum pressure. In the second method, two shocks were explicitly taken into account and the second wave velocity was computed using the first wave velocity. The differences between the two methods was less than 1%.

The shock pressure, obtained through Equation 5b, does not depend on the wave shapes.

The numerical and experimental shocked states can be safely compared regardless of the transient dynamic mechanisms involved. Only results obtained with the first method are presented in the following.

Experimental and numerical results are gathered in Table 2 and the Hugoniot curves are shown Figure 7.a and Figure 7.b. Our experimental work is labeled here GCLT. The direct shocks with NEMD simulations were carried out with the two BKS and CHIK models.

The Hugoniostat method was only tested with the CHIK potential. Other data taken from works of reference [START_REF] Wackerle | ShockWave Compression of Quartz[END_REF][START_REF] Alexander | Changes to the Shock Response of Fused Quartz Due to Glass Modification[END_REF] were also added on the Hugoniot curves. The pressure range studied here lies between 1.35 GPa and 8.71 GPa as the Hugoniot Elastic Limit (HEL) is around 8 GPA [START_REF] Sugiura | Dynamic Response of Fused Quartz in the Permanent Densification Region[END_REF] and the scope of the article is to investigate the structure under elastic shocks.

First of all, the experimental and numerical Hugoniot curves are in agreement with the literature. The fact that our experimental results are very close to the literature and simulated 1D shock proves our 1D hypothesis acceptable in this pressure range. Both NEMD and Hugoniostat simulations are able to reproduce the behavior of fused silica under elastic shock loading although the simulated thermodynamic path is different.

Interestingly, the simulations with the CHIK potential give superior results compared to the BKS potential. The relation between the pressure and the particle velocity is correct but the pressure-density formulation shows a significant discrepancy for the BKS potential.

This result is not a surprise since the BKS fails to reproduce several silica glass properties like the density at ambient temperature [START_REF] Takada | Molecular Dynamics Simulations of Vitreous Silica Structures[END_REF] .

The differences between experimental and numerical front velocities are large. However, comparing the two data sets is a difficult task especially with decaying wave front.

The front shapes are different and the two-waves structure was not observed in our simula- [START_REF] Horbach | Molecular Dynamics Computer Simulation of Amorphous Silica under High Pressure[END_REF] with liquid silica under isotropic compression. In the following, we chose to investigate the structural properties using the CHIK potential only.

Structural Properties

One of the main advantage of molecular dynamics simulation remains in the high level of detail on the microscopic structure. In silica glasses, the molecular structure presents a two-level organization. At short range (closest neighbors), oxygen atoms around silicon atoms tend to form tetrahedra at equilibrium. At medium range (two or three neighbors), the tetrahedra are linked all together to form rings with an average size of 6 tetrahedra per ring. The ring size will be in this work defined by the number of tetrahedra (or number of silicon atoms) in one ring. Additionally, a significant amount of free volume can also be observed inside a glass structure. Those cavities correspond to interstitial voids between rings created by defects (like dangling oxygen atoms for example). The tetrahedra, the ring statistic and the cavities distribution were analyzed to understand the overall behavior of silica glasses. In this section, the structure was extracted from 6 cases:

15

Page 15 of 52

ACS Paragon Plus Environment

The Journal of Physical Chemistry with the Hugoniostat method (6.7 GPa) and the equilibrium configuration at rest (labeled 0 GPa). The properties were computed behind the shock front to analyze the structure of shocked states. We have checked that the structure after the material relaxation was identical to the equilibrium configuration, except of course for the 8.7 GPa shock, which is slightly above the HEL. In this latter case, the final density is 2.336 g.cm -3 which is 6% higher than the initial density. The 8.7 GPa shock is shown only to illustrate the difference between the elastic and the densification regime. The mechanisms involved in the glass densification will be studied in a future work. Again, all results here were obtained with the CHIK potential.

The tetrahedral short-range structure was investigated by computing the radial distribution function (RDF) of the Si-O pair and the OSiO angle distribution, respectively in The RDF shows a narrow peak around 1.62 Å which is on average the experimental Si-O distance [START_REF] Chagarov | Application of Design of Experiments Methodology to Optimization of Classical Molecular Dynamics Generation of Amorphous SiO 2 Structure[END_REF] . The angle distribution has a maximum peak at 109 • which is typical of a tetrahedral organization. The maximum value seems to slightly decrease with the pressure and the distribution becomes wider. However this effect is rather weak and it can be hardly seen as a real structural change. The Hugoniostat method yields very similar results and the tetrahedral organization is very well reproduced.

Using the Si-O RDF, a geometrical criterion of 2 Å has been defined to investigate the coordination number of silicon atoms. The results are gathered in Table 3. The equilibrium configuration is mainly composed of 4-coordinated silicon atoms (99.6%) with a very low amount of defects (few number of 3 and 5-coordinated silicon atoms). The coordination does not change with pressure until 6.1 GPa. The 7.1 GPa and 8.7 GPa systems exhibit an increasing number of 5-coordinating silicon atoms (hexahedra) with pressure. Several octahedra are also observed in the 8.7 GPa simulation. It is well known that the presence of hexahedra is strongly correlated with a change in silica glass overall behavior 39 making 16

Page 16 of 52

ACS Paragon Plus Environment

The Journal of Physical Chemistry the material more ductile. It seems that an elastic compression below the HEL has a rather weak impact on the tetrahedral structure compared to a shock just above the HEL where a small amount of hexaedra and octahedra is observed. To understand the medium-range structure changes the ring organization was investigated. The algorithm developed by Franzblau was used to find the rings [START_REF] Franzblau | Computation of Ring Statistics for Network Models of Solids[END_REF] . A ring is identified when the path linking several silicon atoms is minimal (i.e. no short cut inside the ring) as defined by the short-path criterion [START_REF] Franzblau | Computation of Ring Statistics for Network Models of Solids[END_REF] . The rings organization is a characteristic feature of silica glasses. The 3-and 4-membered rings can also be observed by Raman spectroscopy as they are strongly correlated to two peaks, labeled D1 and D2, respectively 495 and 606 cm -1 41 . Many works including numerical simulations have stated that the silica glass ring distribution presents a maximum of 30% for 6-membered rings assuming the quenching rate is small enough [START_REF] Rahmani | Signature of Small Rings in the Raman Spectra of Normal and Compressed Amorphous Silica: A Combined Classical and ab initio Study[END_REF] . This is partly due to the fact that in crystalline phases (quartz, β-cristoballite), the structure is composed of 6-membered rings only. In amorphous phases, the equilibrium distribution is broader and ring size can range from 3-to 14-membered but the maximum for a size of 6 remains. The evolution of the ring distribution over pressure is shown in the Hugoniostat method and the NEMD simulations proving that this method is able to reproduce correct shocked states for glass materials.

Finally, the free volume was investigated by dividing the system in voxels and checking whether these voxels were empty or not. A voxel is declared empty if no particle is present within a radius of 2.1 Å from the voxel center of mass. This geometric criterion is the Van Der Waals radius of a silicon atom and was already used to probe cavities in silica glass [START_REF] Chan | Theoretical study of the Interstice Statistics of the Oxygen Sublattice in Vitreous SiO 2[END_REF][START_REF] Malavasi | Void Size Distribution in MDmodelled Silica Glass Structures[END_REF] .

The cavity volume is therefore defined as the empty space around silicon atoms which are approximated as hard spheres. In order to probe only the space between rings, a virtual particle was added inside each ring to model the ring center of mass and to ensure that no empty voxel could be found inside a ring. We were therefore able to distinguish, for each configurations, the volume occupied by the rings from the free volume. In order to analyze the ring shrinking and the cavity closing effects, the values were normalized here by the total volume of the material at rest, P=0 GPa. Results are gathered in Table 4. The total volume variation over pressure is of course related to the densities variations obtained in the previous section. Interestingly, 22% of the equilibrium configuration is occupied by cavities and 78% by the rings. This important fraction of empty volume inside silica glass is supposed to be the reason why the densification process occurs. Behind the shock front, the volume occupied by both rings and cavities decreases with the shock pressure. The volume variation is more important for the cavities. For example at 8.7 GPa, nearly half of the cavities were closed by the shock compression whereas the ring shrinking represents only 10% of the volume reduction. A shock compression induces two mechanisms : a ring shrinking and a cavity closing, with an effect more pronounced for the free volume closing mechanism.

After the material relaxation, the final volume occupation is identical to the equilibrium configuration in all elastic cases. In the system undergoing a 8. Several important conclusions may be drawn from the evolution of the glass structure in the elastic regime. First, the tetrahedral and the ring structure (i.e. bonds, coordination and angles) remains unaffected during the shock. Second, the interstitial voids present in the silica glass structure is strongly impacted by the elastic shock. The free volume is closing down with the increasing shock pressure. It is therefore concluded that the elastic behavior of silica glass is highly correlated with the interstitial voids initially present inside the structure. In addition, no significant differences were observed between the Hugoniostat method and the non-equilibrium approach. The two methods are equivalent and the final structures are identical within the pressure range tested here. The Hugoniostat method is able to reproduce the shocked states of glasses and is therefore an interesting alternative to NEMD simulations with a very low computational cost.

One simulation was also carried out just above the HEL to compare the elastic and densification regime. The tetrahedral structure is impacted and a significant amount of hexahedra and octahedra is observed. Regarding the rings distribution, the largest rings are broken into smaller ones. Above the HEL, the silica sample experiences deep structure modifications and the analysis of the relaxed states has shown a permanent densification after shock release related to a cavity closing mechanism. This phenomenon will be thoroughly studied in a future work.

3 Page

 3 used molecular dynamic simulations to investigate shock-induced ejection with high velocity impacts

4 Page

 4 presents the experimental setup. One laser shot consists in focusing, on the front face of the sample, the 1 mm in diameter spot of a laser

5 Page

 5

  in order to improve the performance of the BKS model for amorphous systems. This new CHIK model was adjusted to reproduce the radial distribution functions (RDF) of O-O, Si-O and Si-Si pairs obtained by Car-Parinello Molecular Dynamics (CPMD)

have shown that 6 Page

 6

10 Page

 10 the dynamic phenomenon involved in a shock propagation are not available through the Hugoniostat method. As we seek to investigate the structural properties of the shocked state of silica glass, we carried out simulations with the Hugoniostat method along with NEMD simulations. All the simulations in this work were performed with the MD code STAMP developed by the CEA.

13 Page

 13

Figure 8 .

 8 Figure 8.a and Figure 8.b. The configurations studied here presents very similar results.

Figure 9 . 17 Page 6 . 1 18 Page

 9176118 Figure 10.b. The RDF shows a progressive shift with increasing shock pressure. The magnitude of the second and third peaks, respectively observed around 4.8 Å and 7.2 Å, are identical, however the first peak magnitude tends to decrease slightly with pressure.The global structure does not change over pressure but the Si-Si distance is becoming smaller as the peaks are shifted toward lower distance values. This diminishing distance over pressure is linked to the ring collapsing mechanism. This phenomenon is also visible in the SiOSi angle distribution. As the shock pressure increases, the position of the maximum value is slightly shifted toward lower angles, from 145 • down to 130 • , while the magnitude remains the same. In the 8.7 GPa system, the second peak around 100 • is more pronounced than any other case. This particular angle is related to 3-membered rings which are increasing in number with the shock pressure as seen in Figure9.b. While the tetrahedral organization remains unaffected by the shock pressure, an overall ring closing is observed and this effect increases with the shock pressure as indicated by

Figure 1 :Figure 2 :

 12 Figure 1: Experimental setup of the GCLT laser and its instrumentations.

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Modified CHIK potential for pair interactions. Dotted lines represents the original Buckingham potential 19 diverging for low distances. In this work, these interactions were corrected with a two-order polynomial. Page 29 of 52

Figure 6 :

 6 Figure 6: Simulated particle velocity inside the material at different times. Velocities were extracted from NEMD simulations with the CHIK potential. a) Below 3 GPa (here 2.07 GPa), a ramp front is observed. The ramp propagates at a constant velocity (6272 m/s) and does not decay over time. b) Above 3GPa (here 7.01 GPa), a decaying shock front is obtained. The shock front becomes wider over time. At t=35 ps the front is nearly 300 nm wide.

Figure 7 :Figure 8 :Figure 9 :Figure 10 :

 78910 Figure 7: Comparison between experimental 1,5 and MD silica glass Hugoniot curves. a) Shock pressure as function of particle velocity. b) Shock pressure as function of density. Experimental (GCLT) error on shock pressure is around 10%. Error estimated in molecular simulations (not shown for clarity) is below 4%. Page 33 of 52

Table 2 : Shock states properties

 2

: Densities ρ, shock pressure P, particle velocities u p and shock front velocities u s .

Maximum error observed on shock front velocities is around 10% for GCLT experimental data and 4% in molecular simulations.

		ρ (g.cm -3 ) P (GPa) u p (m.s -1 ) u s (m.s -1 )
		2.200	0	0	0
	GCLT	2.274	1.81	164	5017
		2.276	2.38	190	5694
		2.282	2.63	207	5775
		2.288	3.19	236	6144
		2.315	3.47	275	5736
		2.334	3.87	318	5532
		2.350	4.15	348	5421
		2.346	4.74	366	5887
		2.356	5.61	411	6204
		2.203	0	0	0
	Hugoniostat	2.300	3.05	244	5918
	CHIK	2.400	5.63	460	6016
		2.500	6.71	603	5654
		2.203	0	0	0
	NEMD	2.249	1.35	100	6228
	CHIK	2.272	2.07	150	6272
		2.290	2.73	200	6396
		2.332	3.94	300	6262
		2.380	5.09	400	6176
		2.433	6.12	500	6056
		2.500	7.07	600	5949
		2.675	8.71	800	5742
		2.260	0	0	0
	NEMD	2.280	1.39	100	6250
	BKS	2.329	2.62	200	5996
		2.380	3.81	300	5919
		2.449	4.89	400	5809
		2.533	5.56	500	5420
		2.686	6.23	600	5194
			14

Page 14 of 52 ACS Paragon Plus Environment The Journal of Physical Chemistry tions

  . Experimental space and time scales are far beyond the scales of molecular dynamics simulations. Discrepancies between numerical and experimental works is not a surprise especially on complex dynamic processes. In addition, molecular dynamics simulations relies heavily on the potential. Complex reactive potentials may provide more accurate results, however, such work were never carried out to our knowledge.In order to sum up, considering the experimental error estimated around 10%, the pressures and densities provided by molecular simulations are in a very good agreement with experiments. Despite the discrepancy on the shock front velocities, NEMD simulations and the Hugoniostat method are able to reproduce the shocked states of silica in the elastic regime and to some extent the non-linear behavior (ramp and decaying waves) of glassy systems. The CHIK model especially gives a fair description of the amorphous silica behavior with pressure, as already observed by Horbach

Table 3 : Evolution of silicon atoms coordination number (number of oxygen atoms around every silicon atoms) observed in the shocked state just behind the shock front.

 3

			Coordination Number
	Pressure (GPa)	3	4	5	6
	0.0 GPa	0.14% 99.6% 0.26% 0.00%
	2.0 GPa	0.12% 99.5% 0.38% 0.00%
	2.7 GPa	0.04% 99.5% 0.46% 0.00%
	3.9 GPa	0.17% 99.3% 0.53% 0.00%
	5.1 GPa	0.04% 99.0% 0.96% 0.00%
	6.1 GPa	0.04% 98.5% 1.46% 0.00%
	7.1 GPa	1.78% 95.9% 2.32% 0.00%
	8.7 GPa	0.01% 92.5% 7.39% 0.06%

Page 19 of 52 ACS Paragon Plus Environment The Journal of Physical Chemistry entirely

  due to a cavity closing. The rings occupy the same volume as before the shock. It seems the densification process is strongly correlated with a permanently cavity closing mechanism inside the structure. Higher shock pressures will be studied to analyze in detail the densification in a future work.

7 GPa shock loading, a residual densification is observed. The final volume after relaxation is 6% smaller and is 19

Table 4 : Rings and cavities volumes inside the glass structure observed in the shocked state just behind the shock front.

 4

		Volume Ratio
	Pressure (GPa) Rings Cavities Total Volume
	0.0 GPa	0.7801 0.2199	1.0000
	2.0 GPa	0.7792 0.2119	0.9911
	2.7 GPa	0.7718 0.2020	0.9738
	3.9 GPa	0.7629 0.1929	0.9558
	5.1 GPa	0.7562 0.1804	0.9366
	6.1 GPa	0.7483 0.1678	0.9161
	7.1 GPa	0.7376 0.1544	0.8920
	8.7 GPa	0.7108 0.1225	0.8333
	Conclusion
	In this work, experimental and numerical studies were combined to improve the under-
	standing of the silica glass behavior under elastic shock compression. The elastic regime
	was investigated by laser-driven shocks up to 5 GPa and by molecular dynamics simula-
	tions up to 8.7 GPa. Direct shocks induced by a mobile wall and a Hugoniostat method
	were used to reach silica shocked states. In our simulations, two pair potentials, the BKS
	and the CHIK models, were tested. Fused silica exhibits nonlinear properties and NEMD
	simulations are able to reproduce the particular behavior (ramp waves and decaying
	wave propagation) under elastic shocks. The Hugoniot curves obtained with the CHIK
	potential were in a very good agreement with the experimental curves. In the opposite,

the BKS potential fails to reproduce the shocked states of amorphous silica. Surprisingly, the CHIK potential is able to give a very good description of the shocked states in spite of its simplicity.

20 Page 20 of 52 ACS Paragon Plus Environment The Journal of Physical Chemistry

ACS Paragon Plus EnvironmentThe Journal of Physical Chemistry

Acknowledgment

The authors are grateful to the French Agence Nationale de la Recherche for its financial support through the program GLASS.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

ACS Paragon Plus Environment

The Journal of Physical Chemistry

ACS Paragon Plus Environment

The Journal of Physical Chemistry


Jean-Claude Bajard
email: jean-claude.bajard@lip6.fr

Julien Eynard
email: jeynard@uwaterloo.ca

Anwar Hasan
email: ahasan@uwaterloo.ca

Paulo Martins

Leonel Sousa

Vincent Zucca
email: vincent.zucca@lip6.fr

Efficient reductions in cyclotomic rings -Application to Ring-LWE based FHE schemes

Keywords: Polynomial Reduction, Number Theoretic Transform, Residue Number Systems, Ring-Learning With Errors, Homomorphic Encryption

With Fully Homomorphic Encryption (FHE), it is possible to process encrypted data without access to the private-key. This has a wide range of applications, most notably the offloading of sensitive data processing. Most research on FHE has focused on the improvement of its efficiency, namely by introducing of schemes based on Ring-Learning With Errors (RLWE), and techniques such as batching, which allows for the encryption of multiple messages in the same ciphertext. Much of the related research has focused on RLWE relying on power-of-two cyclotomic polynomials. While it is possible to achieve efficient arithmetic with such polynomials, one cannot exploit batching. Herein, the efficiency of ring arithmetic underpinned by non-power-of-two cyclomotic polynomials is analyzed and improved. Two methods for polynomial reduction are proposed, one based on the Barrett reduction and the other on a Montgomery representation. Speed-ups of up to 2.66 are obtained for the reduction operation using an i7-5960X processor when compared with a straightforward implementation of the Barrett reduction. Moreover, the proposed methods are exploited to enhance homomorphic multiplication of Fan-Vercauteren (FV) and Brakerski-Gentry-Vaikuntantahan (BGV) encryption schemes, producing experimental speed-ups of up to 1.37.

Introduction

There is an increasing conflict between the convenience provided by cloud services and privacy concerns. Privacy can be achieved by encrypting data before uploading it to the cloud. However, with traditional cryptosystems, one is not able to process encrypted data [START_REF] Goluch | The development of homomorphic cryptography[END_REF], nullifying the benefits of cloud computing. A solution to this problem is the application of FHE, which allows for the creation of malleable cryptograms [START_REF] Goluch | The development of homomorphic cryptography[END_REF]. Homomorphic operations can afterwards be applied to these cryptograms. Despite its wide range of applicability, FHE has seldom been applied in practice due to its low computational performance.

Most recent research on FHE has focused on improving its efficiency [START_REF] Martins | Enhancing Data Parallelism of Fully Homomorphic Encryption[END_REF]. LWE [START_REF] Regev | On Lattices, Learning with Errors, Random Linear Codes, and Cryptography[END_REF] and its ring-variant RLWE [START_REF] Lyubashevsky | On Ideal Lattices and Learning with Errors over Rings[END_REF] has been suggested in this context as a framework for improving the complexity of FHE. The benefits brought forth by RLWE are twofold. First, operations are executed in a cyclotomic ring, and therefore benefit from its algebraic structure. Second, the plaintext space is isomorphic to the Cartesian product of multiple smaller spaces, designated batching slots, allowing for multiple messages to be encrypted and processed in a single ciphertext in parallel. Most of the related research has focused on cyclotomic rings of the form Z[X]/(X N + 1), for N a power of two, due to the simple arithmetic associated with the resulting ring. However, they only allow for a single batching slot. While RLWE modulo other cyclotomic polynomials, providing for a large number of batching slots has been previously considered [START_REF] Gentry | Homomorphic Evaluation of the AES Circuit[END_REF][START_REF] Smart | Fully homomorphic simd operations[END_REF], it remains not used much because of its less efficient arithmetic. In [START_REF] Lyubashevsky | On Ideal Lattices and Learning with Errors over Rings[END_REF], Lyubashevsky et al. proposed algorithms using a multivariate/tensored representation of the elements of cyclotomic rings resulting in a more efficient arithmetic in the general case than those using the classical univariate representation. The implementations of their work [START_REF] Crockett | λ • λ: Functional lattice cryptography[END_REF][START_REF] Mayer | Implementing a toolkit for ring-lwe based cryptography in arbitrary cyclotomic number fields[END_REF] are, at the best of our knowledges, the only publicly available implementations optimized for general cyclotomic rings.

In this work we analyse and improve the arithmetic associated with the univariate representation in general cyclotomic rings with efficient reduction algorithms of same asymptotic complexities than those of [START_REF] Lyubashevsky | On Ideal Lattices and Learning with Errors over Rings[END_REF]. Barrett reduction, which is a generic technique to perform polynomial reduction, had been considered in [START_REF] Dai | cuHE: A Homomorphic Encryption Accelerator Library[END_REF] to perform such reductions. However, the algorithms used in [START_REF] Dai | cuHE: A Homomorphic Encryption Accelerator Library[END_REF] do not take into account the characteristics of cyclotomic polynomials, and hence are sub-optimal. Herein, the degree of the polynomials is reduced using cyclotomic properties before applying Barrett's algorithm, decreasing the complexity of the reduction and leading to theoretical speed-ups of up to 2.06. We also show that a Montgomery representation leads to more efficient reductions than generic Barrett algorithms, leading to theoretical speed-ups of up to 2.63.

These gains have been confirmed in practice, with experimental speed-ups of up to 1.95 and 2.55 for the improved Barrett and Montgomery reductions, respectively, when compared with a straightforward Barrett algorithm in an i7-5960X processor. Moreover we have tested the applicability of our algorithms to two of the most currently used homomorphic encryption schemes, BGV [START_REF] Brakerski | Fully homomorphic encryption without bootstrapping[END_REF] developed in HELib [START_REF] Halevi | Algorithms in helib[END_REF] and FV [START_REF] Fan | Somewhat Practical Fully Homomorphic Encryption[END_REF] developped in SEAL 2.0 [START_REF] Laine | Simple encrypted arithmetic library -seal (v2.0)[END_REF] resulting in speed-up of up to 1.37 for homomorphic multiplication. These results should reduce the performance differences between the univiate/multivariate representations and therefore bring more flexibility to future implementations.

Background

Throughout the paper, φ m (X ) ∈ Z[X] will denote the m-th cyclotomic polynomial of degree n = ϕ(m), where ϕ is Euler's totient function. The ring R = Z[X]/(φ m (X )) is the main structure of R-LWE-based schemes such as FV and BGV. An element of R can be thought of as a polynomial with integer coefficients and a degree strictly smaller than n. Unless mentioned otherwise, polynomials are represented in the power-basis {1, X, . . .

X n-1 }. For a = n-1 i=0 a i X i ∈ Z[X],
we denote a ∞ = max{|a i |, 0 i < n}. Finally, the underlying space for ciphertexts is R q = R/qR = Z q [X]/φ m (X ), which is composed of elements of R with coefficients reduced modulo q. The notation | • | q is used to denote the classical residue modulo q in [0, q), while the centered residue in [-q/2, q/2) is denoted by [•] q . Moreover • denotes flooring while • denotes rounding to the nearest integer.

Residue Number System (RNS)

In practice, the value of q underpinning R q is chosen to be the product of small different prime numbers q = q 1 • • • q k . Therefore, thanks to the Chinese Remainder Theorem (CRT), R q splits into a cartesian product of smaller rings through the following isomorphism:

RNS q=q 1 ...q k : R q → R q 1 × . . . × R q k a → (a mod q 1 , . . . , a mod q k ) (1)
The RNS exploits (1) to transfer the arithmetic modulo q of the coefficients to k smaller arithmetics modulo each q i . Thus, better performance is achieved due to smaller arithmetics and parallel computations over each R q i .

Product of elements in R q

Since the degree n of the polynomials is large in practice, the polynomial product is a major bottleneck for the efficiency of R-LWE based schemes. However, when n is close or equal to a power of two, it can be performed efficiently thanks to the Number Theoretic Transform (NTT). Let N 2 be the function defined over N such that for any n ∈ N N 2 (n) is the smallest power of two greater than or equal to n. In our context the products have a degree strictly smaller than 2n, therefore we denote N = N 2 (2n). For a prime q, a primitive N th -root of unity ω ∈ Z q exists if and only if q ≡ 1 mod N. For a Z q equipped with ω, the following ring isomorphism is verified:

NTT q, N,ω : Z q [X]/(X N -1) → Z q [X]/(X -1) × • • • × Z q [X]/(X -ω N -1 ) a → a = a(1), a(ω), . . . , a(ω N -1 ) (2)
Once ( 2) is applied to obtain the NTT representation of polynomials, their product can be performed coordinate-wise. The time-complexity of arithmetic becomes linear in N and the bottleneck changes to the evaluation of (2), as well as of its inverse whose evaluation require O(N log(N )) multiplications in Z q by state-of-the-art algorithms ( [START_REF] Harvey | Faster arithmetic for number-theoretic transforms[END_REF], [START_REF] Longa | Speeding up the Number Theoretic Transform for Faster Ideal Lattice-Based Cryptography[END_REF]). When the context is clear, we will denote an NTT transformation of degree N by NTT N instead of NTT q, N,ω .

In order to efficiently compute the product of elements a and b of R q , seen as polynomials over Z q [X] of degree at most n -1, one has to compute the NTT representation of c = a × b of degree 2n -2 through:

NTT N (c) = NTT N (a) NTT N (b) (3)
where denotes the component-wise multiplication in Z q . To obtain the value of c = a × b ∈ R q , a second step is needed, which consists of reducing the result of (3) modulo φ m (X ). Notice that when m is a power of two (φ m (X ) = X m/2 + 1), one can use a negatively-wrapped convolution for the evaluation of the NTT and its inverse [START_REF] Longa | Speeding up the Number Theoretic Transform for Faster Ideal Lattice-Based Cryptography[END_REF]. This allows us to use an NTT of size N 2 (n) instead of N 2 (2n) for the evaluation of (3) and also to recover the polynomial reduced modulo X m/2 + 1 at the end of (3) just with an inverse NTT of size N 2 (n). However, for other values of m this method cannot be applied and a more complex reduction has to be carried out after applying (3). Barrett's strategy for modular reduction over the integers [START_REF] Barrett | Implementing the Rivest Shamir and Adleman Public Key Encryption Algorithm on a Standard Digital Signal Processor[END_REF] can be adapted to polynomial modular arithmetic to reduce a polynomial c of degree n+α by φ m of degree n. The quotient of the Euclidean division c/φ m is computed through multiplications by precomputed constants and shifts. We present this strategy in Algorithm 1. The performance of the algorithm is directly related to the size of the polynomial to be reduced: the algorithm is more efficient when α is small. By denoting ñ = N 2 (n) and A = N 2 (2α + 1) the cost of the algorithm is:

C NTT (N ) + 2C NTT (ñ) + 2C NTT ( A) + (ñ + A)MMult q
where C NTT (x) denotes the cost for evaluating (2) (or its inverse) of size x and MMult q the cost of a modular multiplication modulo q. One may also notice that Barrett's reduction used in [START_REF] Dai | cuHE: A Homomorphic Encryption Accelerator Library[END_REF] uses NTT of size N = N 2 (2n) to perform the second product while in fact it only requires NTT of size N 2 (n). For the sake of completeness a proof of the correctness of Alg. 1 is given in Appendix A.1.

Algorithm 1 NTTBarr(P, Z q , φ m ): NTT based Barrett reduction in Z q [X], for a prime q ≡ 1 mod N Require:

c NTT = NTT N (c) ∈ Z N q with deg(c) = n + α < 2n with q prime, n = deg(φ m ), ñ = N 2 (n), A = N 2 (2α + 1); precomputed NTT ñ (φ m ) and NTT A ( X n+α /φ m ). Ensure: c mod (q, φ m ) in power-basis. 1: c ← NTT -1 N (c NTT ) 2: f ← c/X n 3: r ← NTT -1 A (NTT A ( f ) NTT A ( X n+α /φ m )) 4: r ← r/X α 5: d ← NTT -1 ñ (NTT ñ (r ) NTT ñ (φ m )) 6: c ← c mod X ñ -1 7: return c -d 2.

RNS variant of the FV and BGV Encryption Schemes

Fan and Vercauteren [START_REF] Fan | Somewhat Practical Fully Homomorphic Encryption[END_REF] have adapted Brakerski's scale invariant FHE scheme [START_REF] Brakerski | Fully homomorphic encryption without modulus switching from classical GapSVP[END_REF] to the RLWE framework. More recently, Bajard et al. have provided a full RNS variant of FV [START_REF] Bajard | A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes[END_REF]. We briefly recall how this RNS variant works.

We first need to introduce the two functions ξ q : R q → R q 1 × • • • × R q k and P RNS,q : R q → R k q such that for any a ∈ R q , ξ q (a) = (| a • q i /q| q i ) i ∈[1,k] and P RNS,q (a) = (| a • q/q i | q ) i ∈ [1,k] . It is straightforward to notice that for any (a, b) ∈ R 2 q , ξ q (a), P RNS,q (b) ≡ ab mod q. This scalar product occurs in R q , and so it is composed of polynomial products.

In the context of the FV scheme, the secret-key s ∈ R is defined as a "small" polynomial drawn from a distribution χ key . An encryption of m ∈ R t corresponds to a pair of polynomials ct = (c 0 , c 1 ) ∈ R 2 satisfying:

[c 0 + c 1 s] q = [ q/t • [m] t + v] q (4)
where v is a noise term that is originally introduced during encryption (which is related to a distribution χ err ) and that grows as homomorphic operations are applied. Decryption works correctly so long as this noise is below a certain bound, which limits the amount of operations one can perform. The homomorphic addition of two ciphertexts corresponds to the pairwise addition of the ciphertexts' polynomials. Regarding homomorphic multiplication, it is useful to see ciphertexts as polynomials of degree 1 with coefficients in R. In this context, homomorphic multiplication takes place in two steps. First,

ct mult ← t q c 1 0 c 2 0 q , t q (c 1 0 c 2 1 + c 1 1 c 2 0 ) q , t q c 1 1 c 2 1 q
is computed with a Karatsuba like algorithm. During this procedure, the RNS representations of the input polynomials are extended to bases with larger dynamic ranges so as to compute the products over R instead of over R q . Moreover, the division operation is achievable using [START_REF] Bajard | A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes[END_REF]Section 4.4]. Finally, one has to convert the three-element ciphertext back to a two-element ciphertext, through a process called relinearization. This is done by multiplying ξ q (ct mult ) by pseudo-encryptions of P RNS,q (s 2 ) (designated --→ rlk), and adding the result to the other two elements:

ct r elin ← c 0 + ξ q (c 2 ), --→ rlk 0 q , c 1 + ξ q (c 2 ), --→ rlk 1 q ∈ R 2 q . (5)
The scheme introduced by Brakerski, Gentry and Vainkuntanathan [START_REF] Brakerski | Fully homomorphic encryption without bootstrapping[END_REF] shares many features of FV, and can be similarly adapted to the techniques in [START_REF] Bajard | A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes[END_REF]. A secret-key is also defined to be a "small" polynomial s ∈ R q , and ciphertexts correspond to pairs (c 0 , c 1 ) ∈ R 2 , but messages are encrypted in the Least Significant Bits (LSBs) of ( 6):

[c 0 + c 1 s] q = [[m] t + tv] q (6)
The change in the positioning of the message bits leads to simpler homomorphic multiplications. First, we compute the degree 2 ciphertext ct mult ←

c 1 0 c 2 0 q , (c 1 0 c 2 1 + c 1 1 c 2 0 ) q , c 1 1 c 2 1 q
. Since operations are performed modulo q, no RNS base extension is required. Afterwards, an operation similar to ( 5) is applied, so as to convert the three-element ciphertext to a classical two-element ciphertext. Finally, a noise management technique is applied to reduce the growth rate of the norm of v in [START_REF] Joppew | Improved Security for a Ring-Based Fully Homomorphic Encryption Scheme[END_REF]. This technique consists of scaling the ciphertext to a smaller ring R q with an appropriate rounding, and is performed in two steps:

δ i ← t • [-ct/t] q/q for i = 0, 1 ct ← ([q /q • (c 0 + δ 0 )] q , [q /q • (c 1 + δ 1 )] q )
In certain steps of the aforementioned schemes, one needs to add the result of multiple polynomial products. In this case, it is possible to combine the modular reduction of all products involved in the same sum, reducing the overall computational complexity. First, ( 3) is computed for all products. Then, the NTT representations of the products are added. Finally, a single polynomial reduction method can be applied to the result.

Batching

A common way to improve the efficiency of RLWE schemes is to encrypt several plaintexts in a single ciphertext, through a technique called batching [START_REF] Smart | Fully homomorphic simd operations[END_REF]. Under some conditions, φ m splits modulo t into distinct irreducible polynomials f 1 , . . . , f of degree n/ . This leads to the following ring isomorphism of the plaintext space

R t : R t Z t [X]/( f 1 ) ו • •×Z t [X]/( f ).
In this manner, plaintexts m 1 , . . . , m can be compactly represented as a single polynomial m ∈ R t . Afterwards m is encrypted and homomorphic operations applied to this ciphertext operate on each slot individually. This technique is mostly used when evaluating Boolean circuits, i.e. with t = 2, to pack -bits in a single ciphertext. However, since X m/2 + 1 ≡ (X + 1) m/2 mod 2, this technique cannot be used when m is a power of two. Thus the efficient arithmetic associated with power-of-two cyclotomic polynomials has limited applicability.

Improving polynomial reduction modulo φ m

In this section, we propose two efficient methods to compute polynomial reductions. The first method takes advantage of the properties of the cyclotomic polynomials to improve the efficiency of the Barrett algorithm. The second reduction rests on an adaptation of the Montgomery modular reduction [START_REF] Peter | Modular multiplication without trial division[END_REF] and can be adapted for other polynomials besides cyclotomics.

Improving Barrett's reduction for cyclotomic polynomials

As explained in section 2, the performance of Barrett's algorithm is sensitive to the difference between the degree of the polynomial to be reduced and that of the polynomial we want to reduce by. The smaller the difference, the more efficient the algorithm will be. Hence, herein we propose an efficient method to reduce this difference.

Polynomials to be reduced modulo φ m have a degree of at most 2n -2. Let c be such a polynomial. If c were reduced by a polynomial Q sp of degree n + α + 1, the difference between the degree of the polynomial and the degree of φ m would be reduced to α. However, in order to obtain the correct value of c mod φ m in the end, φ m has to divide Q sp . For this reduction to be efficient, Q sp should be sparse enough so that its reduction can be efficiently handled through few operations in Z q . Thanks to the cyclotomic property d/m φ d (X ) = X m -1, Q sp can be taken as the product of φ m and some φ d for d dividing m. Good candidates can be found by recursively using the fact that if p is a prime not dividing m then φ m •p (X ) • φ m (X ) = φ m (X p ). If the polynomial Q sp is found in this way with less than 2 recursions it has coefficients in {-1, 0, 1} since it will correspond to a cyclotomic with at most two distinct odd prime factors. In this case, the reduction modulo Q sp only requires shifts and additions in Z q and can be done very efficiently.

In addition, when m < 2n -2, c can initially be reduced by X m -1 with 2nm -1 additions in Z q . Since φ m (X )|Q sp (X )|X m -1 the strategy remains correct, and the complexity of the reduction by Q sp (X ) is further reduced. Let HW(Q sp ) be the Hamming weight of Q sp . The cost of the reduction of c by

Q sp is (HW(Q sp ) -1)(m -deg(Q sp )) additions in Z q . At this point, we obtain c = c mod Q sp (with deg(c ) ≤ n + α) and c ≡ c mod φ m .
The final algorithm is depicted in Algorithm 2. It starts by recovering c in power basis from the NTT representation outputted by [START_REF] Bajard | A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes[END_REF]. Then it consecutively reduces c of degree 2n-2 by X m -1 and the sparse polynomial Q sp . This allows to recover c = c mod Q sp of degree n + α very efficiently. Afterwards, a procedure identical to the one described in steps 2 to 7 in Alg. 1 is applied to c to get c = c mod φ m .

Algorithm 2 modBt φ m : NTT-based Barrett reduction in Z q [X], for q = q 1 . . . q k , with prime integers q i , q

i ≡ 1 mod N, N = N 2 (2n), ñ = N 2 (n) and A = N 2 (2α+1). Require: c NTT = NTT N (c) with deg(c) 2n -2 Ensure: c = c mod φ m in power-basis. 1: c ← NTT -1 N (c NTT ) 2: if m < 2n -2 then 3: c ← c mod X m -1 4: c ← c mod Q sp Reduction by Q sp of degree n + α + 1 5: f ← c/X n 6: r ← NTT -1 A (NTT A ( f ) NTT A ( X n+α /φ m )) 7: r ← r/X α 8: d ← NTT -1 ñ (NTT ñ (r ) NTT ñ (φ m )) 9: c ← c mod X ñ -1 10: return c -d
The impact of this sparse reduction is illustrated in Table 1, where polynomials Q sp are presented for different cyclotomic polynomials. Cyclotomics have been chosen with a degree n = ϕ(m) close or equal to a power of two. The number of batching slots associated with each cyclotomic is also presented. The degree of Q sp is n + α + 1 thus NTTs of size A = N 2 (2α + 1) are required to compute the first polynomial product in Alg. 2. This is in contrast with N = N 2 (2n) which would have been the size required for the Barrett algorithm without using the sparse reduction. In order to highlight the sparsity of Q sp we give HW(Q sp ) which is the number of non-zero coefficients of Q sp .

Complexity Since the complexity of computing multiplications in Z q is much higher than additions, the cost of the reduction by the sparse polynomial can be neglected. Moreover, with the RNS, each multiplication in R q , with q = q 1 . . . q k can be decomposed into k independent and smaller multiplications. Therefore the cost to reduce the polynomial c outputted by ( 3) is essentially k times the cost of Alg. 1:

k (N log 2 (N ) + 2A log 2 ( A) + 2ñ log 2 (ñ) + A + ñ)
While the cost of the method by using directly Barrett's algorithm, i.e. without performing the reduction by the sparse polynomial, is:

k (3N log 2 (N ) + 2ñ log 2 (ñ) + N + ñ)
Based on this analysis we also provide in Table 1 the theoretical speed-up obtained with the use of the sparse reduction. 1. Sparse polynomials used for partial reduction with their related parameters.

m n Q sp deg(Q sp ) α HW(Q sp ) A N Speed-

NTT-based Montgomery's reduction

We propose a Montgomery like reduction of a polynomial given in NTT representation inspired by [START_REF] Bajard | Arithmetic Operations in Finite Fields of Medium Prime Characteristic Using the Lagrange Representation[END_REF]. The bottleneck of our optimized Barrett algorithm is the computation of the inverse NTT of size N of (3). Our Montgomery reduction takes advantage of the presence of the NTT basis of size N/2 (seen as an RNS basis in [START_REF] Bajard | Arithmetic Operations in Finite Fields of Medium Prime Characteristic Using the Lagrange Representation[END_REF]) in the basis of size N allowing then to perform all the computations, in particular the inverse NTT evaluation, in the basis of size N/2 instead of N.

The NTT representation of a polynomial of size N was defined in (2) as the tuple of values {c mod X -ω j |0 ≤ j < N }. This representation can be seen as a polynomial-RNS representation of c mod X N -1 since X N -1 = 0≤ j<N (Xω j ) mod q, with respect to the NTT-basis:

B ω, N = {|X -1| q , |X -ω| q , . . . , |X -ω N -1 | q }
As X N -1 splits in (X N/2 -1)(X N/2 + 1) when N is even, half of the NTT N representation of c corresponds to its NTT N/2 representation. Hence, the basis B ω, N is split along even and odd powers of ω. We can then define two sub-bases defining two polynomials:

             B (e) ω, N = {[X -ω 2j | q , 1 ≤ j ≤ N/2} ; Ψ (e) = | N/2 j=1 (X -ω 2j )| q B (o) ω, N = {|X -ω 2j+1 | q , 0 ≤ j ≤ N/2 -1} ; Ψ (o) = | N/2-1 j=0 (X -ω 2j+1 )| q (7)
It is straightforward to notice that Ψ (e) ≡ X N/2 -1 mod q and Ψ (o) ≡ X N/2 + 1 mod q. We also note that since N is a power of two, one has Ψ (o) ≡ φ N mod q. Thanks to Lemma 1, whose first point is a direct consequence of Lemma 2 in [START_REF] Filaseta | On coverings of the integers associated with an irreducibility theorem of A. Schinzel. Number theory for the millennium[END_REF], we can choose X N/2 + 1 as the Montgomery factor.

Lemma 1. Let φ m be the m-th cyclotomic polynomial of degree n and N be the smallest power of two greater than or equal to 2n. If m is not a power of two then:

there exists

(U, V ) ∈ Z[X] 2 such that U (X )φ m (X ) + V (X )φ N (X ) = 1.
for any prime p, φ m is coprime with X N -1 in Z p . In particular φ m is a unit in Z p [X]/(φ N ).

One can extract from the coordinates of c in B ω, N the representation c (e) of c in B (e) ω, N (resp. c (o) in B (o) ω, N ). So, given c (o) and c (e) , we can use the NTT operator to get:

NTT -1 N/2 ( c (e) ) = c mod (q, X N/2 -1) Definition 1. We define the following function, which takes in as input the residues of the polynomial c (deg(c) < N) modulo a prime p:

modMong φ m ,Ψ (o) , p (c) = c + φ m × | -c/φ m | Ψ (o) Ψ (o) mod p. (8)
The modMong φ m ,Ψ (o) , p function defined in ( 8) is a classical Montgomery reduction with factor Ψ (o) consisting in an exact polynomial division. It always outputs a polynomial congruent to |c/Ψ

(o) | φ m but when deg(c) N/2 + n -1 the output is exactly |c/Ψ (o) | φ m . Lemma 2. If deg(c) N/2+n-1, then modMg φ m ,Ψ (o) , p ( c) = c/Ψ (o) mod (p, φ m ).
First the degree of the numerator in ( 8) is bounded max(deg(c), deg(φ m ) + deg(Ψ (o) ) -1) n+ N/2-1. Thus, the degree of the resulting quotient is bounded by n -1 < N/2. Therefore, the output is |c/Ψ (o) | φ m and the computation of ( 8) can be made modulo X N/2 -1, i.e. in an NTT representation of size N/2 when using primes q i ≡ 1 mod N. Algorithm 3 details the computation of [START_REF] Brakerski | Fully homomorphic encryption without bootstrapping[END_REF]. The following precomputations are used therein:

         W (o) : (-1/φ m ) mod (q, Ψ (o) ) in base B (o) ω, N Y (e) : (1/Ψ (o) ) mod (q, Ψ (e) ) in base B (e) ω, N Z (e) : (φ m /Ψ (o) ) mod (q, Ψ (e) ) in base B (e)
ω, N Algorithm 3 modMg φ m ,Ψ (o) : NTT-based Montgomery reduction in Z q [X], for q = q 1 . . . q k , with prime integers q i , and q i ≡ 1 mod N 2 (2n) BaseConv : (c(ω), c(ω 3 ), . . . , c(ω N -1 )) → (c(1), c(ω 2 ), . . . , c(ω N -2 )) (modq).

Require: c = NTT q, N (c) (ie c in base B (o) ω, N ∪ B (e) ω, N ), with N = N 2 (2n) and deg(c) 2n -2 < N/2 + n -1. Ensure: R = (c/Ψ (o) ) mod (q, φ m ) in power-basis. 1: ( c (e) , c (o) ) ← Split( c) Split the NTT coeff. wrt parity of indexes 2: Q (o) ← c (o) W (o) 3: Q (e) ← BaseConv( Q (o) ) base conversion from B (o) ω, N to B (e)
(9) In [START_REF] Bajard | Arithmetic Operations in Finite Fields of Medium Prime Characteristic Using the Lagrange Representation[END_REF], ( 9) is computed with a classical Lagrange interpolation. Our context is more specific, because the points in which polynomials are evaluated are powers of a N th root of unity ω. With this purpose, Alg. 4 implements such base conversion by only using NTTs of degree N/2, with ω 2 as a primitive N/2 th root of unity. The proof of correctness of Alg. 4 is given in Appendix A.3.

Complexity The total cost of Alg. 3 in terms of modular multiplications is:

k 3N 2 (n) log 2 (N 2 (n)) + 4N 2 (n)
One can find in Table 2 the predicted speed-up of the proposed Montgomery reduction. Despite its better complexity, the Montgomery algorithm suffers from Algorithm 4 BaseConv Require: (c(ω), c(ω 3 ), . . . , c(ω N -1 )) mod q i , for deg(c) < N/2, N a power of 2 and ω a primitive N th root of unity in Z q . Ensure: (c(1), c(ω 2 ), . . . , c(ω N -2 )) mod q.

1:

(c 0 , c 1 , . . . , c N/2-1 ) ← NTT -1 N/2,ω 2 (c(ω), c(ω 3 ), . . . , c(ω N -1 )) 2: (p 0 , p 1 , . . . , p N/2-1 ) ← (c 0 , c 1 , . . . , c N/2-1 ) (1, ω -1 , ω -2 , . . . , ω -(N/2-1) ) 3: (r 0 , r 1 , . . . , r N/2-1 ) ← NTT N/2,ω 2 (p 0 , p 1 , . . . , p N/2-1 ) 4: return (r 0 , r 1 , . . . , r N/2-1 )
one main drawback which is the presence of the Montgomery factor in the output. In the following section, modifications to the BGV and FV cryptosystems are proposed to handle this factor.

Adaptation of FV and BGV to the Montgomery representation

In this section we show how the Montgomery representation impacts the BGV and FV schemes, and propose modifications to handle these changes. For simplicity, we denote by M the Montgomery factor X N/2 + 1 mod φ m . Thanks to Lemma 1 we also know that M -1 exists in R. We assume that ciphertexts ct are given in Montgomery representation such that ct = (c 0 M, c 1 M). The conversion to the Montgomery domain can be integrated in the encryption procedure for increased efficiency, and the M factor can be removed during decryption by applying a Montgomery reduction to [c 0 + c1 s] q . The Montgomery reduction only impacts procedures involving multiplications in R q by multiplying the product by M -1 . Therefore the Montgomery representation is stable with respect to multiplication. Homomorphic additions are not affected by the change in representation. Thus the only impact one has to consider is on homomorphic multiplication. Moreover, we recall that the expansion factor of the ring R is the quantity defined by

δ R = sup{ ab ∞ / a ∞ b ∞ (a, b) ∈ R -{0}}.

Impact of the Montgomery representation in FV

We note that the first step of the FV homomorphic multiplication corresponds to an extension of the polynomials of the ciphertexts to a larger RNS basis, so that multiplications are computed over R instead of R q . In order to improve efficiency, an approximate extension is used [START_REF] Bajard | A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes[END_REF] and thus the norm of the polynomials is bounded by q 2 (1 + ρ) for a parameter ρ > 0 [START_REF] Bajard | A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes[END_REF]. A bound on the noise growth of ct mult ←

t q M c 1 0 c 2 0 q , t q M (c 1 0 c 2 1 + c 1 1 c 2 0 ) q , t q M c 1 1 c 2 1 q
is given in Proposition 1 whose proof can directly be derivated from the one of [START_REF] Bajard | A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes[END_REF].

Proposition 1. Let (c 0 , c1 , c2 ) = ct mult , r ∞ = 1+ρ 2 (1+δ R s ∞ ) +δ R M ∞ and v i be the inherent noise of ct i . Then (c 0 +c 1 s+c 2 s 2 )M -1 ≡ ∆ [m 1 m 2 ] t +ṽ mult [q] with: ṽmult ∞ < δ R t(δ R M -1 ∞ r ∞ + 1 2 )( v 1 ∞ + v 2 ∞ ) + δ 2 min v i ∞ + δt|q| t (r ∞ + 1) + δ R M -1 ∞ 2 (1 + δ R s ∞ (1 + δ R s ∞ )) + |q | t 2 + 1. (10)
Now, we assume that we need to relinearize the ciphertext ct mult = ( c 0 , c 1 , c 2 ). Before, the following dot products were computed over R q : ξ q (c 2 ), rlk i , where

evk 0 = [P R N S,q (s 2 ) + - → a s + - → e ] q and evk 1 = [-- → a ] q .
The goal of relinearisation is to obtain ξ q (c 2 ), P R N S,q (s 2 ) ≡ c 2 s 2 mod q with a limited increase of the noise. Indeed, modulo q we can write:

ξ q (c 2 ), evk 0 ≡ c 2 s 2 + ξ q (c 2 ), - → a s + ξ q (c 2 ), - → e ≡ c 2 s 2 + a s + e ξ q (c 2 ), evk 1 ≡ -ξ q (c 2 ), - → a ≡ -a (11)
Now, we need to obtain the Montgomery representation of the output of this relinearization, i.e. a cryptogram like ((c 2 s 2 + a s + e )M, -M a ).

When the Montgomery representation is used, c 2 replaces c 2 in [START_REF] Fan | Somewhat Practical Fully Homomorphic Encryption[END_REF]. Hence, the relinerization key has to be modified as follows:

evkM 0 = [M 2 (P R N S,q (s 2 /M) + - → a s + - → e )] q , evkM 1 = [-M 2-→ a ] q
In the following equations, we simulate the effect of the Montgomery reduction by introducing a factor M -1 (modφ m ). Then, we can easily establish the following equalities:

ξ q ( c 2 ), evkM 0 M -1 = c 2 s 2 M + ξ q ( c 2 ), - → a sM 2 + ξ q ( c 2 ), - → e M 2 M -1 = c 2 s 2 + (a s + e )M = (c 2 s 2 + a s + e )M
Similarly, we get ξ q ( c 2 ), evkM 1 M -1 = -a M. Hence, we have obtained the Montgomery representation of the output of the relinearization step at no extra cost -both computationally and in terms of noise growth.

Impact of the Montgomery representation in BGV

For the first step of the BGV homomorphic multiplication, no scaling operation is required, and hence noise is not affected by a change in representation. Next, relinearization is applied. An analysis similar to the one in Section 4.1 can be performed, with minor adaptations to the relinearization key. Similarly, one concludes that the Montgomery reduction introduces no cost neither in terms of computation nor in noise growth. Finally, one needs to apply scaling so as to manage noise growth. We consider the ciphertext ( c 0 , c 1 ) in Montgomery representation encrypting m. Let q | q and δ i = [c i /t] q/q × t. Then the BGV-scaling function applied to c i outputs

c i = ( c i + δ i ) × q q . Lemma 3. If [c 0 +c 1 •s] q ∞ < q 2 -δ R M -1 ∞ q q t 2 (1+δ R • s ∞ ) and q = q mod t, then [( c 0 + c 1 • s)M -1 ] q = [c 0 + c 1 • s] q mod t (12)
[(

c 0 + c 1 • s)M -1 ] q ∞ q q [c 0 + c 1 • s] q ∞ + δ R M -1 ∞ t 2 (1 + δ R s ∞ )(13)
Proof. The proof is similar to the proof of lemma 4 of BGV original paper. By definition of c i , we have:

[( c 0 + c 1 • s)M -1 ] q = [c 0 + c 1 • s] q = c 0 + c 1 • s -qu.
By definition of c i , we can write:

( c 0 + c 1 • s)M -1 = q q ( c 0 + c 1 • s + δ 0 + δ 1 • s)M -1 = q q (c 0 + c 1 • s) + q q (δ 0 + δ 1 • s)M -1 = q q [c 0 + c 1 • s] q + q u + q q (δ 0 + δ 1 • s)M -1 . (14)
Moreover, since δ i ∞ qt 2q , we get the following bound

(δ 0 + δ 1 • s)M -1 ∞ δ R qt 2q M -1 ∞ (1 + δ R s ∞ ).
Thus, from the above and the hypothesis on [(

c 0 + c 1 • s)M -1 ] q ∞ = [c 0 + c 1 • s] q ∞ ,
we deduce that

( c 0 + c 1 • s)M -1 -q u ∞ < q /2
and then that

( c 0 + c 1 • s)M -1 -q u = [( c 0 + c 1 • s)M -1 ] q
Hence, from this previous equality and by bounding the norm of last member of ( 14), we obtain [START_REF] Gentry | Homomorphic Evaluation of the AES Circuit[END_REF].

Finally, we get [START_REF] Filaseta | On coverings of the integers associated with an irreducibility theorem of A. Schinzel. Number theory for the millennium[END_REF] by:

[( c 0 + c 1 • s)M -1 ] q = ( c 0 + c 1 • s)M -1 -q u = ( c 0 + c 1 • s)/M -qu mod t ( c i = c i mod t; q = q mod t) = c 0 + c 1 • s -qu mod t (def. of c) = [c 0 + c 1 • s] q mod t
From this lemma, we can see that the Montgomery representation of the ciphertext impacts the scaling by adding an extra factor δ R M -1 ∞ to the last term on the bound of the hypothesis and of (13).

Overall impact on noise growth

For both BGV and FV, the norms M ∞ , M -1 ∞ and the expansion factor δ R are involved in the noise growth due to the scaling steps performed with the Montgomery reduction. When m is a power of two, δ R is equal to n. However, for other m it can be larger than that. Let us consider

F m : Q 2n-2 [X] → Q[X]/(φ m (X )
), so that F m (a) = a mod φ m for every a ∈ Q[X] of degree lesser than or equal to 2n -2. These three parameters are given in Table 3 for the different cyclotomic polynomials considered in this paper. Assuming that distributions χ key and χ err output elements whose infinite norms are bounded by B key and B err = 6σ err , we analyse which depth can be reached in a multiplicative tree.

FV : The initial noise of a ciphertext is at most V init = B err (1 + 2δB key ) [START_REF] Fan | Somewhat Practical Fully Homomorphic Encryption[END_REF]. We recall that r ∞ = 1+ρ 2 (1 + δ R B key ) + δ R M ∞ , the output of a tree of depth L has a noise bounded by [START_REF] Joppew | Improved Security for a Ring-Based Fully Homomorphic Encryption Scheme[END_REF], Lem. 9) with:

C L 1 V + LC L-1 1 C 2 (cf.
         C 1 = 2δ R t(δ R M -1 ∞ r ∞ + 1 2 ) + δ R 2 C 2 = δ R t|q| t (r ∞ + 1) + δ R M -1 ∞ 1+δ R B k ey (1+δB k ey ) 2 + |q | t 2 + 1 +k (1 + δ R B key (1 + δ R B key )) + δ R k B err 2 ν+1
We denote by

L max = max{L ∈ N | C L 1 V + LC L-1 1
C 2 B dec } the depth allowed by the homomorphic multiplication where B dec corresponds to the decryption bound given by the full RNS version of FV [START_REF] Bajard | A Full RNS Variant of FV like Somewhat Homomorphic Encryption Schemes[END_REF].

BGV : As long as a ciphertext satisfies the condition on Lemma 3, one can perform a scaling operation and thus an homomorphic multiplication. Initially c 0 +c 1 s ∞ V init = t/2+t B err (2δ R B key +1). By assuming that after each scaling, the size of q is reduced by ω bits and let l ω,q = log ω q , then the growth of the size of c 0 + c 1 s, denoted V , after one multiplication can be expressed with:

q q (δ R V 2 + B relin ) + δ R M -1 ∞ t 1+δ R B k ey 2
where B r elin = δ R 2 ω,q ωB err B key represents the noise caused by the relinearisation step.

In Table 3 we present the maximal theoretical depths for BGV and FV with or without the use of the Montgomery reduction. For these computations we have taken parameters B key = 1, σ err = 2 √ n and a number k of 62-bits moduli to get the largest size for q ensuring at least 80-bits of security according to [START_REF] Albrecht | On the concrete hardness of Learning with Errors[END_REF].

We notice that the depths of BGV are almost unchanged with the usage of the Montgomery reduction. However for FV the depths are far smaller with a Montgomery representation. This behavior has been confirmed in practice.

m n k M ∞ M -1 ∞ δ R L BGV L M BGV L FV L M FV 4369

Mixing Optimized Barrett and Montgomery reductions

Considering the non-negligible impact of the Montgomery representation on the multiplicative depth of FV, a more robust strategy for this cryptosystem corresponds to a mixed Barrett/Montgomery approach. Alg. 2 is used during the first stage of homomorphic multiplication, with ciphertexts not exploiting a Montgomery representation. This avoids the noise growth caused by the Montgomery factor. Nonetheless, the Montgomery reduction can still be used during the relinearization stage, since we have seen that this does not cause a larger noise growth. To obtain a valid result, the relinearisation key needs to be modified, by replacing the factor M 2 of the Montgomery approach by M.

Experimental Results

The proposed methods for polynomial reduction were implemented using C++, and compiled with gcc using the optimization flag -O3. All the experimental results presented herein were measured on an i7-5960X, running at 3.0 GHz with 32 GB of main memory. No parallelism was exploited.

In Figure 1, one can find the execution timings of polynomial reduction, using NFL for power-of-two cyclotomics [START_REF] Aguilar-Melchor | Topics in Cryptology -CT-RSA 2016: The Cryptographers[END_REF]; the unoptimized and optimized Barrett reductions and the Montgomery reduction for non-power-of-two cyclotomics. In order to highlight the gain brought by our algorithms compare to generic ones we also compare with NTL's reduction using preconditioning [START_REF] Halevi | Implementing bp-obfuscation using graph-induced encoding[END_REF]. All timings were normalized based on the number of batching slots , and executed for a single moduli of 62-bits. One find the straightforward application of Barrett reduction, to be more efficient than the preconditioned methods employed in the NTL library. Moreover, speed-ups of up to 1.95 and 2.55 were achieved for the optimized Barrett and Montgomery algorithms when compared with the unoptimized Barrett reduction. The figure suggests that using power-of-two cyclotomics is not scalable with respect to the throughput. In contrast, the remaining approaches present very little variation when considering the execution timing per batching slot for the different m given in Table 3. It should be noted that using larger values of m enables FHE parameters with a larger multiplicative depth.

The aforementioned reduction methods were used to implement the homomorphic multiplication operations of the FV and BGV schemes. One can find in 3. The experimental results for NFL are omitted due to its low performance with respect to the timing per batching slot. Unlike with Figure 1, the execution time of homomorphic multiplication increases significantly with increasing m. This trend is explained by the relinearization procedure, which requires a number of NTTs that increases quadratically with log 2 q. Nevertheless, the employed reduction procedure plays a preponderant role in the efficiency of the homomorphic multiplication. Indeed, one achieves speed-ups of up to 1.37 and 1.24 when comparing the homomorphic multiplication exploiting the optimized Barrett reduction with the one exploiting the unoptimized Barrett method for the FV and BGV schemes, respectively. Since with the mixed Barrett/Montgomery approach, required by the FV scheme, one is not able to fully take advantage of the gains brought forth by the Montgomery representation, one achieves speed-ups similar to those of the optimized Barrett reduction. In contrast, for BGV, one can exploit the Montgomery representation throughout the whole procedure, leading to speed-ups of up to 1.32.

The speed-up of the proposed methods decreases as the degree n of the cyclotomic, and thus log 2 q, get larger due to the increasing complexity of the relinearization procedure. This suggests that they are most beneficial when one needs to homomorphically evaluate small circuits. Since most of practical applications of FHE [START_REF] Gilad-Bachrach | Cryptonets: Applying neural networks to encrypted data with high throughput and accuracy[END_REF][START_REF] Alhassan Khedr | SHIELD: scalable homomorphic implementation of encrypted data-classifiers[END_REF][START_REF] Naehrig | Can homomorphic encryption be practical[END_REF] have circuits with small depth, the proposed methods have a wide range of applicability.

Conclusion

In this paper, the arithmetic of non-power-of-two cyclotomics has been considered and improved. Two methods for polynomial reduction have been proposed: one based on the Barrett reduction and the other on a Montgomery representation. The optimized Barrett algorithm does not offer a better computational complexity than the Montgomery reduction. However, since it does not require changes in the representation of ciphertexts, it provides for a slower noise growth, making it more amenable to application in the FV homomorphic scheme than the Montgomery approach. In contrast, the Montgomery approach is more suitable for the BGV scheme. Moreover, since the Montgomery reduction does not rely on the properties of cyclomic polynomials, it can be used on schemes relying on other kinds of polynomials. Experimental results have shown the proposed methods to be more efficient than those employed in the NTL library. Furthermore, speed-ups of up to 1.95 and 2.55 have been obtained in an i7-5960X when comparing the optimized Barrett and Montgomery reductions with the unoptimized Barrett reduction, respectively. Finally, the polynomial reductions have been incorporated into the homomorphic multiplication procedures of FV and BGV, producing speed-ups of up to 1.37.

ω, N 4 :

 4 T (e) ← c (e) Y (e) + Q (e) Z (e) 5: R ← NTT -1 N/2 ( T (e) ) 6: return R In line 4 of Alg. 3, we require an operator which takes in as input a vector of coefficients in base B(o) ω, N . This vector defines a unique polynomial c with deg(c) < N/2. Then the operator must output the vector of coefficients of c in base B (e) ω, N . More precisely, the function BaseConv works as follows, for any c with deg(c) < N/2:

Lemma 4 .

 4 Let m be a positive integer and letR = Z[X]/(φ m (X )), with deg(φ m ) = n. If δ R denotes the expansion factor of the ring R, then δ R ≤ n • F m ∞ .

Fig. 1 .

 1 Fig. 1. Execution time per batching slot T/ [µs] for multiple reduction strategies and mth cyclotomic polynomials. The y-axis is in logarithmic scale

Fig. 2 .

 2 Fig. 2. Execution time per batching slot T/ [µs] for the homomorphic multiplication operation of FV with several reduction strategies and mth cyclotomic polynomials.

Fig. 3 .

 3 Fig. 3. Execution time per batching slot T/ [µs] for the homomorphic multiplication operation of BGV with several reduction strategies and mth cyclotomic polynomials.

Table 2 .

 2 Theoretical speed-up of Alg. 3 when compared with Alg. 1

	m	3855 4369 13107 21845 32767 65535
	n	2048 4096 8192 16384 27000 32768
	Speed-up 2.62 2.62 2.63 2.63 2.63 2.63

Table 3 .

 3 Theoretical depths with and without Montgomery reduction. Values in parenthesis are the depths observed in practice.

	4096 2	1	1	35n	1(1) 1(1) 1(4) 1(3)
	13107 8192 5	2	1	205n 4(4) 4(4) 6(17) 4(10)
	21845 16384 11	2	1	739n 10(10) 6(10) 13(40) 8(22)
	32767 27000 18	1	9	2621n 8(17) 7(17) 19(66) 12(39)
	65535 32768 22	4	1	9886n 7(21) 7(21) 22(80) 14(45)

A Proofs

A.1 Correctness of Alg. 1

Since ring F q [X] is Euclidean, we can write Euclidean division of c by Φ m , with a unique remainder r of degree at most n -1: c = c/Φ m Φ m + r. Let a, b 0 be integers. From the previous equation we get the following equivalent equations over the field of rational fractions of F q [X]:

and deg(R ) < a + b. By choosing b 0 and a α, the right term of left member of last equation above have a degree smaller than 0. In this case, we obtain an equality between the two floored polynomials. So, we take b = 0 and a = α, and we get that c Φ m is equal to the flooring of the left part of last equation, which is precisely what Alg. 1 computes. As the product X n+α /Φ m • c/X n is of degree strictly smaller than 2α + 1, the computation can be done with an NTT of size A = N 2 (2α + 1).

To finish, we notice that the result of the computation of r = cc/Φ m ×Φ m has a degree strictly smaller than n. Moreover, the polynomial c at line 5 is nothing but c/Φ m × Φ m mod X ñ -1. Indeed, the reduction modulo X ñ -1 is a consequence of the NTT based polynomial product in dimension ñ. Thus, at the end we have that c

The last equality holds precisely because the degree of c mod Φ m is strictly smaller than ñ.

A.2 Proof of Lemma 1

The first point is a direct consequence from Lemma 2 in [START_REF] Filaseta | On coverings of the integers associated with an irreducibility theorem of A. Schinzel. Number theory for the millennium[END_REF]. Since m is not a power of two, m cannot divide N. By denoting m = 2 r m with m > 1 an odd integer we have n = 2 r-1 ϕ(m ), thus 2n = 2 r ϕ(m ) and then if N divides m, N 2 (ϕ(m )) = 1 which is not possible since m 3. Therefore N and m do not divide each others and we can apply Lemma 2 from [START_REF] Filaseta | On coverings of the integers associated with an irreducibility theorem of A. Schinzel. Number theory for the millennium[END_REF].

Let α be a root of φ m in the algebraic closure of Z p . If α is also a root of X N -1 then α N = 1, since α is of order m by definition of φ m it implies that m divide N which is impossible since N is a power of two and m is not. So, φ m and X N -1 are coprime on the algebraic closure of Z p thus in Z p . The second point comes from Bezout equality in Z p and from the fact that X N -1 ≡ (X N/2 -1)(X N/2 +1) mod p.

A.3 Correctess of Alg. 4

First we notice that since c is a polynomial of degree smaller than N/2 we have (c(ω), . . . , c(ω N -1 )) = NTT N/2,ω 2 (ξ N/2,ω (c)) where ξ N/2,ω (c 0 , c 1 , . . . , c N/2-1 ) = (c 0 , c 1 ω, . . . , c N/2-1 ω N/2-1 ). Therefore the polynomial c recovered at the first line of Alg. 4 is nothing but ξ N/2,ω (c). The second line of the algorithm is just the computation of ξ N/2,ω -1 (c ) to recover c from c . Once c recovered we just need to compute and return NTT N/2,ω 2 (c) which is precisely the final step of the algorithm.

A.4 Proof of Lemma 4

Let a and b two elements of R-{0}. They naturally embed in Z

We can write ab ∞ ≤ n a ∞ b ∞ . As the product ab has degree at most 2n -2 with coefficients in Z, it belongs to Q 2n-2 [X]. Since F m is a linear map between two vector spaces of finite dimension it is continuous, then we obtain