Jun Zhang
email: junzhang@xjtu.edu.cn

Xiang Xu
email: xuxiang@stu.xjtu.edu.cn

Jose ́outeiro
email: jose.outeiro@ensam.eu

Hongguang Liu

Wanhua Zhao
email: whzhao@xjtu.edu.cn

Speed Machining of OFHC Copper Using Cellular Automata Method

Keywords: high-speed machining, cellular automaton, dynamic recrystallization, grain refinement, machining processes, modeling and simulation

Simulation of Grain Refinement

Induced by High-

Introduction

High-speed machining (HSM) has been widely used due to its advantages of low cutting forces, small thermal deformation, high material removal rates, and low cost [START_REF] Rahman | A Review on High-Speed Machining of Titanium Alloys[END_REF][START_REF] Dudzinski | A Review of Developments Towards Dry and High Speed Machining of Inconel 718 Alloy[END_REF]. High-speed machining is a highly nonlinear thermal-mechanical coupling process, in which the microstructure in the deformation zone is altered under high strain rate, complex state of stress, and relatively high temperature [START_REF] Outeiro | Experimental and Numerical Assessment of Subsurface Plastic Deformation Induced by OFHC Copper Machining[END_REF][START_REF] Courbon | Further Insight Into the Chip Formation of Ferritic-Pearlitic Steels: Microstructural Evolutions and Associated Thermo-Mechanical Loadings[END_REF][START_REF] Arrazola | Recent Advances in Modelling of Metal Machining Processes[END_REF]. Consequently, the mechanical behavior of the work material is affected, which in turn affects the quality of machined surface, the wear resistance, and fatigue life of workpiece [START_REF] Outeiro | Experimental and Numerical Assessment of Subsurface Plastic Deformation Induced by OFHC Copper Machining[END_REF][START_REF] Saoubi | A Review of Surface Integrity in Machining and Its Impact on Functional Performance and Life of Machined Products[END_REF][START_REF] Pu | Grain Refined and Basal Textured Surface Produced by Burnishing for Improved Corrosion Performance of AZ31B Mg Alloy[END_REF]. With the rapid improvement of material microscopic testing technique and experimental methods, metallurgical structure and grain size are found to be the key factors influencing the mechanical behavior of the work material in machining, thus the chip formation process and surface quality. So, the microstructure is the indispensable link to establish the relationship between machining deformation parameters and physical and mechanical properties of materials [START_REF] Arisoy | Prediction of Machining Induced Microstructure in Ti-6Al-4V Alloy Using 3-D FE-Based Simulations: Effects of Tool Micro-Geometry, Coating and Cutting Conditions[END_REF][START_REF] Denguir | A Physical-Based Constitutive Model for Surface Integrity Prediction in Machining of OFHC Copper[END_REF][START_REF] Melkote | A Physically Based Constitutive Model for Simulation of Segmented Chip Formation in Orthogonal Cutting of Commercially Pure Titanium[END_REF]. However, the mechanism of microstructural evolution has not been fully revealed at present, so the model that can accurately reflect and predict microstructural evolution is of great significance for characterizing variation of materials properties and building connection between surface integrity and machining process.

By means of advanced material testing techniques, such as electron backscatter diffraction (EBSD), X-ray diffraction, and transmission electron microscope, several researches have investigated the evolution of grain structure and grain size with different deformation conditions and machining parameters [START_REF] Saoubi | Application of the EBSD Technique for the Characterisation of Deformation Zones in Metal Cutting[END_REF][START_REF] Shekhar | Multimodal Ultrafine Grain Size Distributions From Severe Plastic Deformation at High Strain Rates[END_REF][START_REF] Huang | Study of Machining-Induced Microstructure Variations of Nanostructured/Ultrafine-Grained Copper Using XRD[END_REF]. A great number of studies have shown that the grain refinement induced by machining process has remarkable influence on the strength and hardness of materials [START_REF] Arisoy | Prediction of Machining Induced Microstructure in Ti-6Al-4V Alloy Using 3-D FE-Based Simulations: Effects of Tool Micro-Geometry, Coating and Cutting Conditions[END_REF][START_REF] Denguir | A Physical-Based Constitutive Model for Surface Integrity Prediction in Machining of OFHC Copper[END_REF][START_REF] Melkote | A Physically Based Constitutive Model for Simulation of Segmented Chip Formation in Orthogonal Cutting of Commercially Pure Titanium[END_REF]. Moreover, the grain refinement of machined surface under extrusion of tool clearance face could generate ultrafine grains even nanoscale grains, which results in the appearance of white layers [START_REF] Rancic | Microstructural Investigations of the White and Deformed Layers Close to the Turned Surface of Ti-6Al-4V[END_REF][START_REF] Xu | Grain Refinement Mechanism Under High Strain-Rate Deformation in Machined Surface During High Speed Machining Ti6Al4V[END_REF][START_REF] Nie | Plastic Strain Threshold Determination for White Layer Formation in Hard Turning of AISI 52100 Steel Using Micro-Grid Technique and Finite Element Simulations[END_REF]. The microstructure alterations, including grain refinement, deformation twining, and phase transformation, would lead to the increasing of microhardness and flow stress in machining process [START_REF] Wang | Microhardness Prediction Based on a Microstructure-Sensitive Flow Stress Model During High Speed Machining Ti-6Al-4V[END_REF]. Besides that, the tribology behavior in machining process also strongly depends on microstructure property of materials [START_REF] Courbon | The Contribution of Microstructure and Friction in Broaching Ferrite-Pearlite Steels[END_REF]. The application of microstructure-based flow stress model has significant improvement of numerical simulation accuracy in predicting mechanical process [START_REF] Cheng | Microstructure Integrated Modeling of Multiscan Laser Forming[END_REF]. In order to accurately predict microstructural evolution and study its effect on cutting process, the development of microstructure numerical model related to machining process has made great progress in recent years. Based on Estrin and Kim's dislocation density-based grain refinement material model [START_REF] Estrin | Modelling Microstructure Evolution Toward Ultrafine Crystallinity Produced by Severe Plastic Deformation[END_REF], a finite element (FE) model that can predict dislocation density and grain size distribution in deformation zone has been developed and verified to be suitable for machining simulation of many kinds of materials [START_REF] Ding | Modeling of Grain Refinement in Aluminum and Copper Subjected to Cutting[END_REF][START_REF] Ding | Dislocation Density-Based Grain Refinement Modeling of Orthogonal Cutting of Titanium[END_REF], including the oxygen-free high-conductivity (OFHC) [START_REF] Denguir | A Physical-Based Constitutive Model for Surface Integrity Prediction in Machining of OFHC Copper[END_REF]. Recent studies have identified that the grain refinement of many materials is mainly resulted from the dynamic recrystallization (DRX) behavior during machining process [START_REF] Tabei | Dynamic Recrystallization of Al Alloy 7075 in Turning[END_REF][START_REF] Behnagh | Experimental Analysis and Microstructure Modeling of Friction Stir Extrusion of Magnesium Chips[END_REF]. Based on Johnson-Mehl-Avrami-Kolmogorov DRX model, the evolution of grain size and microhardness induced by DRX can be simulated and predicted [START_REF] Pan | Prediction of Machining-Induced Phase Transformation and Grain Growth of Ti-6Al-4 V Alloy[END_REF][START_REF] Wang | Evolutions of Grain Size and Micro-Hardness During Chip Formation and Machined Surface Generation for Ti-6Al-4V in High-Speed Machining[END_REF][START_REF] Imbrogno | A Physically Based Constitutive Model for Predicting the Surface Integrity in Machining of Waspaloy[END_REF]. The development of models for grain refinement in machining process has made advances. However, the current models are based on the assumption of the homogenization of microstructural parameters. Although these models have access to present the distribution of dislocation density and grain size in macroscale, the disadvantage of these models is that they cannot reflect the grain refinement evolution process in mesoscale and show the final deformed and refined microstructure characteristics.

In studying of microstructure evolution under large plastic deformation, Goetz and Seetharaman established a dislocation density evolution model to simulate the grains nucleation and growth process induced by DRX [START_REF] Goetz | Modeling Dynamic Recrystallization Using Cellular Automata[END_REF]. Ding and Guo studied the DRX phenomenon of OFHC copper under low strain rate thermoplastic deformation and proposed a theoretical model of grain nucleation and growth and grain boundary migration to analyze the effect of temperature and initial grain size on the DRX process [START_REF] Ding | Coupled Quantitative Simulation of Microstructural Evolution and Plastic Flow During Dynamic Recrystallization[END_REF]. Based on these dislocation density and grain evolution theories, cellular automata (CA) method has been applied successfully in many researches to simulate the DRX process for various materials [START_REF] Wang | Modeling of Dynamic Recrystallization of Magnesium Alloy Using Cellular Automata Considering Initial Topology of Grains[END_REF][START_REF] Li | Simulation of Dynamic Recrystallization in AZ80 Magnesium Alloy Using Cellular Automaton[END_REF][START_REF] Hallberg | Simulation of Discontinuous Dynamic Recrystallization in Pure Cu Using a Probabilistic Cellular Automaton[END_REF]. Therefore, based on the physical mechanism of dislocation and grain evolution, the CA method is an available way to conduct simulation of microstructure evolution under thermoplastic deformation in mesoscale. Few researches have investigated the microstructure evolution in machining process by using CA method. In addition, existing CA models are generally based on the discontinuous dynamic recrystallization theory [START_REF] Huang | A Review of Dynamic Recrystallization Phenomena in Metallic Materials[END_REF], which has a lot of limitations for HSM with high strain rate and large deformation. Studies have shown that the dDRX behavior could be inhibited by high strain rate, and the cDRX begins to play a more important role gradually [START_REF] Meyers | Microstructural Evolution in Adiabatic Shear Localization in Stainless Steel[END_REF][START_REF] Liu | Effect of Microstructure Evolution on Chip Formation and Fracture During High-Speed Cutting of Single Phase Metals[END_REF]. So far, a great number of studies have shown that the grain refinement induced by machining process has remarkable influence on the strength and hardness of materials. However, there are few models developed to predict the evolution of the grain refinement of machining process in mesoscale with sufficient accuracy. Few researches investigated the microstructure evolution in cutting process by using CA method. In this work, the application of CA method in cutting leads a new approach to simulate the microstructure morphology during machining OFHC copper. In addition, considering the effect of high strain rate of HSM on microstructure evolution, continuous dynamic recrystallization model is added for the first time to conduct CA simulation and the model combining cDRX and dDRX can simulate the microstructure of deformation zone of HSM more accurately.

To investigate the effect of DRX on grain refinement of materials during HSM of OFHC copper, a typical face-centered cubic single phase material is selected to study microstructure evolution without considering the effect of phase transformation. In this paper, the CA model combing cDRX and dDRX is established to study the effect of high temperature, large deformation, and high strain rate of HSM on microstructure evolution. FE simulation was first carried out to obtain the deformation parameters of HSM, and then based on these parameters CA method was used to simulate the microstructure evolution process. Besides that, the HSM experiments and microstructure observations were conducted to verify accuracy of the CA model. Finally, the influence of deformation parameters on microstructure evolution and final grain characteristics is analyzed with the cutting speed increasing.

2 Experimentation and Finite Element Simulation of High-speed Machining 2.1 Experimental Cutting Tests. Dry cutting condition and cemented carbide inserts (model: SandvikN331.1A) were used in the machining tests. A milling cutter with a diameter of 80 mm and a helix angle of 0 deg was used by up milling to reach orthogonal high cutting speeds. Tool geometry was represented by a rake angle γ 0 of 0 deg, a clearance angle α 0 of 15 deg, and a tool edge radius r n of 10 μm. The thickness and width of workpiece were 3 mm and 20 mm, respectively. The width of the workpiece is designed to be only a quarter of the tool diameter to avoid significant changes in cutting thickness, so the cutting thickness here is assumed to be constant. The axis of the cutter is kept aligning with the left edge of workpiece. As shown in Fig. 1, the geometry of the workpiece is specially designed by the Buda's method [START_REF] Buda | New Methods in the Study of Plastic Deformation in the Cutting Zone[END_REF]. The workpiece including several holes with a dimension of 3 mm × 3 mm was adopted to realize quick-stop and obtain chip roots. The total cutting length for one feed is about 31 mm. As shown in Fig. 2, the cutting experiments were carried out on a five-axis machine center (Model: DMU 50). The milling tools were redesigned to realize orthogonal cutting. The main cutting edge of inserts was vertical to cutting plane and the rake face of inserts coincided with the axis of cutter. During the cutting process, cutting forces were measured by a piezoelectric dynamometer (Model: Kistler 9265B), and chip root samples were obtained to measure the shear angle. The thermal, mechanical, and physical properties of both workpiece and tool are listed in Table 1.

Table 2 shows the cutting parameters with the rotation speed varying from 3000 to 12000 rpm, a constant feed per tooth of 0.15 mm/z, the cutting feed varying from 450 to 1800 mm/min, and the cutting speed varying from 750 to 3000 m/min. Chips obtained at various cutting speeds are collected for further microstructure analysis by EBSD.

Finite Element Simulation.

To predict the distributions of strain, strain rate, and temperature, a 2D FE model of orthogonal cutting for plane strain condition was developed and schematically represented in Fig. 4. The constitutive material model proposed by Johnson-Cook [START_REF] Holmquist | A Computational Constitutive Model for Glass Subjected to Large Strains, High Strain Rates and High Pressures[END_REF] describes the material flow stress in function of the strain, strain rate, and temperature was used and is given by the following equation: where σ is flow stress, ɛ and ε, respectively, are equivalent plastic strain and strain rate, T is instantaneous temperature. ε0 is the reference strain rate, T r and T m are room and melted temperature, respectively. A, B, n, C, and m are material constants, which is strongly dependent on the condition of split Hopkinson pressure bar test.

σ = (A + Bε n ) 1 + Cln ε ε0 1 - T -T 0 T m -T 0 m (1)
The values of the coefficients of this model for the OFHC copper are specified in Table 3.

In order to realize the material failure and chip separation process of machining, the node separation and element deletion technique are specified. In this study, an energetic-based ductile fracture criterion is applied, in which two stages, i.e., damage initiation and damage evolution, are included.

The damage initiation is described by the variable w, which is defined as

w = Δ ε pl ε pl 0 (2)
where ε pl 0 is the equivalent plastic strain of damage initiation, Δ ε pl is equivalent plastic strain. And ε pl 0 can be calculated by the following equation, which is dependent on stress triaxiality, strain rate, and temperature

ε pl 0 = D 1 + D 2 exp D 3 σ p σ 1 + D 4 ln ˙ ε ˙ ε 0 1 + D 5 T -T r T m -T r (3)
where parameters D 1 -D 5 are the JC damage constant, which is given in Table 4.

As shown in Fig. 3, w increases continuously with plastic deformation, and when w reaches a critical value of 1, damage evolution is activated. In the stage of damage evolution, the fracture energy is obtained by integrating yield stress with plastic deformation and it can be described as

G f = ε pl f ε pl 0 Lσ y d ε pl = u pl f 0 σ y d u pl ( 4
)
where L is the characteristic length, u pl is equivalent displacement, and u pl f is equivalent displacement at failure. According to the research of Mabrouki et al. [START_REF] Mabrouki | Numerical and Experimental Study of Dry Cutting for an Aeronautic Aluminium Alloy (A2024-T351)[END_REF], the fracture energy G f can be calculated by fracture toughness K C as

G f = 1 -v 2 E K 2 C ( 5
)
The stiffness degradation of materials in condition of damage evolution can be described by variable D, and when D increase from 0 to 1, the damage evolution process completes where equivalent displacement u = L ε and equivalent displacement at failure u pl f = 2G f /σ y0 . The cutting process with several cutting speeds (v c = 750, 1500, 2250, 3000 m/min) is simulated by finite element analysis in ABAQUS/Explicit. Tool geometry and material is the same with the experiments. The mesh approach and boundary conditions of the orthogonal cutting model are shown in Fig. 4. Four-node quadrilateral continuum elements (CPE4RT) with reduced integration and hourglass control features are used to mesh the workpiece for a coupled displacement-temperature analysis. A constant uncut chip thickness of 0.15 mm and cutting width of 3 mm were used in the simulations. Three-node linear displacement and temperature elements (CPE3T) with plane strain feature are chosen to mesh the tool. Reducing mesh size gradually is considered as an effective way to balance calculation accuracy and efficiency. The mesh size of the uncut chip and sacrificial layer is determined to be 5 μm. Besides that, the mesh size of tool in tool-chip contact zone is set less than that of workpiece. The total number of elements of tool and workpiece are 48,000 and 1339. The workpiece is fixed on the bottom and left side and the tool is moving in the X-direction only, from right to left at constant cutting speed. The length of the cutting path is set to be 2.2 mm to achieve the steady state. The initial temperature of the workpiece and tool was set to 20 °C.

D = L ε u pl f = u u pl f (6)
The stick-slip friction contact model proposed by Zorev is used [START_REF] Zorev | Interrelationship Between Shear Processes Occurring Along Tool Face and on Shear Plane in Metal Cutting[END_REF], which can be represented by the following equation:

τ f = τ Y , μσ n ≥ τ Y (sticking region) μσ n , μσ n ≥ τ Y (sliding region) ( 7
)
where μ is the Coulomb friction coefficient with a value of 0.3 adopted in this study, σ n is the normal stress. The yield shear stress, τ Y , is calculated by the estimated yield stress of the material of the chip, given by τ Y = σ Y / 3 √ . To verify the accuracy of the FE machining model, the cutting tangential force and chip morphology obtained from simulation and experiments under different cutting speeds were selected to compare. When analyzing the experimental cutting force, to avoid the vibration of cutting force caused by the instability at the beginning and end of cutting process, the stable periods of cutting process were selected to calculate the average tangential forces. The comparison of simulated and experimental cutting force is listed in Table 5, from which it can be summarized that cutting forces obtained by experiment is slightly larger than the simulated value, and the error is less than 10% and within the reliable range.

The shear angle and chip compression ratio (CCR) is also taken into account to verify the accuracy of plastic deformation from simulation. The CCR was calculated by the ratio of the chip thickness to the uncut chip thickness [START_REF] Astakhov | The Assessment of Plastic Deformation in Metal Cutting[END_REF]. It can be seen from Table 6 that the shear angle and CCR obtained by this model show a well agreement with experiments.

Figure 5 show the distributions of the equivalent plastic strain, strain rate, and temperature in the deformation zone at a cutting speed of 1500 m/min. The evolution of these variables is analyzed along both chip free surface, middle area of chip, and chip surface in contact with the tool. The strain in chip surface in contact with the tool reaches about 2.0, which is higher than 1.5 in the middle area of the chip and 1.0 of chip free surface due to the influence of tool-chip friction. The maximum of strain rate distributed along the shear plane is close to 10 6 /s. In the second deformation zone, the temperature goes up significantly to about 480 K due to extrusion and friction of the cutting tool and reaches the critical temperature of DRX, which has great influence on the microstructure evolution. Generally, the severe plastic deformation has the characteristics of high strain rate and high temperature during HSM. Then, the maximum values of strain, strain rate, and temperature obtained along the chip surface in contact with the tool (path 3) are extracted and represented in Table 7. The CA simulation of microstructure evolution from Samanta et al. [START_REF] Samanta | Cellular Automaton Simulation of Microstructure Evolution for Friction Stir Blind Riveting[END_REF] showed that the simulated microstructure depended on the thermomechanical loading condition and the accurate acquisition of strain rate and temperature was the key to conduct a reliable CA simulation. So these values in Table 7 are used as input boundary conditions in the CA model for predicting DRX during HSM.

Theoretical Model for Dynamic Recrystallization

Dynamic recrystallization occurs only when the dislocation density in a deforming matrix reaches a critical level, which depends on physical deformation parameters, such as strain, strain rate, and temperature during HSM. DRX includes two types: dDRX and cDRX. Two important aspects of dDRX are nucleation and grain growth based on grain boundary migration while cDRX appears as lattice rotation and grain subdivision, as shown in Fig. 6. In order to simplify the analysis, two assumptions are proposed here.

(1) Initial dislocation density ρ ini in the matrix material is uniform and when it reaches a critical value, nucleation of dDRX occurs. For the new recrystallization grains (R-grains), the dislocation density evolves from initial dislocation density. The grains are set with an initial dislocation density of 10 9 /m 2 , representative for the annealed and undeformed state [START_REF] Hordon | Dislocation Density and Flow Stress of Copper[END_REF]. Corresponding to the grain boundary bulging nucleation mechanism, dDRX nucleate at both initial grain boundaries and R-grain boundaries. (2) For OFHC copper, dislocation cross-slip and low angle boundaries are easily observed during hot deformation, especially at high strain rate and large strain, and cDRX occurs when the critical dislocation density is achieved. Equiaxed recrystallization grains with high orientation angle are transformed from deformation substructures by lattice rotation [START_REF] Meyers | Shear Localization in Dynamic Deformation of Materials: Microstructural Evolution and Self-Organization[END_REF]. In this paper, assume that a grain is divided into three grains with 120 deg.

3.1 Dislocation Density. Dislocation density plays an important role in microstructure evolution, which is temperature and strain rate dependent. The dependence of the dislocation density ρ x,y at the site (x, y) on the shear strain ɛ is given by

dρ x,y dε = k 1 ρ x,y + k 2 ρ x,y (8)
which is the Kocks and Mecking (KM) law [START_REF] Mecking | Kinetics of Flow and Strain-Hardening[END_REF], where k 1 is a constant that represents hardening and k 2 is the softening parameter representing the recovery of dislocations. The high-temperature flow stress is related to the average dislocation density ρ by

σ = αμb ρ √ ( 9
)
where α is a dislocation-interaction term, which is 0.5 for most metals, μ is the shear modulus, and b is Burger's vector. The values of k 1 and k 2 are required from stress-strain curve. k 2 is strongly dependent on the condition of strain rate and temperature. k 1 is related to k 2 and saturated stress σ s . k 1 and k 2 can be determined by following equations [START_REF] Lee | Numerical Modeling of Dynamic Recrystallization During Nonisothermal Hot Compression by Cellular Automata and Finite Element Analysis[END_REF]: Fig. 6 Mechanisms of (a) dDRX process [START_REF] Cram | Modelling Discontinuous Dynamic Recrystallization Using a Physically Based Model for Nucleation[END_REF] and (b) cDRX process [START_REF] Mishra | Microstructural Evolution in Copper Subjected to Severe Plastic Deformation: Experiments and Analysis[END_REF]

σ = σ s [1 -exp(-k 2 ε)] 1/2 (10)
k 1 = k 2 σ s αμb 2 (11)
According to the condition of strain rate and temperature and stress-strain curve at different cutting speeds, the values of k 1 and k 2 at different cutting speed are listed in Table 8.

The accumulation of dislocations could lead to DRX. Roberts and Ahlblom [START_REF] Roberts | A Nucleation Criterion for Dynamic Recrystallization During Hot Working[END_REF] have proposed that the critical dislocation density ρ cr depends on deformation conditions and can be calculated after considering the free energy change according to the following equations:

ρ cr = 20γε 3blMτ 2 1/3 (12) τ = c 2 μb 2 ( 13
)
σl μb = K 1 ( 14
)
where γ is the grain boundary energy, ε is strain rate, l is the dislocation mean free path, M is the grain boundary mobility, τ is the dislocation line energy, c 2 is a constant of 0.5, and K 1 is a constant which is about 10 for most metals.

3.2 Nucleation and Growth of dDRX. The nucleation rate P d for dDRX is a function of both the temperature T and the strain rate ε and can be described by the following equation proposed by Kugler and Turk [START_REF] Kugler | Modeling the Dynamic Recrystallization Under Multi-Stage hot Deformation[END_REF]

P d = Cε exp - Q a RT Δt ( 15
)
where C is a constant of 1, Q a is the activation energy, R is the gas constant, and Δt is time-step. The growth driving force of an R-grain comes from the reduction of storage energy, or in other words the difference of dislocation densities between the R-grain and the matrix. It is assumed that the growth velocity V i of the ith R-grain is directly proportional to the mobility of the boundary and driving force per unit area, thus

V i = λMF i 4πr 2 i ( 16
)
where λ is a constant of 1, r i is the radius of the ith R-grain. The mobility of grain boundary, M, is given by the following equation proposed by Stüwe and Ortner [51]

M = δD ob b KT exp - Q b RT ( 17
)
where δ is the characteristic grain-boundary thickness, D ob is the boundary self-diffusion coefficient, K is Boltzmann's constant, T is the temperature, Q b is the boundary-diffusion activation energy, and R is gas constant.

If the R-grain is assumed spherical, the driving force F i can be derived from the energy change of surface and volumetric energy given by the following equation [START_REF] Qian | Cellular Automata Simulation of Microstructural Evolution During Dynamic Recrystallization of an HY-100 Steel[END_REF]:

F i = 4πr 2 τ(ρ m -ρ d ) -8πrγ ( 18
)
where γ is the grain boundary energy which is dependent on the misorientation between neighboring grains, τ is the dislocation line energy, ρ m and ρ d are the dislocation density of the matrix and the R-grain.

Modeling of cDRX.

During large deformation, a large number of dislocation cell structures are generated with low dislocation density in cell interiors and high dislocation density at cell walls. The dislocation in cell interior continually migrates to cell walls resulting from high temperature of machining, and the increasing of dislocation density at cell walls gives rise to the formation of subgrains. Furthermore, the misorientations of subgrains gradually increase due to the rotation of grains, which can generate the fully recrystallized grains. According to Estrin and Vinogradov's model [START_REF] Estrin | Extreme Grain Refinement by Severe Plastic Deformation: A Wealth of Challenging Science[END_REF], the average cell size d cr is inversely proportional to the square root of the total dislocation density ρ, and combining with Eq. ( 9), the relationship between the average cell size and flow stress can be expressed by

d cr = H 0 ρ √ = αμbH 0 σ ( 19
)
where H 0 is a constant equal to 10 for OFHC copper [START_REF] Holt | Dislocation Cell Formation in Metals[END_REF]. This equation is based on the dislocation theory, in which the pile-up of dislocation can result in increasing of local dislocation density in grain interior, grain distortion, and grain refinement. It is worth noting that the propagation of dislocation in refined grains is of certain limitations so that the saturation of cDRX behavior will cause the grain size to reach the critical value. With the increasing of dislocation density, when the relationship between grain size and flow stress cannot satisfy Eq. ( 19), which means that current grain size d i is greater than the critical grain size d cr , cDRX process is carried out to conduct the segmentation and refinement of grains. The current grain size is calculated by the following equation [START_REF] Chen | Modeling the Dynamic Recrystallization in Austenitic Stainless Steel Using Cellular Automaton Method[END_REF]:

d i = 4N i a 2 π ( 20
)
where N i is the lattice number that makes up the current grain, a is the lattice size. Considering the Hall-Petch relation between material strength and grain size [START_REF] Hall | The Deformation and Ageing of Mild Steel: II Characteristics of the Lüders Deformation[END_REF][START_REF] Petch | The Cleavage Strength of Polycrystals[END_REF], with the decreasing of grain size, the intense distribution of grain boundary with high energy induces the increasing of material strength, which makes grains harder to refine. Therefore, grain segmentations caused by cDRX was carried out with a probability (P c ), which can be calculated by following equation:

P c = ωd 2 i d 2 ini ( 21
)
where ω is constant equal to 0.5, d i and d ini are current and initial grain size, respectively. The used material parameters in the CA simulation are listed in Table 9. FORTRAN by following the CA transformation rules. The CA method is an algorithm describing the discrete spatial and/or temporal evolution in a complex physical system by applying a local deterministic or probabilistic transformation rule [START_REF] Raabe | Computational Material Science: The Simulation of Materials, Microstructures and Properties[END_REF]. For the selection of lattice size and space, smaller lattice size requires more lattice number to describe the same spatial area, which helps improve the calculation accuracy but reduce the calculation efficiency. And the lattice size should be five times smaller than the final grain size in general and the whole spatial size should be at least twice bigger than the initial grain size. Based on these boundary conditions, several tests with different lattice numbers are conducted, and finally a 2D square lattice space with a size of 300 × 300 cells are determined to have a balance between calculation accuracy and efficiency, in which size of every lattice is 1 μm × 1 μm. Every lattice site has state variables: one dislocation density variable, one orientation variable, one status variable that indicates whether it is the matrix or R-grain, one status variable that judge whether it is the grain boundary, and one color variable that represents different grains. The von Neumann's neighboring rule which considers the nearest neighbors is selected and the value of a state variable at an arbitrary site (x, y) at a time (t + Δt) can be expressed as

ψ t+Δt x,y = f (ψ t x-1,y , ψ t x,y-1 , ψ t x,y , ψ t x+1,y , ψ t x,y+1 ) (22)
where ψ t x,y represents the value of variable ψ at site (x, y) at time t. Function f is the transformation rule depending on the characteristics of the system, in which the relation indicates that the lattice state variable at next step is determined by the state variable of the lattice and of its four neighbors, as shown in Fig. 7.

In order to simulate the equiaxed growth of R-grains, the following probabilistic transformation rules are applied as

F i > 0 t 0 v i dt ≥ a (23)
When the driving force F i for the growth of the ith R-grain is positive, and the fraction of R-grain t 0 V i dt/a at time t is greater than or equal to one, the distance of grain boundary migration is greater than or equal to lattice size, where V i is the growth velocity determined by Eq. ( 16), dt is the time-step, and a is the side length of one lattice. The time-step is the shortest growth time of a cell, which is calculated by the following equation:

dt = a V max = k 2 2 a Mτk 2 1 (24)
The transformation probability is defined as P z = n z /4, where n z is the number of transformed neighboring sites and 4 is the total number of von Neumann's neighboring sites. The average grain size is obtained by dividing the total area by the number of grains. The detailed steps are shown in Fig. 8 and described as follows:

(1) First step: generation of initial grains of matrix material With the reference of the natural grain growth mode, the initial grains with random orientations are generated through seeding and grain growth process, in which the simulated annealing algorithm is adopted to reduce the energy of grain boundary in space [START_REF] He | Computer Simulation of 2D Grain Growth Using a Cellular Automata Model Based on the Lowest Energy Principle[END_REF]. The total energy is provided by the grain boundary in the CA space, and the energy of the interfacial cells depends on the number of the same state of central cell and neighboring cells. The reduction of the interfacial energy after the state transition is the driving force for the natural grain growth. The average grain size of initial grains of OFHC copper is determined by EBSD results, and the morphology of grains is set to be hexagon and equiaxed.

(2) Second step: evolution of dislocation density

The critical dislocation density, shear modulus, dislocation line tension, and grain boundary mobility were calculated based on the deformation parameters obtained from FE simulation, i.e., the strain, strain rate, and temperature. After total calculation time-step is determined, the dislocation density of each lattice will be updated in every time-step in turn according to the KM model in Sec. 3.1.

(3) Third step: grains nucleation and growth process of dDRX After scanning the whole lattice space, when the dislocation density of a lattice is found to increase to the critical value, this lattice may transform as dDRX grain nuclei with probability P d . If so, this lattice is marked as dDRX grain, and then its orientation changes randomly and its dislocation density is set to be initial value ρ ini . As for the grain growth process of dDRX, when the status variables of a lattice satisfy the transformation rules in Eq. ( 23), the lattice may transform with probability P z and marked as dDRX grain, and its orientation and dislocation density will be set as same as neighboring grain. The nucleation and growth process of dDRX is shown as Fig. 9. (4) Fourth step: grains partition process of cDRX By scannning the whole lattice space, when the size of a grain is greater than the critical grain size of cDRX, this grain may be patitioned with probability P c , and the orientation of new produced grains is set as random value. The grain partition process of cDRX is displayed in Fig. 10.

(5) Fifth step: the loop and iteration process of CA simulation Cellular automata simulation of microstructure evolution is conducted with the increasing of strain, in which time-step is Δt = ε/(N ε), where N is the total iteration step, equal to 300 in this study. During each time-step, the whole lattice space is searched and related judgement and calculation are carried out to update the lattice state at the next step. The steps (2), (3), and (4) loop in every step until the predetermined deformation is achieved. After that, the results of dislocation density, grain size, grain distribution, and DRX characteristics will be output.

Microstructure Evolution at Various Cutting Speeds.

Without the effect of phase transformation, as a typical single phase material, OFHC copper is chosen for modeling to investigate the effect of DRX process on grain refinement. Strain, strain rate, and temperature obtained from finite element simulation of orthogonal cutting are used as input parameters of the CA model. Figure 11 shows the initial grain distribution with different orientations obtained by a natural grain growth algorithm. The coarse grains are just equiaxed with an average grain size of 107 μm, determined experimentally by EBSD. The final microstructure at the cutting speed of 3000 m/min is illustrated in Fig. 12. It is obviously seen that the final grain is no longer the initial coarse crystal after DRX, and almost all of it has been transformed into fine equiaxed grains. The original grain boundaries have disappeared and a large number of new grain boundaries are formed, which reflects the effect of DRX on refining grains. At the same time, the grains formed by different evolution mechanisms can be distinguished. There is not enough time for the nuclei at grain boundaries after dDRX to grow due to the high strain rate, so that these dDRX grains keep quite fine and are just composed of a few cells. The cDRX process would lead to elongated subgrain breakage and grain boundary rotation, which shows coarse grains are divided fine grains in CA simulation. Most of the cDRX grains are equiaxed and only a few shows elongated. Note that cDRX nuclei are set at the boundaries of not only initial grains but also dDRX and cDRX grains. Above all, it could be deduced that the dominating mechanism of microstructure evolution during HSM of OFHC copper is a mixed mechanism of dDRX and cDRX.

Figure 13 illustrates microstructure evolution process with the increasing strain ɛ at different cutting speed (v c = 750, 1500, 2250, 3000 m/min). At the beginning of simulation, dislocation density is too small to cause grain evolution and grains are still coarse. When the critical dislocation density is reached, nucleation of dDRX and grain division of cDRX occur and considerable microstructure refinement is achieved. Strain increases with the deformation and the effect on grain refinement is more significant. Grains of dDRX grow a little while cDRX continues to happen. Finally, the grain has been fully refined when the given strain value is reached.

Moreover, as the cutting speed increases, the final grain size is further reduced, and the number of grain boundaries is gradually increased. At 750 m/min and 1500 m/min, extremely fine grains are formed by dDRX mechanism because the nuclei do not have enough time to grow. In addition, the grain refinement generated by sub crystalline accumulation and grain boundary rotation is starting to play an important role. Nevertheless, the nucleation rate increases with the increase of strain rate at high cutting speed (v c = 2250 and 3000 m/min), there are more dDRX grains that nucleate but could not grow. The cDRX process is gradually complete and it is shown that more and more coarse grains are divided into fine grains. In conclusion, with increasing cutting speed, the formation mechanism of fine grains are transformed from dDRX mechanism guiding at relatively low cutting speed to the mixed mechanism of dDRX and cDRX at high cutting speed.

Figure 14 illustrates average grain sizes with the increasing strain at different cutting speeds (v c = 750, 1500, 2250, 3000 m/min). During microstructure evolution process, the rate of decrease in grain size is very fast at the beginning and then slows down gradually as the increasing of the strain. A large number of subgrains and low angle boundaries are formed by dislocation accumulation. Meanwhile, subgrain size decreases rapidly until it reaches the size corresponding to maximum strain, as a result of high dislocation density under large strain and high strain rate [START_REF] Belyakov | Grain Refinement in Copper Under Large Strain Deformation[END_REF]. As shown in Fig. 14, the average grain size is larger at the same strain with higher cutting speed. Higher temperature at higher cutting speed would enhance the annihilation of dislocations and results in slower dislocation generation rate, which leads to larger grain size. The effect of temperature on grain evolution is comparatively complex. On one hand, the dislocation density decreases with the increasing temperature and dynamic recovery and grain growth may occur. On the other hand, the dislocations in the cell boundaries rearrange and neutralize and the dislocation cells translate into subgrains at a certain temperature. High angle boundaries and fine grains are formed with the increase of misorientation angle. It is also shown that the rate of grain size decreasing slows down gradually due to a balance between grain refinement and grain growth. In general, the sizes of grains in chip deformation zone have all reached micron scale.

As shown in Fig. 15, the percentage of cDRX grains increases from 7.1% to 24.7% with the cuttitng speed increases from 750 m/min to 3000 m/min. This phenomenon results from high strain rate because the time for the formation of dDRX grains by boundary migration is far longer than that of cDRX grains by boundary rotation. During cutting process, dislocations transfer from the cell interior to its wall due to thermal deformation and lots of dislocation cells with low angle boundaries translate into cDRX grains with high angle. According to finite element results, the strain rate is up to 10 6 and the deformation time is quite short. Therefore, cDRX could meet these characteristics in kinetics and play an important role in microstructure evolution during HSM.

Experimental Observation of Microstructure in Chip and Verification of Cellular Automata Model

Figure 16 shows the distributions of grain morphology, size, and misorientation at four different cutting speeds. EBSD images with over 80% resolution are obtained. It can be seen that the grain size is greatly reduced from the initial coarse grain to the fine equiaxed grain. From the EBSD image at different cutting speeds in Fig. 16, it can be seen that larger grains are surrounded by many smaller grains at cutting speed of 750 m/min and 1500 m/min, while more elongated and fine equiaxed grains are observed. Besides that, finer grains can be found at higher cutting speed. It is indicated that grain size reaches micron scale due to the deformation during HSM, which is coincident approximately with simulation results. In the cutting deformation process, grains as elongated along the shear direction and are then refined by During HSM, many low angle boundaries form in deformation and migrate or rotate inside initial coarse grains, converting to high angle boundaries. At 750 m/min and 1500 m/min, necklacelike structures composed of fine grains around a few wavy boundaries are observed, which explains the main mechanism grain refinement is dDRX. As shown in misorientation distribution of Fig. 16, more low angle boundaries and broken elongated subgrains occur due to high strain rate at 2250 and 3000 m/min, and low angle boundary with less than 15 deg accounts for more than 40%. A large number of low angle boundaries are generated means that the cDRX process becomes more intense, which is beneficial for grain refinement. And then, more high angle boundaries and finer grains are transformed. Above all, it could be deduced that grain refinement can be interpreted in terms of dDRX and cDRX at higher cutting speed.

In this study, the average grain sizes are obtained by the mean linear intercept method. As for the experimental observed EBSD results in Fig. 15, the grain size was measured by using Channel 5 EBSD analysis software. And the CA simulated microstructures in Fig. 11 were analyzed by using an open source MATLAB program "linecut". As shown in Table 10, the results of average grain size from experimental observation and CA simulation are collected and compared. It can be seen from this table that the grains in the chip surface in contact with the tool are severely refined with average size of 2-4 μm, compared with the average initial grains of 100 μm. Moreover, the relative errors between the simulated and experimental observed grain size is less than 7%, which verifies the accuracy of the CA simulation. It can be found that the experimental and simulated results show that the grain size slightly decreases with the increase of cutting speed. In this study, microstructure evolution of OFHC copper during HSM has been thoroughly investigated. The process is modeled using a CA method with DRX mechanism based on the deformation parameters obtained from FE simulation, which has been compared with the cutting experiments and EBSD analysis. It could be concluded as following.

(1) The CA model considering both dDRX and cDRX mechanism was established and the DRX process in the chip during HSM of OFHC copper was simulated. The CA simulation results show well agreement with the experimental observation results, which verifies the feasibility of using CA method to conduct simulation of microstructure evolution in HSM of OFHC copper. (2) From the CA simulation results, it is identified that the success of CA simulation during HSM of OFHC copper depends on the supplement of cDRX model and deformation parameters (strain, strain rate, and temperature). Due to the addition of the cDRX model, the effect of high strain rate and large deformation on DRX process will be adequately considered. The high strain rate promotes an increasing dislocation density and grain refinement while grain recovery and growth under high temperature will occur, of which this CA model considering both dDRX and cDRX mechanism gives us a further understanding of the refinement process.

(3) The grain refinement mechanism during HSM of OFHC copper was analyzed by comparing the CA simulated and EBSD measured grain sizes. With increasing cutting speed, under the effect of higher strain rate and larger deformation, the dDRX process based on the grain boundary migration theory is inhibited and dDRX grain nuclei does not have enough time to grow so that the mechanism of DRX transforms with the enhance of the cDRX process. As a result, the size of equiaxed grains decreases and more low angle boundaries and deformed substructures are generated, which is thought to result in the evolution of microhardness.

Fig. 1 Fig. 2

 12 Fig.1Schematic of Buda's experimental method[START_REF] Buda | New Methods in the Study of Plastic Deformation in the Cutting Zone[END_REF]

Fig. 3

 3 Fig. 3 Stress-strain curve of material damage degradation

Fig. 5

 5 Fig. 5 Distributions of (a) strain, (b) temperature, and (c) strain rate along both chip free surface (path 1), middle area of chip (path 2), and chip surface in contact with the tool (path 3) for a cutting speed equal to 1500 m/min

Fig. 7

 7 Fig. 7 Neighboring relationship of lattices in the cell space

Fig. 8

 8 Fig. 8 Flowchart of the loop and iteration process in CA simulation

Fig. 11

 11 Fig. 11 The initial grain distribution of CA model

Fig. 13

 13 Fig. 13 Simulated microstructure evolution with strain ɛ at different cutting speeds: (a) v c = 750, (b) v c = 1500, (c) v c = 2250, and (d) v c = 3000 m/min

Fig. 14

 14 Fig. 14 Evolution of grain sizes with increasing strains at different cutting speeds

Table 1

 1 Material properties of OFHC copper and cutting tool

	Material	Density (kg/m 3 )	Elastic modulus (GPa)	Poisson's ratio	Specific heat (J/(kg K))	Thermal conductivity (W/(m K))	Thermal expansion (K -1 )
	Tool OFHC	11,900 8960	534 124	0.22 0.34	400 380	50 401	-1.65 × 10 -5

Table 2

 2 Cutting parameters of HSM tests

	Rotation	Feed (f z )	z (number of	Cutting feed (V f )	v c
	speed (rpm)	(mm/z)	teeth)	(mm/min)	(m/min)
	3000			450	750
	6000	0.15	1	900	1500
	9000			1350	2250
	12,000			1800	3000

Table 3 JC

 3

parameters of the OFHC copper used in the simulations

[START_REF] Holmquist | A Computational Constitutive Model for Glass Subjected to Large Strains, High Strain Rates and High Pressures[END_REF]

Table 4

 4

	D 1	D 2	D 3	D 4	D 5
	0.54	4.89	-3.03	0.014	1.12

JC damage constant of OFHC copper

[START_REF] Johnson | Fracture Characteristics of Three Metals Subjected to Various Strains, Strain Rates, Temperatures and Pressures[END_REF]

Table 6

 6 Shear angle and CCR obtained from experiments and simulations

	Table 5 Comparison of simulated and experimental cutting
	force
		Simulated force F t	Experimental force F t
	Cutting speed (m/	(N)		(N)		Error
	min)	Average		Average		(%)
	V = 750	612.4		634.4		3.5
	V = 1500	639.0		663.3		3.7
	V = 2250	723.0		714.5		1.2
	V = 3000	817.9		748.7		9.2
		Average		Average simulation
		experimental results	results
	Cutting speed	Shear angle		Shear angle		Error
	(m/min)	(deg)	CCR	(deg)	CCR	(%)
	V = 750	13.6	4.13	15.1	3.71	10.2
	V = 1500	15.8	3.53	16.9	3.29	6.8
	V = 2250	17.7	3.1	18.3	3.02	2.6
	V = 3000	18.3	3	19.2	2.87	4.3

Table 7

 7 Maximum values of strain, strain rate, and temperature along path 3 in Fig.5

	v c (m/min)	Strain	Temperature (K)	Strain rate (1/s)
	750	1.6	435	6 × 10 4
	1500	1.9	478	1.2 × 10 5
	2250	2.2	486	3.6 × 10 5
	3000	2.8	512	7.8 × 10 5

Table 8

 8 The values of k 1 and k 2 at different cutting speeds

	Cutting speed (m/min)	750	1500	2250	3000
	k 1	9.21 × 10 8 8.87 × 10 8 8.54 × 10 8 7.92 × 10 8
	k 2	4.62	4.55	4.44	4.22

Table 9

 9 Material parameters used in CA simulation Detailed Description of Cellular Automata Model for Dynamic Recrystallization. Based on the above theoretical model of DRX, the routine has been developed in

	Variable	Q a (kJ/mol)	Q b (kJ/mol)	μ (MPa)	b (m)	γ (J/m 2 )	R (J/(mol K))	δD ob (m 3 /s)	K (J/K)	ρ ini (1/m 2 )
	Value	261	104	4.2 × 10 4	2.56 × 10 -10	0.5	8.314	5 × 10 -15	1.38 × 10 -23	1 × 10 9

Table 10

 10 Comparison between experimental and simulated average grain sizes in chips of OFHC at different cutting speeds

	v c	CA simulated	Experimental observed	Relative
	(m/min)	grain size (μm)	grain size (μm)	error (%)
	750	3.48	3.32	4.8
	1500	3.12	3.15	0.6
	2250	2.88	2.92	1.4
	3000	2.52	2.71	7.0

Acknowledgment

This work was supported by the National Science Foundation of China (51675417), Shaanxi Science Foundation for Distinguished Young Scholars (2019JC08), and China Scholarship Council (No. 201906280422). We also thank Dr. Zijun Ren at Instrument Analysis Center of Xi'an Jiaotong University for his assistance with EBSD test. Acknowledgements.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper. Data provided by a third party listed in


Eidan Tzdaka
email: eidantz@gmail.com

Camille Benaroch
email: camille.benaroch@inria.fr

Camille Jeunet
email: camille.jeunet@univ-tlse2.fr

Fabien Lotte
email: fabien.lotte@inria.fr

Assessing The Relevance Of Neurophysiological Patterns To Predict Motor Imagery-based BCI Users' Performance

Keywords: Brain-Computer Interfaces (BCI), electroencephalography (EEG), neurophysiological predictor

Motor Imagery-based Brain-Computer Interfaces (MI-BCI) allow users to control a computer for various applications using their brain activity alone, which is usually recorded by an electroencephalogram (EEG). Although BCI applications are numerous, their use outside laboratories is still scarce due to their poor accuracy. Some users cannot use BCIs, a phenomenon sometimes called "BCI illiteracy", which impacts around 10% to 30% of BCI users, who cannot produce discriminable EEG patterns. By performing neurophysiological analyses, and notably by identifying neurophysiological predictors of BCI performance, we may understand this phenomenon and its causes better. In turn, this may also help us to better understand and thus possibly improve, BCI user training. Therefore, this paper presents statistical models dedicated to the prediction of MI-BCI user performance, based on neurophysiological users' features extracted from a two minute EEG recording of a "relax with eyes open" condition. We consider data from 56 subjects that were recorded in a 'relax with eyes open' condition before performing a MI-BCI experiment. We used machine learning regression algorithm with leave-onesubject-out cross-validation to build our model of prediction. We also computed different correlations between those features (neurophysiological predictors) and users' MI-BCI performances. Our results suggest such models could predict user performances significantly better than chance (p ≤ 0.01) but with a relatively high mean absolute error of 12.43%. We also found significant correlations between a few of our features and the performance, including the previously explored µ µ µ-band predictor, as well as a new one proposed here: the µ µ µ-peak location variability. These results are thus encouraging to better understand and predict BCI illiteracy. However, they also require further improvements in order to obtain more reliable predictions.

I. INTRODUCTION

Motor Imagery-based Brain-Computer Interfaces (MI-BCIs) allow users to control computers by producing brain activity patterns that are usually measured with Electroencephalography (EEG) [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF]. MI leads to specific rhythmic macroscopic EEG oscillations (mostly in the µ and β bands) in sensorimotor cortices, which can be detected using a BCI [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF]. Thus, by discriminating MI patterns, the recognized patterns can be associated with a specific command for the computer [START_REF] Pfurtscheller | EEGbased discrimination between imagination of right and left hand movement[END_REF], e.g., imagining feet movements to make a cursor move towards the bottom of the screen. Although BCI applications are numerous, e.g., for assistive technologies or entertainment, [START_REF] Kübler | Brain-computer communication: Unlocking the locked in[END_REF] [START_REF] Allison | Towards Practical Brain-Computer Interfaces[END_REF], their use outside laboratories is rather scarce due to their poor accuracy. On average, one out of four or five user commands is misinterpreted by BCI systems based on 2 MI classes [START_REF] Guger | How many people are able to operate an EEG-based brain-computer interface (BCI)?[END_REF]. The BCI community tried to overcome this problem in different ways such as making users train longer a given MI task, to become more BCI proficient. The larger amount of data collected this way can also be used to refine BCI classifiers. There has also been attempts to change the instructions given or to improve EEG signal processing methods and hardware (e.g., electrodes) [START_REF] Allison | Could anyone use a BCI?[END_REF]. However, we also need to consider the user, who is one of the main components in the BCI loop, and may not be able to produce reliable EEG patterns, which is sometimes called "BCI illiteracy/deficiency1 ". Producing distinct EEG patterns is a skill that must be mastered by the user [START_REF] Neuper | Motor imagery and action observation: Modulation of sensorimotor brain rhythms during mental control of a brain-computer interface[END_REF], [START_REF] Jeunet | Advances in user-training for mental-imagery-based BCI control: Psychological and cognitive factors and their neural correlates[END_REF]. Therefore, user training is also essential and should be studied and improved so that users can become proficient at BCI control. Indeed, if a naive BCI user were unable to perform the desired MI commands due to a lack of ability to produce discriminable EEG patterns, regardless of the signal processing used, the BCI system would struggle and would issue random commands. Insufficient user's attention to the imagery tasks, or frustration due to incorrect feedback may also cause poor performances [START_REF] Guger | How many people are able to operate an EEG-based brain-computer interface (BCI)?[END_REF]. While the approaches presented above have only been partly successful, they could make BCIs work for some previously "illiterate" users [START_REF] Nijholt | Brain-computer interfacing for intelligent systems[END_REF]. However, some users still remain unable to use MI-BCI systems or have poor control.

There was thus some research work dedicated to finding predictors of BCI performance, in order to quickly identify whether a BCI user is likely to be able to control a given BCI [START_REF] Jeunet | Advances in user-training for mental-imagery-based BCI control: Psychological and cognitive factors and their neural correlates[END_REF], [START_REF] Blankertz | Neurophysiological predictor of SMR-based BCI performance[END_REF], [START_REF] Ahn | High theta and low alpha powers may be indicative of BCI-illiteracy in motor imagery[END_REF]. Those studies showed that users' gender [START_REF] Randolph | Not all created equal: Individual-technology fit of braincomputer interfaces[END_REF], daily hand movement [START_REF] Randolph | Individual characteristics and their effect on predicting mu rhythm modulation[END_REF], experience with state-of-theart technology [START_REF] Burde | Is the locus of control of reinforcement a predictor of brain-computer interface performance?[END_REF] and good spatial abilities [START_REF] Jeunet | Advances in user-training for mental-imagery-based BCI control: Psychological and cognitive factors and their neural correlates[END_REF] are all related to user's MI-BCI performance.

These studies have explored users' traits as predictors, but not their neurophysiological correlates, nor other neurophysiological predictors. A different study [START_REF] Blankertz | Neurophysiological predictor of SMR-based BCI performance[END_REF] uncovered the µband predictor, which can be determined from a two-minute recording of a 'relax with eyes open' condition. This predictor had a positive correlation with users' MI-BCI performance. In [START_REF] Ahn | High theta and low alpha powers may be indicative of BCI-illiteracy in motor imagery[END_REF], the task of relaxing with open eyes was used again. They found out that high theta and low alpha EEG patterns could be observed for users with difficulty in using the BCI. Those works identified predictors correlated to BCI users performances. However, how well can they really predict the performances of new unseen users? Can we combine multiple predictors together to improve this prediction? We aim to answer these questions in this paper. In particular, its main objective is to create a computational model that could predict the performances of a BCI user based on multiple users' neurophysiological features. Thus, we first needed to identify predictors that correlate with MI-BCI performances. Our database contains 56 healthy MI-BCI naive participants who were asked to 'relax with eyes open' as a baseline for 3-minute followed by six 7-minute runs during which they had to learn to perform two MI-tasks, i.e., left and right-hand movement imagination. We used an elastic net regularized regression to predict between-subject MI-BCI performance variability. We extracted predictors based on the EEG power spectral density to build those models. We notably re-used the µ-band predictor combined with other new spectral predictors proposed here (e.g., the width of the peak of the µ-band) to create multivariate models to predict MI-BCI performances.

This paper is organized as follows: Section II describes the data set that we used and its paradigm, together with the predictor features and the approach used to build models and the elastic net regression. Then, Section III is dedicated to the results while Sections IV and V discuss the results and conclude the paper, respectively.

II. MATERIALS AND METHODS

A. Data set

We used the data of fifty-six subjects (N=56) out of the fiftynine subjects from the experiment in [START_REF] Roc | Would motor-imagery based BCI user training benefit from more women experimenters?[END_REF] (29 women; age 19-59; X =29; SD=9.318). We discarded three subjects for which the µ-predictor (see the µ-band predictor section) could not be computed due to noisy data: the optimization process failed to converge (note that we used the original authors' code). The main goal of this previous experiment was to assess the impact of the interaction between experimenters' and participants' gender on MI-BCI performances. Each subject participated in one MI-BCI session of 2 hours, comprising a three-minute recording of a 'relax with eyes open' condition followed by six 7-minute blocks (a.k.a runs) of trials in which the subject had to execute MI-tasks, i.e., imagine right or left-hand movements (1 of the 56 subjects did only four runs due to time limits).

Every subject performed 40 trials in a block, 20 per MI-task that appeared in random order. Each trial lasted for 8s. At t=0s a cross was displayed on the screen. After two seconds the subject was provided with an acoustic signal which indicated the appearance of a red arrow, that appeared one second later. The red arrow direction indicated the type of MI-task to perform, e.g., a left pointing arrow to indicate the imagination of a movement of the left hand. From t=4.250s a blue bar was displayed as feedback, representing the classifier output: its length changed according to the classifier output and it appeared only if the instruction matched the recognized task (positive feedback only). The feedback lasted 3.75s.Afterwards the screen turns black and the user had a few seconds to rest before the next trial started.

B. EEG Recordings and Signal Processing

To record the EEG signals, 27 active scalp electrodes, referenced to the left ear, were used (Fz, FCz, Cz, CPz, Pz, C1, C3, C5, C2, C4, C6, F4, FC2, FC4, FC6, CP2, CP4, CP6, P4, F3, FC1, FC3, FC5, CP1, CP3, CP5, P3, 10-20 system). EEG signals were measured using a g.USBAmp (g.tec, Austria), sampled at 256 Hz and processed online using OpenviBE 2.1.0. In order to discriminate MI tasks from EEG data, participant-specific spectral and spatial filters have been used based on the two first runs (calibration runs). First, a participant-specific discriminant frequency band was selected using the heuristic algorithm proposed in [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF]. It selects the frequency band whose power in the sensorimotor channels maximally correlates with the class labels. In this study, channels C3 and C4 were used after spatial filtering with a Laplacian filter. The discriminant frequency band was selected within the interval from 5 Hz to 35 Hz, with 0.5Hz large bins. Then EEG was filtered in that selected band with a Butterworth filter of order 5. To build the classifier, the Common Spatial Pattern (CSP) algorithm [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF] was first used to optimize 3 pairs of spatial filters. CSP maximizes the signal band power difference between the two MI conditions. The log-transformed band power of the CSP filtered signals was used as input of a Linear Discriminant Analysis classifier.

C. µ-band predictor

As predictor, we first computed the µ-band predictor proposed in [START_REF] Blankertz | Neurophysiological predictor of SMR-based BCI performance[END_REF]. We first used the raw EEG data from the three minutes recordings in a 'relax with eyes open' condition, and used band-pass filtering between 2Hz to 40Hz to reduce physiological noise and power line noise that might affect the EEG recording. We then created two Laplacian channels around C3 and C4, from 5 monopolar channels for each. From the three minutes of recording, we only used the middle two minutes to avoid user movement at the beginning and end of the recording. From these data, we computed the power spectral density (PSD) between 2Hz to 30 Hz by using Welch's method with a Fast Fourier transform (FFT) size of 4096. The PSD were computed for each of the two Laplacian channels.

As proposed in [START_REF] Blankertz | Neurophysiological predictor of SMR-based BCI performance[END_REF], we modeled the PSD curve by using the sum of two functions:

g(f ; λ, µ, σ, k) = g 1 (f ; λ, k) + g 2 (f ; µ, σ, k) (1) g 1 (f ; λ, k) = k 1 + k 2 f λ (2) g 2 (f ; µ, σ, k) = k 3 ϕ(f ; µ 1 , σ 1 ) + k 4 ϕ(f ; µ 2 , σ 2 ) (3)
where

k = (k 1 , k 2 , k 3 , k 4 ) ∈ R 4 , λ ∈ R
max f ∈f g 2 (f; µ, σ, k)-g 1 (f; λ, k).
Here, we force f to be in the range of the µ band frequencies, which means that we search for the maximum difference between g 2 to g 1 only between 5Hz to 15Hz.

D. Variable to predict

In the experiment, MI-BCI performance was assessed in terms of mean classification accuracy, i.e., the mean accuracy measured over all the windows of the feedback periods of the runs 3 to 6 (online feedback runs). We tried to predict each subject performance (i.e., their classification accuracy) using an elastic net model with features extracted from the 2-minute baseline of that subject.

E. Predictor Variables

In our study, we extracted features from the PSD curve of each subject during the recording of the baseline. These features describe the PSD properties, and are used as input of the model to predict the performance. As described in Figure 1, the features we extracted are the following:

• µ µ µ-band predictor that was computed by: max f ∈f g 2 (f; µ, σ, k)-g 1 (f; λ, k), where f is in the µ-band frequency range (5 to 15 Hz). • β β β-band predictor computed similarly to the µ-band predictor, as max f ∈f g 2 (f; µ, σ, k) -g 1 (f; λ, k), where f is in the β-band frequency range (15 to 25 Hz). • Log transformation of µ µ µ-band predictor and log transformation of β β β-band predictor to decrease the variability of the data and make it more normally distributed. • Frequency of the µ µ µ-band peak and frequency of the β β β-band peak. • Width of the µ µ µ-band peak and width of the β β β-band peak first calculated by evaluating the height h eval with h eval = h P eak -P where h P eak is the real height of the peak and P is the peak prominence, [START_REF] Virtanen | SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF]. Then an horizontal line is drawn to evaluate the height of both sides of the peak bases. Starting at the peak current vertical position until the lines either intersect a slope, Fig. 1. In the plot, we can see most of the features that were extracted from the signal. The solid black signal is the original power spectral density of the user, the dashed blue signal is the estimated floor noise (g 1 (f ; λ, k)) while the dotted red signal is the fitted values g(f ; λ, µ, σ, k). The vertical green line represent the base of the peaks (and the width of the peak for the µ-band peak). The two horizontal lines describe the predictors, the black one is for µ-band predictor and the yellow one is for the β-band predictor. The blue 'X' represents the frequency of the peak and the number of peaks in the signal i.e., 2 in this example.

the signal border or crosses the vertical position of the peak base. The width is the distance between the two chosen points. • Variance of the µ µ µ-band amplitudes across the frequencies comprised in the peak width. • Variance of the µ µ µ-band frequencies comprised in the peak width. • Number of peaks, i.e., the amount of distinct peaks in the signal which is either 0, 1 or 2. Furthermore, we also investigated how the µ peak varies across time. Indeed, its variability may explain differences in mean peak amplitude or could be related to BCI performances. For that, we make use of the whole three minutes baseline recording, and estimated the variance of the peak location, width, and amplitude across all 10s epochs with 3s overlap from that baseline. This led to three additional features:

• Variance of the frequency of the µ µ µ-band peak across time: we calculated the frequency of the µ-band peak for each epoch and calculated the variance of those values. • Variance of the width of the µ µ µ-band peak across time:

we calculated the width of the µ-band peak for each epoch and calculated the variance of those values. • Variance of the amplitude of the µ µ µ-band peak across time: we calculated the amplitude of the µ-band peak for each epoch and calculated the variance of those values. If no significant peak was found in the given epoch, all three values that were mentioned above were set to 0 for that epoch, in order to indicate the lack of peak in the µ-band.

These features were extracted from both Laplacian channels C3 and C4, leading to a total of 22 features per subject. They were then used by our model to predict MI-BCI performances. It should be noted that except the µ-band predictor, all other features were explored for the first time in this paper.

F. Correlation analyses

We performed Pearson correlation analyses for each feature, to estimate how well they were correlated to performances. As we performed numerous tests, we performed a correction for multiple testing using the Benjamini-Hochberg (BH) procedure [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF] to control the False Discovery Rate (FDR).

G. Analyses

We first normalized each feature by using Z-score. Then we built our Elastic Net model and compared its predictions to that of a random model to investigate its reliability.

1) Elastic Net: We used Elastic Net regression in order to get an interpretable model that was able to predict MI-BCI users' performance. Elastic Net regression uses l1 and l2-norm regularization with two penalty parameters λ 1 and λ 2 [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. This combination allows us to create a sparse model with some weights that are non-zero like Lasso while still maintaining the regularization properties of Ridge regression [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. Elastic Net is also more robust and stable than Lasso [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] since with Lasso variable selection can be too dependent on data and thus unstable. In addition, Elastic Net can handle correlated variables better than Lasso. Indeed, with such variable, Lasso usually randomly chooses one and ignores the others, i.e., one of the coefficients of the correlated features will be nonzero while the others will be zeros. In contrast, Elastic Net will pick both of them with similar coefficient values. To set up the model we need: y ∈ R n which here is the MI-BCI performance vector, a matrix X ∈ R n×p with p features (here the neurophysiological predictors) for n subjects and a coefficient vector β ∈ R n which is the regression weight. With Elastic Net, the regression weights are estimated by:

βe-net = argmax β∈R p y -Xβ 2 2 + λ 1 β 1 + λ 2 β 2 2 (4)
where,

u 2 2 = n i=1 u i for u ∈ R n , β 1 = p j=1 |β j | and β 2 2 = p j=1 β 2 j .
In order to determine the ratio between λ 1 and λ 2 , we performed inner cross-validation (leave-onesubject-out cross-validation) on the N-2 subjects (our 'inner' training set) to estimate the optimal ratio, i.e., the one that minimizes the mean absolute error on the training set. Then we used this optimal ratio to build the model for the N-1 training subjects. The left out subject is the test set. This process was repeated by using separately each subject as the test set.

2) Random model: The random model operates as an empirical chance level for our data. By rearranging the training set BCI performances randomly and keeping the features the same, we broke the before-known relationship between features and performance. We kept the same Elastic Net parameters (λ 1 and λ 2 ) found in the inner cross-validation. Moreover, we kept the same features value and performance score for the test set that was left out to compare the actual model to the random one. We ran the Elastic Net model for the same test subject with the randomly shuffled training set. We did this process 10,000 times and stored the mean absolute error to obtain the distribution of the chance prediction performances. We then sorted the values in descending order and computed the 95 th and 99 th percentiles, which represent the chance level for the mean absolute error for p = 0.05 and p = 0.01 respectively.

III. RESULTS

A. Correlation analyses

We obtained a correlation of r = 0.426 (p < 0.01, p adj. < 0.012 after BH procedure) (Figure 2) for the µ-band predictor. We also found a correlation of r = 0.27 (p < 0.05) between the mean of the µ peak width over both channels and the performance. In other words, the larger the width of the peak the better the performance. However, after FDR correction, this correlation was not significant anymore (p adj. < 0.11). None of the β-band predictors were significantly correlated to performances though. We did find a correlation between the number of peaks in both channels and the performance (r = 0.351, p < 0.01, p adj. < 0.032). We also found a correlation between the mean (over channels) of the variance of the µband peak frequency over time (Figure 3) and the performance (r = -0.477, p < .0001, p adj. < .0001). This correlation seems mostly driven by the number of epochs with no peak, which is also negatively correlated to the performance (r = -0.42, p < 0.001, p adj. < 0.003). Table I summarizes the results of all correlation analyses. Each feature is the mean of the values from both EEG channels.

B. Model

After using Elastic Net regression and leave-one-subjectout cross-validation to create a model for each test user, we computed the mean absolute error across all the models. The results can be found in Figure 4. This mean error was 12.43%. We compared it to the mean absolute error that we got from the 10,000 random models. The estimated chance level (from the random models) was 12.85% and 13.24% for p = 0.01 and p = 0.05 respectively. Therefore, on average, our prediction models were better than chance at p < 0.01.

IV. DISCUSSION

A. µ-band predictor

From our results, we successfully reconfirmed the value of the µ-band predictor [START_REF] Blankertz | Neurophysiological predictor of SMR-based BCI performance[END_REF] on our database, thus further emphasizing its relevance and usefulness. Since MI-BCI operates by recognizing event-related desynchronization (ERD) [START_REF] Pfurtscheller | Event-related eeg/meg synchronization and desynchronization: basic principles[END_REF] performed voluntarily by the user, the higher the µ amplitude at rest, the more noticeable its decrease (i.e., an ERD) during the MI task. The µ-band predictor represents the µ amplitude and thus, the bigger the µ-band predictor, the better the chance that the user will be successfully gain control of the MI-BCI system. In blue the real performance of each subject and in orange, the predicted performance of each subject generated using our model.

B. Correlations for the new neurophysiological features

In this paper, we also proposed and studied new neurophysiological features, some of which proving to be new significant predictors of performance. Among the significant correlations that were found, the mean width of the µ peak over both C3 and C4 channels was correlated to the users' performance. However, after FDR corrected its significance dropped under significance level which might suggest that it might not be a predictor of performance, or at least not a strong one. We did not to find significant β-band predictors. However, we did find a significant correlation between the number of peaks in both channels with the performance. In other words, our results suggest that MI-BCI users who have two distinct peaks (i.e., one µ and one β peak) perform better than those who only have one peak. The β-band peak is known to be suppressed during MI (ERD) [START_REF] Pfurtscheller | Event-related eeg/meg synchronization and desynchronization: basic principles[END_REF]. As MI-BCIs are based on these changes in the spectral power of the EEG signal in both µ and β, that might explain the correlation between the number of peaks and performances. Additionally, outside of BCI, pre-training β-band ERD in the sensorimotor cortex significantly accounted for motor performance [START_REF] Espenhahn | Cortical beta oscillations are associated with motor performance following visuomotor learning[END_REF]. Moreover, pre-movement resting β-band power in the sensorimotor cor-tex was also significantly enhanced after training a motor task and performing it successfully [START_REF] Espenhahn | Cortical beta oscillations are associated with motor performance following visuomotor learning[END_REF]. Both observations might thus also explain why having two peaks, including a β one, at rest could also explain motor BCI performances.

Finally, we found that the more stable the µ-band peak (in terms of frequency location), the better the chance for the subject to perform well. A possible interpretation for that phenomenon could be that the more variable the µ location, the more its amplitude would be averaged out across time, thus making the detection of an ERD more difficult and lowering BCI performances. Alternatively, this µ peak variability might make it more difficult for machine learning algorithms to identify a clear and stable EEG pattern to use for classification. This would need to be studied in more details in the future. Altogether, these observed correlations should help us to understand and predict MI-BCI performance.

C. Elastic net model

Altogether, the Elastic Net, although significantly better than chance, had a relatively high mean absolute error. With more data we might obtain better predictions. Moreover, with more data, we might use a different cross-validation method which might lead to more accurate predictions. We do plan to add more features to our model, both neurophysiological predictors, e.g., before or during trials, and non-neurophysiological ones, e.g., traits [START_REF] Benaroch | Are users' traits informative enough to predict/explain their mental-imagery based bci performances?[END_REF]. The feedback paradigm that was used in the BCI experiment, although a standard one (which is why we studied it), is also known to be non-optimal [START_REF] Jeunet | Why standard brain-computer interface (BCI) training protocols should be changed: an experimental study[END_REF]. This may explain why several of our subjects had below chance level performance. Perhaps due to the later, the best performers were underestimated by the models.

V. CONCLUSIONS

In this study, we used data of 56 subjects to predict their performances based on features we extracted from their EEG signal at rest. The objective was to create a machine learning model to predict the user BCI performance, and therefore better understand, and possibly improve BCI user training protocols. In order to achieve this goal, we used an Elastic Net that promotes sparse solutions. Such model obtained better than chance prediction performance. We hypothesize that by combining this model with other factors such as psychometric and personality questionnaires data, we might be able to create a more robust and precise model of prediction for MI-BCI performances. We also identified three different neurophysiological predictors that are significantly correlated to the users' BCI performance. Altogether, this study contributed to increase our understanding of BCI illiteracy/deficiency and paves the way for more advanced computational modeling of MI-BCI users performances. On the longer-term, such models should help to better understand and thus optimize BCI designs.

  and ϕ(•; m, s) indicates the probability density function of a normal distribution with mean m and standard deviation s. Function g 1 models the noise spectrum and g 2 models the peaks around the µ (8-13Hz) and β (13-30Hz) bands. Therefore µ 1 and µ 2 represent the µ and β peaks locations, σ 1 and σ 2 their scale and k 1 and k 2 their amplitude. The nine parameters are estimated by minimizing the L 2 -norm of the difference vector P SD(f)g(f; λ, µ, σ, k), where f is the vector of all available frequency values for the PSD. The µ-band predictor is defined as

Fig. 2 .

 2 Fig. 2. Correlation of the µ-band predictor with MI-BCI performance, each dot corresponds to one participant. The linear line is the least-squares polynomial fit of the µ-band predictor and the corresponding MI-BCI performance.

Fig. 3 .

 3 Fig. 3. Correlation of the mean of the variance of the peak frequency over time on both channels with MI-BCI performance, each dot corresponds to one participant. The linear line is the least-squares polynomial fit of the mean of the variance of the peak frequency over time on both channels and the corresponding MI-BCI performance.

Fig. 4 .

 4 Fig. 4. Results of the Elastic net predictions compared to the real performance.In blue the real performance of each subject and in orange, the predicted performance of each subject generated using our model.

There is valid criticism about the broad use of the term "BCI illiteracy". Indeed, this concept seems to rely on an inaccurate presumption that users and their traits are the sole responsible for poor BCI performances[START_REF] Thompson | Critiquing the concept of BCI illiteracy[END_REF].

ACKNOWLEDGMENTS

This work was supported by the European Research Council with project BrainConquest (grant ERC-2016-STG-714567).